Incorrect distance between geography polygons

[Mike Taves]

Looking into ​https://lists.osgeo.org/pipermail/postgis-users/2024-July/046496.html there appears to be issues with a few distance functions on polygons with geographic coords. Running the same query on a few docker images it appears this been broken since version 3.1.

Here is a query that does several distance methods between two valid/simple polygons:

SELECT ST_Distance(ST_Transform(A, 27700), ST_Transform(B, 27700)) AS distance_27700,
  ST_DistanceSpheroid(A, B, 'SPHEROID["WGS 84",6378137,298.257223563]') AS DistanceSpheroid,
  ST_DistanceSphere(A, B) AS DistanceSphere,
  ST_Distance(A::geography, B::geography, true) AS distance_geography_spheroid,
  ST_Distance(A::geography, B::geography, false) AS distance_geography_sphere
FROM (
 SELECT
  ST_GeomFromText('POLYGON ((-1.7485534738529238 52.477989324720184, -1.7485528426344397 52.47797089391474, -1.748563887653701 52.47797046786393, -1.7485727216213993 52.47797048661752, -1.7485733502887397 52.47798936691661, -1.7485739840754728 52.478007348228196, -1.748563677771212 52.47800732634904, -1.7485541076315525 52.47800730603188, -1.7485534738529238 52.477989324720184))', 4326) AS A,
  ST_GeomFromText('POLYGON ((-1.748665164272404 52.47789651546308, -1.748705653223533 52.477896601376905, -1.7487057222918938 52.47788446504611, -1.7487307518185344 52.47788451814952, -1.7487299491521497 52.47789620342481, -1.7487660211267058 52.477896279947124, -1.748767153309673 52.477956065736684, -1.7486648265496951 52.477955848635304, -1.7485639721184771 52.477955634571074, -1.7484653236191374 52.47795587460495, -1.748465659049668 52.47789699092737, -1.7485017310241382 52.477897067530265, -1.7485018001486063 52.47788493119968, -1.7485297717654242 52.47788544009214, -1.7485297103282949 52.477896227941585, -1.748562835091827 52.47789674776739, -1.7485996432291142 52.47789682590328, -1.7485997123266428 52.47788468957258, -1.7486276865033974 52.47788474894817, -1.748627619972449 52.477896435785155, -1.748665164272404 52.47789651546308))', 4326) AS B
) AS data;

With PostGIS 3.0 (via postgis/postgis:10-3.0) here is the expected result:

-[ RECORD 1 ]---------------+-----------------
distance_27700              | 1.64995405783897
distancespheroid            | 1.65055968246111
distancesphere              | 1.64935316
distance_geography_spheroid | 1.65055968
distance_geography_sphere   | 1.64935316

whereas with PostGIS 3.1 (via postgis/postgis:10-3.1) here is the malformed result:

-[ RECORD 1 ]---------------+-----------------
distance_27700              | 1.6499540575588
distancespheroid            | 1.65055968246111
distancesphere              | 0
distance_geography_spheroid | 0
distance_geography_sphere   | 0

and similar results with the latest docker images (PostGIS 3.4.2):

-[ RECORD 1 ]---------------+-------------------
distance_27700              | 1.649954057558797
distancespheroid            | 1.6505596824611086
distancesphere              | 0
distance_geography_spheroid | 0
distance_geography_sphere   | 0

Based on the above, there appears to be some issues with:

• ST_Distance(geography, geography)
• ST_DistanceSphere(geometry, geometry)

And to be clear, there are no issues with:

• ST_Distance(geometry, geometry)
• ST_DistanceSpheroid(geometry, geometry, spheroid)

[Mike Taves]

QGIS screenshot

[Mike Taves]

Here is a screenshot of the geometries in QGIS with a manually measured distance. The base projection is EPSG:27700

[pramsey]

Minimal reproducer

SELECT ST_Distance(ST_Transform(A, 27700), ST_Transform(B, 27700)) AS distance_27700,
  ST_DistanceSpheroid(A, B, 'SPHEROID["WGS 84",6378137,298.257223563]') AS DistanceSpheroid,
  ST_DistanceSphere(A, B) AS DistanceSphere,
  ST_Distance(A::geography, B::geography, true) AS distance_geography_spheroid,
  ST_Distance(A::geography, B::geography, false) AS distance_geography_sphere
FROM (
 SELECT
  ST_GeomFromText('LINESTRING (
-1.7485638876537 52.477970467863,
-1.7485727216213 52.477970486617
)', 4326) AS A,
  ST_GeomFromText('LINESTRING (
-1.7486648265496 52.477955848635,
-1.7485639721184 52.477955634571
)', 4326) AS B
) AS data;

[pramsey]

So the minimal case shows the problem in edge_intersects, showing up because the two edges in this case are parallel and very close together, and they are being erroneously flagged as co-linear edges (distance 0).

[pramsey]

Fixed in 4f33aa5bc/git. Applied back to stable-3.1