ST_DWithin using spheroid does not calculate distance properly for geography

[amc6]

It appears that ST_DWithin is not properly calculating distance for geography points. Or at the very least, it is not calculating distance the same as ST_Distance does. (When forced to use a sphere, both ST_DWithin and ST_Distance match.) I've included a query below to demonstrate the problem. It includes several different calculations to illustrate what is going on, but I will summarize the key points here:

The points are Point(1 2) and Point(1 1). According to a 3rd party calculator (​http://www.flymicro.com/records/distcalc.cfm) the distance between these two points is 110575 meters and this matches with the calculation produced by ST_Distance. ST_DWithin appears to be calculating a value of around 111155.5 meters, which does not match the value calculated by ST_Distance with either a spheroid or sphere. However, the function _ST_DWithinUncached appears to be calculating the correct value of 110575.

Query:

select ST_asEWKT(a), ST_asEWKT(b), pg_typeof(a), pg_typeof(b),
ST_Distance(a, b, true) as dist_true, ST_Distance(a, b, false) as dist_false, 
ST_DWithin(a, b, 111155.5, true) as within_111155_5_true, ST_DWithin(a, b, 111155.6, true) as within_111155_6_true, 
ST_DWithin(a, b, 111195.079, false) as within_111195_079_false, ST_DWithin(a, b, 111195.08, false) as within_111195_08_false, 
_ST_DWithinUncached(a, b, 110576, true) as within_uncached_110576, _ST_DWithinUncached(a, b, 110575, true) as within_uncached_110575,
_ST_DWithinUncached(a, b, 111195.08, false) as within_uncached_111195_08_false, _ST_DWithinUncached(a, b, 111195.079, false) as within_uncached_111195_079_false
from (
	select ST_GeogFromText('SRID=4326;POINT(1.0 2.0)') as a, ST_GeogFromText('SRID=4326;POINT(1.0 1.0)') as b
) as points

Results:

      st_asewkt       |      st_asewkt       | pg_typeof | pg_typeof |    dist_true    |   dist_false    | within_111155_5_true | within_111155_6_true | within_111195_079_false | within_111195_08_false | within_uncached_110576 | within_uncached_110575 | within_uncached_111195_08_false | within_uncached_111195_079_false 
----------------------+----------------------+-----------+-----------+-----------------+-----------------+----------------------+----------------------+-------------------------+------------------------+------------------------+------------------------+---------------------------------+----------------------------------
 SRID=4326;POINT(1 2) | SRID=4326;POINT(1 1) | geography | geography | 110575.06481434 | 111195.07973463 | f                    | t                    | f                       | t                      | t                      | f                      | t                               | f

[amc6]

Version Info:

postgis_full_version(): POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d6" PROJ="Rel. 5.0.1, April 1st, 2018" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.9.7" LIBJSON="0.12.1" RASTER

version(): PostgreSQL 10.3 on x86_64-apple-darwin17.3.0, compiled by Apple LLVM version 9.0.0 (clang-900.0.39.2), 64-bit

OSX 10.13.2

Both postgres and postgis were installed via homebrew

[pramsey]

Hm, so it seems completely contained in the tree distance calculation when used with a non-zero tolerance. So the stopping condition code must have an issue…

select 
ST_asEWKT(a), 
ST_asEWKT(b), 
pg_typeof(a), 
pg_typeof(b),
ST_Distance(a, b, true) as dist_true, 
ST_Distance(a, b, false) as dist_false, 
_ST_DistanceTree(a, b, 0, true) as dist_tree_true,
_ST_DistanceTree(a, b, 0, false) as dist_tree_false,
_ST_DistanceUnCached(a, b, 0, true) as dist_uncached_true,
_ST_DistanceUnCached(a, b, 0, false) as dist_uncached_false,
ST_DWithin(a, b, 111155.5, true) as within_111155_5_true, 
ST_DWithin(a, b, 111155.6, true) as within_111155_6_true, 
ST_DWithin(a, b, 111195.079, false) as within_111195_079_false, 
ST_DWithin(a, b, 111195.08, false) as within_111195_08_false, 
_ST_DWithinUncached(a, b, 110576, true) as within_uncached_110576,
_ST_DWithinUncached(a, b, 110575, true) as within_uncached_110575,
_ST_DWithinUncached(a, b, 111195.08, false) as within_uncached_111195_08_false, 
_ST_DWithinUncached(a, b, 111195.079, false) as within_uncached_111195_079_false
from (
	select ST_GeogFromText('SRID=4326;POINT(1.0 2.0)') as a, ST_GeogFromText('SRID=4326;POINT(1.0 1.0)') as b
) as points

I'm surprised that the spheroid/sphere don't both show the problem, but maybe it's just an artefact of the particular error, which is probably some kind of precision issue we're seeing with just this case. Thanks for the very small test case.

[pramsey]

Aha, it's not the distance calculation, it's the index call, so all the uncached calls work fine, since they avoid the SQL index wrapper, but the calls to ST_DWithin include the index optimization, which is returning the wrong answer:

select
a && _ST_Expand(b,111155.5) as index_b_11155_5,
b && _ST_Expand(a,111155.5) as index_a_11155_5,
a && _ST_Expand(b,111155.6) as index_b_11155_6,
b && _ST_Expand(a,111155.6) as index_a_11155_6
from (
	select ST_GeogFromText('SRID=4326;POINT(1.0 2.0)') as a, ST_GeogFromText('SRID=4326;POINT(1.0 1.0)') as b
) as points

[pramsey]

In 16668:

Ensure index filters on expanded boxes are large enough to encompass
the radii they are searching, closes #4081

[pramsey]

In 16669:

Ensure index filters on expanded boxes are large enough to encompass
the radii they are searching, references #4081