SpatiaLite calcutes differently than SQL Server in simple "shortest distance to coast" example

[Ron Grabowski]

We are trying to migrate a distance to coast solution from SQL Server to SpatiaLite but the distances are off by more than 1% and the two lines being drawn are very different.

SQL Server 2019 (https://dbfiddle.uk/):

DECLARE @section as geography= geography::STGeomFromText('LINESTRING(-70.565641 41.936550,-70.566803 41.937228,-70.567408 41.938029,-70.567885 41.938888,-70.568366 41.939728,-70.569003 41.940581,-70.569727 41.941131,-70.570453 41.941709,-70.571622 41.942286,-70.572600 41.942735,-70.573164 41.943099,-70.573327 41.943555,-70.573931 41.943787,-70.574443 41.944150,-70.575010 41.944433,-70.575882 41.944722,-70.576708 41.944958,-70.577261 41.945253,-70.578227 41.945664,-70.578972 41.945694,-70.579366 41.945698,-70.579506 41.946013,-70.579331 41.946415,-70.579087 41.946683,-70.578901 41.947012,-70.579273 41.947234,-70.579615 41.947320,-70.580101 41.947449,-70.580844 41.947918)',4326);

DECLARE @poi as geography = geography::STGeomFromText('POINT (-70.586729 41.929293)', 4326);
DECLARE @shortest_line as geography = @poi.ShortestLineTo(@section);
DECLARE @nearestPoint as geography = geography::STGeomFromText(@shortest_line.STEndPoint().STAsText(), 4326);
SELECT @nearestPoint.STDistance(@poi) as DistanceInMeters
SELECT @poi.ShortestLineTo(@section).STAsText() as ShortestLINESTRING;

DistanceInMeters: 1872
Shortest Line: LINESTRING (-70.586729 41.929292, -70.567065 41.937575)

SpatiaLite 5.0.1:

SELECT GeodesicLength(ST_ShortestLine(
ST_GeomFromText('LINESTRING (-70.565641 41.936550,-70.566803 41.937228,-70.567408 41.938029,-70.567885 41.938888,-70.568366 41.939728,
-70.569003 41.940581,-70.569727 41.941131,-70.570453 41.941709,-70.571622 41.942286,-70.572600 41.942735,-70.573164 41.943099,-70.573327 41.943555,
-70.573931 41.943787,-70.574443 41.944150,-70.575010 41.944433,-70.575882 41.944722,-70.576708 41.944958,-70.577261 41.945253,-70.578227 41.945664,
-70.578972 41.945694,-70.579366 41.945698,-70.579506 41.946013,-70.579331 41.946415,-70.579087 41.946683,-70.578901 41.947012,-70.579273 41.947234,
-70.579615 41.947320,-70.580101 41.947449,-70.580844 41.947918)'),ST_GeomFromText('POINT(-70.586729 41.929293)',4326)));

DistanceInMeters: 1921
Shortest Line: LINESTRING(-70.586729 41.929294, -70.579366 41.945698)

This tool visualizes the coastline with the two lines from the different solutions:

https://clydedacruz.github.io/openstreetmap-wkt-playground/

GEOMETRYCOLLECTION(POINT(-70.586729 41.929293),LINESTRING (-70.566803 41.937228,-70.567408 41.938029,-70.567885 41.938888,-70.568366 41.939728,
-70.569003 41.940581,-70.569727 41.941131,-70.570453 41.941709,-70.571622 41.942286,-70.572600 41.942735,-70.573164 41.943099,-70.573327 41.943555,
-70.573931 41.943787,-70.574443 41.944150,-70.575010 41.944433,-70.575882 41.944722,-70.576708 41.944958,-70.577261 41.945253,-70.578227 41.945664,

-70.578972 41.945694,-70.579366 41.945698,-70.579506 41.946013),LINESTRING(-70.567071 41.937583,-70.586729 41.929293),LINESTRING(-70.586729 41.929294,-70.579367 41.945697))

We've tried variations of ST_Length including the use_ellipsoid=true version.

Seems like SQL Server considers the space between two POINTs while SpatiaLite snaps to a specific POINT?

For such a small example set shouldn't the two approaches be within a few meters of each other? This is our first GIS project so we're very much beginners and might be missing an obvious concept.

Thanks,

Ron

[a.fu...@lqt.it]

Hi Ron,

It's not at unusual that different sw implementations could give
different
results for the same problem.
As a general rule open source software is usually expected to be more
accurate than proprietary sw simply because anyone can directly correct
any errors or flaws.

let's go into specific details.
there is absolutely nothing wrong or suspect in GeodesicLength():

SELECT GeodesicLength(
GeomFromText('LINESTRING (-70.586729 41.929292, -70.567065
41.937575)', 4326));
----------------
1872

As you can see the shortest line identified by SQL server has the
expected length measured in metres also on SpatiaLite.

So the real problem is in the supposed shortest line.
Yet another step forward:

SELECT ST_Distance(
ST_GeomFromText('LINESTRING (-70.565641 41.936550,
-70.566803 41.937228,-70.567408 41.938029,
-70.567885 41.938888,-70.568366 41.939728,
-70.569003 41.940581,-70.569727 41.941131,
-70.570453 41.941709,-70.571622 41.942286,
-70.572600 41.942735,-70.573164 41.943099,
-70.573327 41.943555,-70.573931 41.943787,
-70.574443 41.944150,-70.575010 41.944433,
-70.575882 41.944722,-70.576708 41.944958,
-70.577261 41.945253,-70.578227 41.945664,
-70.578972 41.945694,-70.579366 41.945698,
-70.579506 41.946013,-70.579331 41.946415,
-70.579087 41.946683,-70.578901 41.947012,
-70.579273 41.947234,-70.579615 41.947320,
-70.580101 41.947449,-70.580844 41.947918)', 4326),
ST_GeomFromText('POINT(-70.586729 41.929293)',4326), 0);
---------------------------
1923

SpatiaLite definitely confirms that the distance between
the line and the point is 1923 metres and not 1872 as
reported by SQL Server.

Last final checks:

SELECT ST_Distance(
ST_GeomFromText('LINESTRING (-70.565641 41.936550,
-70.566803 41.937228,-70.567408 41.938029,
-70.567885 41.938888,-70.568366 41.939728,
-70.569003 41.940581,-70.569727 41.941131,
-70.570453 41.941709,-70.571622 41.942286,
-70.572600 41.942735,-70.573164 41.943099,
-70.573327 41.943555,-70.573931 41.943787,
-70.574443 41.944150,-70.575010 41.944433,
-70.575882 41.944722,-70.576708 41.944958,
-70.577261 41.945253,-70.578227 41.945664,
-70.578972 41.945694,-70.579366 41.945698,
-70.579506 41.946013,-70.579331 41.946415,
-70.579087 41.946683,-70.578901 41.947012,
-70.579273 41.947234,-70.579615 41.947320,
-70.580101 41.947449,-70.580844 41.947918)', 4326),
ST_GeomFromText('POINT(-70.579366 41.945698)', 4326), 0);
----------------------------
0

this is the nearest point on the line identified by
SpatiaLite, and as we were expecting it lays
exactely on the line.


SELECT ST_Distance(
ST_GeomFromText('LINESTRING (-70.565641 41.936550,
-70.566803 41.937228,-70.567408 41.938029,
-70.567885 41.938888,-70.568366 41.939728,
-70.569003 41.940581,-70.569727 41.941131,
-70.570453 41.941709,-70.571622 41.942286,
-70.572600 41.942735,-70.573164 41.943099,
-70.573327 41.943555,-70.573931 41.943787,
-70.574443 41.944150,-70.575010 41.944433,
-70.575882 41.944722,-70.576708 41.944958,
-70.577261 41.945253,-70.578227 41.945664,
-70.578972 41.945694,-70.579366 41.945698,
-70.579506 41.946013,-70.579331 41.946415,
-70.579087 41.946683,-70.578901 41.947012,
-70.579273 41.947234,-70.579615 41.947320,
-70.580101 41.947449,-70.580844 41.947918)', 4326),
ST_GeomFromText('POINT( -70.567065 41.937575)', 4326), 0);
----------------------------
0.00684

this is instead the nearest point on the line identified
by SQL Server.
Unexpected surprise: it doesn't lays at all over the line.

quick conclusiion: SQL Server seems to be affected by
errors in calculating the minimum distance between a
line and a point, presumably caused by rounding issues
in floating point arithmetics.

bye Sandro

[i-s-o]

Ron,

The way I see it, you are comparing apples to oranges (geography to geometry):

• Your calculations in SQL Server are based on geographic algorithms (~3D trigonometry).
• Your calculations in SpatiaLite are based on geometric algorithms (~2D Cartesian arithmetic).

SpatiaLite provides several geographic-enabled functions (e.g. GeodesicLength() as Sandro mentioned). However,  ST_ShortestLine() is not one of them. It only works with geometries: if you use 4326 coordinates as inputs, they will be treated as Cartesian X and Y values on a 2D-plane: this will give a severely distorted view of your area of interest. Consequently, it will give you an incorrect result. To make ST_ShortestLine() work as intended, you need to first project your geography onto a suitable geometry (e.g. EPSG:6491 for Massachusetts):
ST_AsText( ST_ShortestLine( ST_Transform(poi.geometry, 6491), ST_Transform(coast.geometry, 6491 ) ) )
=> LINESTRING(275745.616864 853621.395165, 277366.362188 854558.919153)
ST_AsText( ST_Transform( ST_ShortestLine( ST_Transform(poi.geometry, 6491), ST_Transform(coast.geometry, 6491) ), 4326) )
=> LINESTRING(-70.586729 41.929293, -70.567065 41.937575)

Or, if all you need is the distance:
ST_Distance(ST_Transform(poi.geometry, 6491), ST_Transform(coast.geometry, 6491))
=> 1872.369257

I attached the picture of the area of interest (in Massachusetts State Plane projection). Applying ST_ShortestLine() on 4326 coordinates will give you the shortest angular distance of 0.018 from your POI, but at this latitude, this angle gives you the pink ellipse (rather than a circle), which touches your coastline at a northerly point, at a farther location than the real closest location. The blue circle is based on constant distance of 1,872m in projected units: it touches the coastline at the shortest location from the POI.

Bonus point: note that the shortest distance calculation above is based on Cartesian (planar) formula. The more accurate geodesic calculation is slightly longer, a little over 5cm:
ST_Length( GeomFromText( 'LINESTRING(275745.616864 853621.395165, 277366.362188 854558.919153)' ) ) = 1872.369257
GeodesicLength( GeomFromText(' LINESTRING(-70.586729 41.929293, -70.567065 41.937575)' ) ) = 1872.421805
(the value in the image is 1872.415, which is based on more precise coordinates)