Opened 5 years ago

Closed 5 years ago

#2534 closed defect (fixed)

st_distance is returning incorrect results for large geographies

Reported by: samuelb Owned by: pramsey
Priority: medium Milestone: PostGIS 2.1.2
Component: postgis Version: 2.0.x
Keywords: st_distance, geography Cc:

Description

Following geography distance query:

SELECT st_distance(st_geographyfromtext('POINT(-0 0)'), st_geographyfromtext('POLYGON((-20 -20,-20 60,90 60,90 -20,-20 -20))'))

returns 2226389.81558413 instead of the expected result 0 as you get for the same query but with geometries:

SELECT st_distance(st_geomfromewkt('SRID=4326;POINT(-0 0)'), st_geomfromewkt('SRID=4326;POLYGON((-20 -20,-20 60,90 60,90 -20,-20 -20))'))

I've tried adding extra vertices to my rectangle but that didn't help.

When you slightly decrease to a maximum latitude of 55 then it works:

SELECT st_distance(st_geographyfromtext('POINT(-0 0)'), st_geographyfromtext('POLYGON((-20 -20,-20 55,90 55,90 -20,-20 -20))'))

Changing the order of the coordinates didn't work. In the end I would like to be able to check if a point is within 2km of the bounding box (-179 -85) (179 85) as this is one of the geometries in my database.

Note that st_dwithin also returns false instead of true.

Tested on "PostgreSQL 9.1.10 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu/Linaro? 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS="1.5.3" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.8" USE_STATS"

and on:

"PostgreSQL 9.2.4, compiled by Visual C++ build 1600, 32-bit POSTGIS="2.0.3 r11132" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08 GDAL_DATA not found" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER"

Attachments (2)

2534.patch (715 bytes) - added by pramsey 5 years ago.
prelim sol'n
2534a.patch (1.2 KB) - added by pramsey 5 years ago.
better

Download all attachments as: .zip

Change History (14)

comment:1 Changed 5 years ago by pramsey

It seems to be consistently not finding that the point is in the polygon, so the distance is the distance to the ring.

http://www.gcmap.com/mapui?P=+20S+20W-60N+20W%0D%0A60N+20W-+60N+90E%0D%0A+60N+90E-20S+90E%0D%0A20S+90E-+20S+20W%0D%0A0N+0E&MS=wls&DU=mi

Interesting, the distance POINT(-19 0) and POINT(-21 0) are the same, which means the ring calculation is consistent, but the point-in-polygon calculation is not. Reported on 1.5 and 2.0, still exists in 2.2, so not a regression.

I wouldn't expect your desired use case (creating bounding "boxes" on a sphere) to work particularly well regardless, since you're expecting straight lines where they don't exist (see the great circle image above and look at the top and bottom lines)

comment:2 Changed 5 years ago by samuelb

For the same polygon I also noticed that 0 is returned when a point is near the edge but on the outside of the polygon

e.g. SELECT st_distance(st_geographyfromtext('POINT(-21.4 -21.4)'), st_geographyfromtext('POLYGON((-20 -20,-20 60,90 60,90 -20,-20 -20))'))

SELECT st_distance(st_geographyfromtext('POINT(-21.4 61.8)'), st_geographyfromtext('POLYGON((-20 -20,-20 60,90 60,90 -20,-20 -20))'))

comment:3 Changed 5 years ago by pramsey

Milestone: PostGIS 2.1.1PostGIS 2.1.2

Changed 5 years ago by pramsey

Attachment: 2534.patch added

prelim sol'n

comment:5 Changed 5 years ago by pramsey

I'm really rather shocked this didn't come to light sooner, and now I need to backtrack a bit and see how long this problem has been around. Basically a couple algorithms attempt to test edge crossing without actually computing intersection points. And the code there now mistakes antipodal edges as crossing cases. Amazing. What idiot wrote that? (oh, right)

comment:6 in reply to:  5 Changed 5 years ago by pramsey

Replying to pramsey:

Now I need to backtrack a bit and see how long this problem has been around.

It's been around since 1.5, though it got used a little more widely in 2.1, so more opportunities to showcase the mistake.

comment:7 Changed 5 years ago by pramsey

Need some pretty big tests to make sure this fix doesn't break something else.

Changed 5 years ago by pramsey

Attachment: 2534a.patch added

better

comment:8 Changed 5 years ago by pramsey

I'm more confident of this patch, conceptually.

comment:9 Changed 5 years ago by pramsey

I've patched trunk, so it's easier to run compare/contrast tests. Trunk has all the fixes 2.1 does, plus this one. So it should be a fair test. r12228

comment:10 Changed 5 years ago by robe

okay tested with trunk and get 0 as expected. I'm doing a garden regress against 2.2 / 2.0 / 2.1 to make sure there are no unplanned surprises from the change.

comment:11 Changed 5 years ago by robe

garden tests look about the same as 2.1 except for 3 differences in ST_3DExtent, ST_Within, ST_Intersects, but these I think are just instability issues than anything different since they sometimes give different answers in 2.1 depending on run.

I'l try to find some real world geomtries to exercise with next.

comment:12 Changed 5 years ago by pramsey

Resolution: fixed
Status: newclosed

Applied to 2.1 as well at r12308

Note: See TracTickets for help on using tickets.