lineSlice is very imprecise

[crubier]

On turf 6.3.0: lineSlice is very imprecise:

      // Start with a line to slice
      const line = {
        "type": "Feature",
        "properties": {
          "index": 0
        },
        "geometry": {
          "type": "LineString",
          "coordinates": [
            [
              3.69140625,
              51.72702815704774
            ],
            [
              -5.3173828125,
              41.60722821271717
            ]
          ]
        }
      }

     // The start point of the line (Just the first coordinate of the line)
     const pointStart = {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Point",
          "coordinates": [
            3.69140625,
            51.72702815704774
          ]
        }
      }

      // The end point of our slice (obtained by intersecting the line with another line with lineIntersect(line, splitter) )
      const pointEnd = {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Point",
          "coordinates": [
            0.31936718356317106,
            47.93913163509963
          ]
        }
      }
      
      // Double check that this point is on the line :
      pointToLineDistance(point, segment, { units: "meters" }) < 0.001 ;
      // true
      
      // Now run lineSlice
      lineSlice( pointStart, pointEnd, line )
      // It returns this result:
      {
        "type": "Feature",
        "properties": {
          "index": 0
        },
        "geometry": {
          "type": "LineString",
          "coordinates": [
            // The first coordinate is the start point, normal
            [
              3.69140625,
              51.72702815704774
            ],
            ////////////////////////////////////////////////////////////////////////////////////
            // The second coordinate below is VERY far 
            // from the slice point we asked for, even though we 
            // showed it was on the line
            [
              -0.13132027083960374,
              47.4328630517935
            ]
            // It should be this instead
            // [
            //  0.31936718356317106,
            //  47.93913163509963
            // ]
            ////////////////////////////////////////////////////////////////////////////////////
          ]
        }
      }

[crubier]

This is caused by inaccuracy in nearestPointOnLine

• length and nearestPointOnLine - accuracy #1726
• nearestPointOnLine deviation #1440
• nearestPointOnLine uses wrong calculation #1878

[JamesLMilner]

Let me know if there's anything I can do to help with helping improve the inaccuracy of nearestPointOnLine

[crubier]

@JamesLMilner as I said in #1440 (comment) , this is the big blocker to use turf for anything serious. I fell back to python / qgis for my application, even if it's less convenient in a JS app.

I think a good starting point would be to add my code above as a test case, and then modify nearestPointOnLine until the test passes. A similar case is explained here #1726 (comment) .

