turf icon indicating copy to clipboard operation
turf copied to clipboard

Changed line-offset algorithm to account for spherical coordinates

Open PereUbu7 opened this issue 5 years ago • 6 comments

Since the current algorithm is only correct at the equator i'm suggesting the following change:

  • Use spherical coordinate and dot product between the two points to get 'real' distance
  • Transform a local cartesian coordiante system from the midpoint that accounts for the fact that longitudes are squashed closer to the poles.

PereUbu7 avatar Sep 24 '20 07:09 PereUbu7

Hi @PereUbu7. Understand it's been a while since you proposed this. Are you still interested in getting this merged with some assistance from us?

smallsaucepan avatar Nov 10 '23 04:11 smallsaucepan

@smallsaucepan Hi there! I'd love to. Just got to refresh and recapitulate what this was all about :)

PereUbu7 avatar Nov 10 '23 05:11 PereUbu7

@smallsaucepan Hi again. I've written some unit tests for the process segment method in vanilla js. Could you guys help putting this in the test framework you're actually using and maybe remove any old tests if there are any? :)

` "use strict"

const IsWithin = (value1, value2, variance) => { return value1 > (value2 - variance) && value1 < (value2 + variance); }

const MeterInKm = 0.001; const KmInDegrees = 0.00899;

const TestAbsoluteDistance = (point1, point2) => { console.log("Testing absolut distance with points:"+ point1 + ";" + point2);

/* Offset points by one km */
const offsetPoints = turf.processSegment(point1, point2, KmInDegrees);

const offsetPoint1 = offsetPoints[0];
const offsetPoint2 = offsetPoints[1];

/* Distances from original points to new points */
const distance1 = turf.distance(point1, offsetPoint1);
const distance2 = turf.distance(point2, offsetPoint2);

/* Distances should be one km within half a meter (0.05% error) */
console.assert(IsWithin(distance1, 1, 0.5*MeterInKm), distance1);
console.assert(IsWithin(distance2, 1, 0.5*MeterInKm), distance2);

};

const TestDistanceFromPoints = (point1, point2) => { console.log("Testing distance with points:"+ point1 + ";" + point2);

/* Offset points by one km */
const offsetPoints = turf.processSegment(point1, point2, KmInDegrees);

const offsetPoint1 = offsetPoints[0];
const offsetPoint2 = offsetPoints[1];

/* Distances from original points to new points */
const distance1 = turf.distance(point1, offsetPoint1);
const distance2 = turf.distance(point2, offsetPoint2);

/* Distances should be same within half a meter (0.05% error) */
console.assert(IsWithin(distance1, distance2, 0.5*MeterInKm), distance1, distance2);

};

const TestBearing = (point1, point2) => { console.log("Testing bearing with points:" + point1 + ";" + point2);

/* Offset points by one km */
const offsetPoints = turf.processSegment(point1, point2, KmInDegrees);

const offsetPoint1 = offsetPoints[0];
const offsetPoint2 = offsetPoints[1];

const bearing1 = turf.bearing(point1, point2);
const bearing2 = turf.bearing(offsetPoint1, offsetPoint2);

/* Bearing between new points should be same as between old points within 0.01 degree */
console.assert(IsWithin(bearing1, bearing2, 0.01));

};

TestAbsoluteDistance([0.0, 0.0], [0.0, 0.0 + KmInDegrees]); TestAbsoluteDistance([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]); TestAbsoluteDistance([0.0, 0.0], [0.0 + KmInDegrees, 0.0]); TestAbsoluteDistance([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestAbsoluteDistance([0.0, 45.0], [0.0, 45.0 + KmInDegrees]); TestAbsoluteDistance([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]); TestAbsoluteDistance([0.0, 45.0], [0.0 + KmInDegrees, 45.0]); TestAbsoluteDistance([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestAbsoluteDistance([0.0, 65.0], [0.0, 65.0 + KmInDegrees]); TestAbsoluteDistance([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]); TestAbsoluteDistance([0.0, 65.0], [0.0 + KmInDegrees, 65.0]); TestAbsoluteDistance([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestAbsoluteDistance([30.0, 65.0], [30.0, 65.0 + KmInDegrees]); TestAbsoluteDistance([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]); TestAbsoluteDistance([30.0, 65.0], [30.0 + KmInDegrees, 65.0]); TestAbsoluteDistance([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]);

TestDistanceFromPoints([0.0, 0.0], [0.0, 0.0 + KmInDegrees]); TestDistanceFromPoints([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]); TestDistanceFromPoints([0.0, 0.0], [0.0 + KmInDegrees, 0.0]); TestDistanceFromPoints([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestDistanceFromPoints([0.0, 45.0], [0.0, 45.0 + KmInDegrees]); TestDistanceFromPoints([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]); TestDistanceFromPoints([0.0, 45.0], [0.0 + KmInDegrees, 45.0]); TestDistanceFromPoints([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestDistanceFromPoints([0.0, 65.0], [0.0, 65.0 + KmInDegrees]); TestDistanceFromPoints([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]); TestDistanceFromPoints([0.0, 65.0], [0.0 + KmInDegrees, 65.0]); TestDistanceFromPoints([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestDistanceFromPoints([30.0, 65.0], [30.0, 65.0 + KmInDegrees]); TestDistanceFromPoints([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]); TestDistanceFromPoints([30.0, 65.0], [30.0 + KmInDegrees, 65.0]); TestDistanceFromPoints([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]);

TestBearing([0.0, 0.0], [0.0, 0.0 + KmInDegrees]); TestBearing([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]); TestBearing([0.0, 0.0], [0.0 + KmInDegrees, 0.0]); TestBearing([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestBearing([0.0, 45.0], [0.0, 45.0 + KmInDegrees]); TestBearing([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]); TestBearing([0.0, 45.0], [0.0 + KmInDegrees, 45.0]); TestBearing([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestBearing([0.0, 65.0], [0.0, 65.0 + KmInDegrees]); TestBearing([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]); TestBearing([0.0, 65.0], [0.0 + KmInDegrees, 65.0]); TestBearing([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestBearing([30.0, 65.0], [30.0, 65.0 + KmInDegrees]); TestBearing([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]); TestBearing([30.0, 65.0], [30.0 + KmInDegrees, 65.0]); TestBearing([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]); `

