[postgis-users] Inside / Outside of geography polygons

Mon Aug 12 09:55:20 PDT 2019

You are exposing a limitation in the postgis geodetic implementation, namely that we don’t handle things larger than a hemisphere well. This is not unlike the SQL Server 2008 limitation, except we don’t toss errors at you.

https://alastaira.wordpress.com/2012/01/27/ring-orientation-bigger-than-a-hemisphere-polygons-and-the-reorientobject-method-in-sql-server-2012/ <https://alastaira.wordpress.com/2012/01/27/ring-orientation-bigger-than-a-hemisphere-polygons-and-the-reorientobject-method-in-sql-server-2012/>

The second part of the blog post is instructive, in that it shows that “fixing” the problem will just migrate a new problem out to users with small, but mis-oriented polygons. What they think of as local objects (parcels, cities, whatever) will, to their surprise, suddenly become nearly-global.

There’s some interesting stuff in the SQL Server documentation, in particular the idea of using circular envelopes for geographic objects, which has a certain appeal, but would necessitate a pretty major re-write of geography.  

There’s actually two “hemisphere problems” to deal with in any geodetic database: 

- what is the correct interpretation of an edge from A to B, since there are two paths between A and B over the great circle that joins them
- what is the correct interpretation of a ring

The PostGIS geodetic attempts to use “take the smaller thing”, but imperfectly, as you have demonstrated. The question is whether to continue to patch that up, or attempt to reimplement with something more explicit, like “ring orientation determines enclosure” and accept all the other knock-on issues that arise as people try and load mis-oriented rings into the system.

It’s a very interesting and challenging topic, thanks for investigating, I particularly like your point-field method of visualizing the issues with the algorithm, which are associated with the point-in-polygon routine. Point-in-polygon is one of those algorithms that is challenging on the sphere, and more so when what is “in” the polygon is not necessarily clear.

ATB,

P

> On Aug 9, 2019, at 9:53 AM, Christian Pschierer <christian.pschierer at gmx.net> wrote:
> 
> Hi Darafei,
> 
> yes, the MakeBox example was not a good one. Even though the basic question -- which rule determines inside and outside of a large geography polygon -- remains the same. 
> 
> But the problem is a different one: 
> SELECT ST_Distance(ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452,
> 130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105,
> 172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745,
> -36.5625 40.9798980696201))')::geography  , ST_Point(0, 0)) / 1000
> return 0 in PostGis 2.5.2. (as expected) 
> 
> However: ST_Intersect or ST_DWithin on the same geometries return false. Distance 0 and Intersects=false does not make sense. 
> The culprit is the bounding-box check in the latter two function ( OPERATOR(public.&&) ). 
> 
> Test Case: 
> Create a sample dataset with one point per 1x1 degree: 
> CREATE MATERIALIZED VIEW public.tmp_points_1x1 AS ( 
> SELECT row_number() over() AS eid, ST_Translate(point, j, i) AS geom 
> FROM 
> generate_series(-89, 89) AS i, 
> generate_series(-180, 179) AS j, 
> (SELECT ('POINT(0 0)')::geometry AS point) AS b )
> 
> Now select from this dataset: 
> SELECT eid,ST_SetSRID(geom, 4326) FROM tmp_points_1x1 
> WHERE ST_Distance(geom::geography, ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452, 
> 130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105, 
> 172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745, 
> -36.5625 40.9798980696201))')::geography ) <= 0.0
> 
> returns the expected result: 
> <image.png> 
> 
> ST_Intersects leaves a hole around the point 0, 0! 
> SELECT eid,ST_SetSRID(geom::geography, 4326) FROM tmp_points_1x1 
> WHERE ST_Intersects(geom::geography, ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452, 
> 130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105, 
> 172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745, 
> -36.5625 40.9798980696201))')::geography )
> <image.png> 
> 
> && leaves 2 holes and select a lot of additional points outside of the expected bounding box. 
> SELECT eid,ST_SetSRID(geom, 4326) FROM tmp_points_1x1 
> WHERE geom::geography OPERATOR(public.&&) ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452, 
> 130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105, 
> 172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745, 
> -36.5625 40.9798980696201))')::geography
> <image.png> 
> 
> Greetings 
> Christian 
> 
>> On 09 August 2019 at 14:56 "Darafei \"Komяpa\" Praliaskouski" < me at komzpa.net <mailto:me at komzpa.net>> wrote:
>> 
>> 
>> Hi Christian,
>> 
>> To understand the query please try plotting your polygon on a globe and
>> then imagine axis aligned box that will contain it (one side parallel to
>> equator, other parallel to 0 and third to 90th meridian). You will see that
>> your "large" geometry is really 40 degrees wide spot around antimeridian.
>> To visualize it better using common planar GIS tools, feed it into
>> geography ST_Segmentize - it will produce the actual 'curved' geometry used
>> in calculation. To see closer to behavior you expect, make sure there are
>> no points connected by longer, not shorter, part of Great Circle, as in
>> your example.
>> 
>> Hope this helps.
>> 
>> On Thu, Aug 8, 2019 at 2:46 AM Christian Pschierer <
>> christian.pschierer at gmx.net <mailto:christian.pschierer at gmx.net>> wrote:
>> 
>>> Hi,
>>> 
>>> we found some unexpected results when doing spatial queries on very
>>> large geography polygons. For example
>>> 
>>> SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(160,
>>> 60),ST_Point(-160,-60)), 4326)::geography,ST_SetSRID(ST_Point(0, 0),
>>> 4326)::geography)/1000
>>> 
>>> returns 13130km instead of 0 as the point 0,0 should be inside this polygon.
>>> 
>>> Queries on smaller search areas, or geometry queries return 0.
>>> 
>>> SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(60,
>>> 60),ST_Point(-60,-60)), 4326)::geography,ST_SetSRID(ST_Point(0,
>>> 0),4326)::geography)/1000
>>> 
>>> SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(160,
>>> 60),ST_Point(-160,-60)), 4326),ST_SetSRID(ST_Point(0, 0), 4326))/1000
>>> 
>>> It seems like PostGis switches inside/outside on geographies if they
>>> exceed a certain size. Is this correct? Is there a way to control this
>>> behaviour?
>>> 
>>> Greetings
>>> Christian
>>> 
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
>>> https://lists.osgeo.org/mailman/listinfo/postgis-users <https://lists.osgeo.org/mailman/listinfo/postgis-users> 
>> 
>> 
>> 
>> --
>> Darafei Praliaskouski
>> Support me: http://patreon.com/komzpa <http://patreon.com/komzpa> 
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
>> https://lists.osgeo.org/mailman/listinfo/postgis-users <https://lists.osgeo.org/mailman/listinfo/postgis-users> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190812/ab85969b/attachment.html>