PostGIS 3, ST_DWithin Query - "lwgeom_distance_spheroid returned negative!" by changing third arg? Postgres 11, Ubuntu 20.04

[dchun]

The unexpected error:

lwgeom_distance_spheroid returned negative! By making a tiny change in the third argument (distance) to ST_DWithin, we get an error in one version but not another. In other words, by changing the distance by a tiny fraction: 3754219.5 → 3754219.6 we get the error.

This distance where it starts failing is different for different POINTs.

Additionally, geometry or geography makes no difference, it only changes the value where the distance error occurs.

It does not happen on Ubuntu 16.04, Postgres 11, PostGIS 2.5.

I've tried to isolate it here:

FAILING ON Ubuntu 20.04 PostGIS 3.3:

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 20.04.5 LTS
Release:        20.04
Codename:       focal

psql (11.19 (Ubuntu 11.19-1.pgdg20.04+1))

# SELECT VERSION();
                                                               version                                                               
-------------------------------------------------------------------------------------------------------------------------------------
 PostgreSQL 11.19 (Ubuntu 11.19-1.pgdg20.04+1) on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0, 64-bit
(1 row)

# SELECT postgis_full_version();
                                                                        postgis_full_version                                                                         
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="3.3.2 4975da8" [EXTENSION] PGSQL="110" GEOS="3.9.1-CAPI-1.14.2" PROJ="6.3.1" LIBXML="2.9.10" LIBJSON="0.13.1" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"
(1 row)

# SELECT PostGIS_Liblwgeom_Version();
 postgis_liblwgeom_version 
---------------------------
 3.3.2 4975da8
(1 row)

# SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 3754219.5, true) as is_within;
 is_within 
-----------
 f
(1 row)

# SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 3754219.6, true) as is_within;
ERROR:  lwgeom_distance_spheroid returned negative!
CONTEXT:  SQL function "st_dwithin" statement 1

SUCCESSFUL ON Ubuntu 16.04 PostGIS 2.5:

lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 16.04.6 LTS
Release:        16.04
Codename:       xenial

(12.1 (Ubuntu 12.1-1.pgdg16.04+1), server 11.6 (Ubuntu 11.6-1.pgdg16.04+1))
You are now connected to database "mydatabasewithpostgisinit" as user "postgres".
=# SELECT VERSION();
                                                                   version
---------------------------------------------------------------------------------------------------------------------------------------------
 PostgreSQL 11.6 (Ubuntu 11.6-1.pgdg16.04+1) on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609, 64-bit
(1 row)

=# SELECT postgis_full_version();
                                                                                                   postgis_full_version
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="2.5.3 r17699" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" LIBPROTOBUF="1.2.1" RASTER
(1 row)

=# SELECT PostGIS_Liblwgeom_Version();
 postgis_liblwgeom_version
---------------------------
 2.5.3 r17699
(1 row)

=# SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 3754219.5, true) as is_within;
 is_within
-----------
 f
(1 row)

=# SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 3754219.6, true) as is_within;
 is_within
-----------
 f
(1 row)

# SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 111113754219.6, true) as is_within;
 is_within 
-----------
 t

You can see, for a huge distance, in the working version, it still doesn't fail on 2.5.

Is this a bug? Or has something changed? I haven't been able to find any explanation for it.

Thanks!

[robe]

I think this might be an issue with your proj version.

On My PostgreSQL 11, 3.3.2 I do not get any errors.

PostgreSQL 11.17, compiled by Visual C++ build 1914, 64-bit
POSTGIS="3.3.2 3.3.2" [EXTENSION] PGSQL="110" GEOS="3.11.1-CAPI-1.17.1" PROJ="7.2.1" GDAL="GDAL 3.4.3, released 2022/04/22" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" RASTER

For the following queries output

#  SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 111113754219.6, tru
e) as is_within;
 is_within
-----------
 t
(1 row)


#  SELECT ST_DWithin('SRID=4326;POINT(-74.8456 40.1752)','SRID=4326;POINT(-118.493 34.0172)', 3754219.5, true) as
 is_within;
 is_within
-----------
 f

The difference between 2.5 and 3 series I think is when we flipped to using the proj.db that shipped with proj instead of relying on what spatial_ref_sys told us.

[pramsey]

I'd expect if this is a geography distance issue to see things like variable results for:

SELECT ST_Distance('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)'::geography, true);
SELECT _ST_DistanceTree('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)'::geography);
SELECT _ST_DistanceUncached('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)'::geography, true);

Unfortunately for you, I'm not seeing any errors testing on 3.4 on my amazing Apple macbook, but I wouldn't be entirely surprised if you saw some.

[robe]

FWIW this is what I get on my PG 11, PostGIS 3.3.2 install

db=# SELECT ST_Distance('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)'::geograph
y, true);
   st_distance
------------------
 3900084.69544365
(1 row)


db=# SELECT _ST_DistanceTree('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)'::geo
graphy);
 _st_distancetree
------------------
 3900084.69544365
(1 row)


db=# SELECT _ST_DistanceUncached('SRID=4326;POINT(-74.8456 40.1752)'::geography,'SRID=4326;POINT(-118.493 34.0172)':
:geography, true);
 _st_distanceuncached
----------------------
     3900084.69544365
(1 row)

[robe]

From IRC reports, this seems to be an issue with mixing apt.postgresql.org (postgis) and then getting proj and other stuff from ubuntugis, that even though PostGIS is using 6.3, it's hitching a ride on the proj 7 data from ubuntugis, which seems incompatible.

So going to mark this as a wontfix and a PostGIS packaging issue.