So, the TestAbsoluteDistance checks that the actual offset is correct, and the TestDistanceFromPoints checks that the two new points are at the same distance from their orignal points (on a local scale), and the TestBearing checks that the line that crosses the two new points has the same bearing as the old line (on a local scale again). Do you guys think we need some more tests?

Is is possible to add some comment to the lineOffset function that is is only correct locally (up to some 10s kms)?

Greetings, Martin

PereUbu7 avatar Nov 19 '23 20:11 PereUbu7

I think it's a good idea, but in that case we should mark the old method as deprecated and then try to implement the geodesic version correctly.

This PR contains the planar version. The old (current) method seems to be planar too, but seems only to be correct close to the equator.

PereUbu7 avatar Aug 06 '24 11:08 PereUbu7

@PereUbu7 this should get you started with the tests, it should copy/paste into packages/turf-line-offset/test.ts cleanly. You'll have to add @turf/distance and @turf/bearing as devDependencies to turf-line-offset. The only question I have is whether we want to move processSegment into its own file so that we can export it for testing, or if you can just rework the tests to use the public API of lineOffset instead of hitting processSegment directly. I'd prefer to just test the public API if we can get away with it.

const IsWithin = (value1: number, value2: number, variance: number) => {
  return value1 > value2 - variance && value1 < value2 + variance;
};

const MeterInKm = 0.001;
const KmInDegrees = 0.00899;

const TestAbsoluteDistance = (
  point1: [number, number],
  point2: [number, number]
) => {
  console.log("Testing absolut distance with points:" + point1 + ";" + point2);

  /* Offset points by one km */
  const offsetPoints = processSegment(point1, point2, KmInDegrees);

  const offsetPoint1 = offsetPoints[0];
  const offsetPoint2 = offsetPoints[1];

  /* Distances from original points to new points */
  const distance1 = distance(point1, offsetPoint1);
  const distance2 = distance(point2, offsetPoint2);

  test(`TestAbsoluteDistance - ${point1}; ${point2}`, (t) => {
    /* Distances should be one km within half a meter (0.05% error) */
    t.ok(IsWithin(distance1, 1, 0.5 * MeterInKm));
    t.ok(IsWithin(distance2, 1, 0.5 * MeterInKm));
    t.end();
  });
};

const TestDistanceFromPoints = (
  point1: [number, number],
  point2: [number, number]
) => {
  /* Offset points by one km */
  const offsetPoints = processSegment(point1, point2, KmInDegrees);

  const offsetPoint1 = offsetPoints[0];
  const offsetPoint2 = offsetPoints[1];

  /* Distances from original points to new points */
  const distance1 = distance(point1, offsetPoint1);
  const distance2 = distance(point2, offsetPoint2);

  test(`TestDistanceFromPoints ${point1};${point2}`, (t) => {
    /* Distances should be same within half a meter (0.05% error) */
    t.ok(IsWithin(distance1, distance2, 0.5 * MeterInKm));
    t.end();
  });
};