Then I am not sure if the root problem is in nearestPointOnLine or in one of its dependencies (bearing, distance, destination, lineIntersects). But I believe there is an inconsistency somewhere in there, where some of these functions work on a spherical referential, and some others expect a cartesian referential, or something like that, (I don't think floating point errors can create a 50km difference on 500km lines intersection like in my test case above).

If I had time, I'd fix it myself by looking at these 5 functions, and implement them "by the book", applying the mathematical formula directly, maybe using Mathematica to derive a general closed-form solution to each of these functions. With a clearly documented and consistent hypothesis on the geoid model being used.

[JamesLMilner]

To clarify I'm not on the Turf.js team - however I also rely on it for some important parts of our project and am currently experiencing issues with the nearestPointOnLinewhich seems to be a (root?) cause for this issue.

One thing I'd like to understand is what nearestPointOnLine is meant to determine i.e. is it the nearest point on a great circle line or a rhumb line or a line in Cartesian coordinates in Web Mercator for example? I feel like there should probably be two separate implementations as 'nearest point on a line' is somewhat ambiguous in relation to geographic coordinates.

I've implemented Chris Vernesses nearestPointOnSegment, which takes nearest point on a great circle line between two points, as a pure function in JavaScript (technically TypeScript), but would be happy to share if it would be useful here.

[stebogit]

@JamesLMilner I'm not familiar with the details on this issue, but just to answer your question, @turf/nearest-point-on-line takes a Point P and a LineString L and returns the closest Point C on/along L that is closest to P.

As far as I can tell it uses internally @turf/distance, ergo it should calculate the closest distance using the Haversine formula to account for global curvature. Now, I would say that the great circle line is the shortest line on a sphere between two points (i.e. represents the distance calculated with the Haversine formula), but I might be wrong.

[crubier]

@stebogit you said this:

Point P and a LineString L and returns the closest Point C on/along L that is closest to P.

The problem is that this is wrong! In particular, the Point C is NOT along LineString L , it is NOT on the line, which is clearly a bug, given that the name of the function is nearestPointONLINE.

See this comment who expresses the problem clearly #1726 (comment)

Basically if we have point from "nearestPointOnLine" and check it with "booleanPointOnLine", we will get "false"...

And this is one aspect of the general problem expressed here #1440 (comment)

[crubier]

And @JamesLMilner yes you are right, I think some parts of Turf use Rhumb, some other use great circle (Haversine). This is implicit and never expressed clearly, so it leads to bugs and inconsistencies...

[rowanwins]

@JamesLMilner if you could share your implementation that would be great. It would be interesting to see how large & different the implementation is from what's currently in Turf. Potentially if it's small enough we could add an option to the existing module as to whether we use one implementation or the other, we'll just need to see what makes sense I guess...

[JamesLMilner]

@rowanwins here's the implementation - it uses [longittude, latatutide] ordering as per the GeoJSON spec and has no dependencies - inputs are coordinates rather than GeoJSON points/lines, but that could easily be changed:

function equals(coord0: [number, number], coord1: [number, number]) {
    if (Math.abs(coord0[1] - coord1[1]) > Number.EPSILON) return false;
    if (Math.abs(coord0[0] - coord1[0]) > Number.EPSILON) return false;

    return true;
}

const toRadians = (lngLat: number) => {
    return (lngLat * Math.PI) / 180;
};

function cross(first: { x: number; y: number; z: number }, v: { x: number; y: number; z: number }) {
    const x = first.y * v.z - first.z * v.y;
    const y = first.z * v.x - first.x * v.z;
    const z = first.x * v.y - first.y * v.x;

    return { x, y, z };
}

function toNvector(coord: [number, number]) {

    const φ = toRadians(coord[1]);
    const λ = toRadians(coord[0]);

    const sinφ = Math.sin(φ),
        cosφ = Math.cos(φ);

    const sinλ = Math.sin(λ),
        cosλ = Math.cos(λ);

    // right-handed vector: x -> 0°E,0°N; y -> 90°E,0°N, z -> 90°N
    const x = cosφ * cosλ;
    const y = cosφ * sinλ;
    const z = sinφ;

    return { x, y, z };
}

function minus(first: { x: number; y: number; z: number }, v: { x: number; y: number; z: number }) {
    return { x: first.x - v.x, y: first.y - v.y, z: first.z - v.z };
}

function dot(first: { x: number; y: number; z: number }, v: { x: number; y: number; z: number }) {
    return first.x * v.x + first.y * v.y + first.z * v.z;
}

function isWithinExtent(coord0: [number, number], coord1: [number, number], coord2: [number, number]) {
    if (equals(coord1, coord2)) {
        return equals(coord0, coord1); // null segment
    }

    const n0 = toNvector(coord0),
        n1 = toNvector(coord1),
        n2 = toNvector(coord2); // n-vectors

    // get vectors representing p0->p1, p0->p2, p1->p2, p2->p1
    const δ10 = minus(n0, n1),
        δ12 = minus(n2, n1);
    const δ20 = minus(n0, n2),
        δ21 = minus(n1, n2);

    // dot product δ10⋅δ12 tells us if p0 is on p2 side of p1, similarly for δ20⋅δ21
    const extent1 = dot(δ10, δ12);
    const extent2 = dot(δ20, δ21);

    const isSameHemisphere = dot(n0, n1) >= 0 && dot(n0, n2) >= 0;

    return extent1 >= 0 && extent2 >= 0 && isSameHemisphere;
}

const toDegrees = function (vector3d: number) {
    return (vector3d * 180) / Math.PI;
};

function nVectorToLatLon(vector: { x: number; y: number; z: number }) {
    // tanφ = z / √(x²+y²), tanλ = y / x (same as ellipsoidal calculation)

    const x = vector.x,
        y = vector.y,
        z = vector.z;

    const φ = Math.atan2(z, Math.sqrt(x * x + y * y));
    const λ = Math.atan2(y, x);

    return [toDegrees(λ), toDegrees(φ)];
}

function vectorLength(v: { x: number; y: number; z: number }) {
    return Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

export function distanceTo(coord1: [number, number], coord2: [number, number], radius = 6371e3) {
    const R = Number(radius);

    const n1 = toNvector(coord1);
    const n2 = toNvector(coord2);

    const sinθ = vectorLength(cross(n1, n2));
    const cosθ = dot(n1, n2);
    const δ = Math.atan2(sinθ, cosθ); // tanδ = |n₁×n₂| / n₁⋅n₂

    return δ * R;
}

/**
 * Returns closest coordinate on great circle segment between lineCoordOne & lineCoordTwo to coord.
 *
 * If this coord is ‘within’ the extent of the segment, the coord is on the segment between coord1 &
 * coord2; otherwise, it is the closer of the endcoords defining the segment.
 */
export function nearestCoordinateOnSegment(
    coord: [number, number],
    lineCoordOne: [number, number],
    lineCoordTwo: [number, number]
) {
    let closestCoords = null;

    const isBetweenLineCoords = isWithinExtent(coord, lineCoordOne, lineCoordTwo);
    const isNotEqualToLineCoords = !equals(lineCoordOne, lineCoordTwo);

    if (isBetweenLineCoords && isNotEqualToLineCoords) {
        // closer to segment than to its endcoords, find closest coord on segment
        const n0 = toNvector(coord),
            n1 = toNvector(lineCoordOne),
            n2 = toNvector(lineCoordTwo);
        const c1 = cross(n1, n2); // n1×n2 = vector representing great circle through p1, p2
        const c2 = cross(n0, c1); // n0×c1 = vector representing great circle through p0 normal to c1
        const n = cross(c1, c2); // c2×c1 = nearest coord on c1 to n0

        closestCoords = nVectorToLatLon(n);
    } else {
        // beyond segment extent, take closer endcoord
        const d1 = distanceTo(coord, lineCoordOne);
        const d2 = distanceTo(coord, lineCoordTwo);
        closestCoords = d1 < d2 ? lineCoordOne : lineCoordTwo;
    }

    return closestCoords;
}

It passes the unit tests in the original code.

Interestingly compared to the Turf implementation here is the difference:

Above implementation results:

1.9000033116244113,
51.00038411380565

Turf

 1.9
 51

I would assume this means that Turf implementation is trying to find the point on a Great Circle as per the above implementation. The precision on the latitude there concerns me slightly as 0.0003 degrees is probably in the region of 33 meters of error.

[smallsaucepan]

Hi @crubier. Testing a fix for this issue, and have some comments to add.

// Double check that this point is on the line :
pointToLineDistance(point, segment, { units: "meters" }) < 0.001 ;
// true

pointToLineDistance suffered from a similar issue to nearestPointOnLine (which impacted lineSlice). So the above wasn't actually true.

Looking at the below close up of your example, the blue pin is the point mentioned in the double check above, the red line is the path of the line you want to split, and the red pin is the second coordinate that was VERY far from the slice point:

With the fix mentioned, the green pin is now the slice point, and the yellow line coming in from the top is the slice you would get. Hope that all makes sense. If that's not what you would expect, please let me know.