const TestBearing = (point1: [number, number], point2: [number, number]) => {
  console.log("Testing bearing with points:" + point1 + ";" + point2);

  /* Offset points by one km */
  const offsetPoints = processSegment(point1, point2, KmInDegrees);

  const offsetPoint1 = offsetPoints[0];
  const offsetPoint2 = offsetPoints[1];

  const bearing1 = bearing(point1, point2);
  const bearing2 = bearing(offsetPoint1, offsetPoint2);

  test(`TestBearing ${point1};${point2}`, (t) => {
    /* Bearing between new points should be same as between old points within 0.01 degree */
    t.ok(IsWithin(bearing1, bearing2, 0.01));
    t.end();
  });
};

TestAbsoluteDistance([0.0, 0.0], [0.0, 0.0 + KmInDegrees]);
TestAbsoluteDistance([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]);
TestAbsoluteDistance([0.0, 0.0], [0.0 + KmInDegrees, 0.0]);
TestAbsoluteDistance([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestAbsoluteDistance([0.0, 45.0], [0.0, 45.0 + KmInDegrees]);
TestAbsoluteDistance([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]);
TestAbsoluteDistance([0.0, 45.0], [0.0 + KmInDegrees, 45.0]);
TestAbsoluteDistance([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestAbsoluteDistance([0.0, 65.0], [0.0, 65.0 + KmInDegrees]);
TestAbsoluteDistance([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]);
TestAbsoluteDistance([0.0, 65.0], [0.0 + KmInDegrees, 65.0]);
TestAbsoluteDistance([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestAbsoluteDistance([30.0, 65.0], [30.0, 65.0 + KmInDegrees]);
TestAbsoluteDistance([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]);
TestAbsoluteDistance([30.0, 65.0], [30.0 + KmInDegrees, 65.0]);
TestAbsoluteDistance([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]);

TestDistanceFromPoints([0.0, 0.0], [0.0, 0.0 + KmInDegrees]);
TestDistanceFromPoints([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]);
TestDistanceFromPoints([0.0, 0.0], [0.0 + KmInDegrees, 0.0]);
TestDistanceFromPoints([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestDistanceFromPoints([0.0, 45.0], [0.0, 45.0 + KmInDegrees]);
TestDistanceFromPoints([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]);
TestDistanceFromPoints([0.0, 45.0], [0.0 + KmInDegrees, 45.0]);
TestDistanceFromPoints([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestDistanceFromPoints([0.0, 65.0], [0.0, 65.0 + KmInDegrees]);
TestDistanceFromPoints([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]);
TestDistanceFromPoints([0.0, 65.0], [0.0 + KmInDegrees, 65.0]);
TestDistanceFromPoints([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestDistanceFromPoints([30.0, 65.0], [30.0, 65.0 + KmInDegrees]);
TestDistanceFromPoints([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]);
TestDistanceFromPoints([30.0, 65.0], [30.0 + KmInDegrees, 65.0]);
TestDistanceFromPoints([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]);

TestBearing([0.0, 0.0], [0.0, 0.0 + KmInDegrees]);
TestBearing([0.0 + KmInDegrees, 0.0], [0.0, 0.0 + KmInDegrees]);
TestBearing([0.0, 0.0], [0.0 + KmInDegrees, 0.0]);
TestBearing([0.0, 0.0 + KmInDegrees], [0.0 + KmInDegrees, 0.0]);

TestBearing([0.0, 45.0], [0.0, 45.0 + KmInDegrees]);
TestBearing([0.0 + KmInDegrees, 45.0], [0.0, 45.0 + KmInDegrees]);
TestBearing([0.0, 45.0], [0.0 + KmInDegrees, 45.0]);
TestBearing([0.0, 45.0 + KmInDegrees], [0.0 + KmInDegrees, 45.0]);

TestBearing([0.0, 65.0], [0.0, 65.0 + KmInDegrees]);
TestBearing([0.0 + KmInDegrees, 65.0], [0.0, 65.0 + KmInDegrees]);
TestBearing([0.0, 65.0], [0.0 + KmInDegrees, 65.0]);
TestBearing([0.0, 65.0 + KmInDegrees], [0.0 + KmInDegrees, 65.0]);

TestBearing([30.0, 65.0], [30.0, 65.0 + KmInDegrees]);
TestBearing([30.0 + KmInDegrees, 65.0], [30.0, 65.0 + KmInDegrees]);
TestBearing([30.0, 65.0], [30.0 + KmInDegrees, 65.0]);
TestBearing([30.0, 65.0 + KmInDegrees], [30.0 + KmInDegrees, 65.0]);

mfedderly avatar Aug 06 '24 15:08 mfedderly