[postgis-users] Inside / Outside of geography polygons
Christian Pschierer
christian.pschierer at gmx.net
Wed Aug 14 07:48:27 PDT 2019
Thanks Paul for the quick reply and the patch!
Regarding reimplementation with something more explicit, like “ring orientation determines enclosure”.
I would prefer such a solution. From a users perspective, some polygons are larger than a hemisphere. E.g. great circle mapper: "show me a polygon with 15000km radius around London" (http://www.gcmap.com/mapui?R=15000km@LHR ). This is different from "5000km around a point somewhere in the south Pacific". As a power user I want to be able to force one or the other interpretation.
To mitigate the knock-on effects, a convenience function might be helpful to automatically reverse the polygon if the area is >>255M square km.
I can't speak for the implementation complexity though. Especially if such a solution could open new issues, or break existing applications.
Greetings
Christian
> On 13 August 2019 at 00:48 Paul Ramsey <pramsey at cleverelephant.ca> wrote:
>
>
> So, I have a patch for this, against 3.0, it’s a little big, but can probably also apply against 2.5 and 2.4… it’s not great, it’s just another layer of plaster on top of a shakey foundation. An error in the computation of the XYZ bounding volume is fixed, which exposes more cases where an external point cannot be generated, so a hack to plaster over that issue is also added. At the end of it, these two test cases now return the right answers, though I have not today yet run your clever point field tests to ensure that things look more correct over the whole point field.
>
> WITH data AS (
> SELECT
> 'SRID=4326;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 AS ply,
> 'SRID=4326;POINT(0 0)'::geography AS pt
> )
> SELECT
> ST_Distance(ply,pt) AS distance,
> ST_DWithin(ply,pt,0) AS dwithin_0,
> ST_DWithin(ply,pt,1) AS dwithin_1,
> _ST_DistanceTree(ply,pt) AS distance_tree,
> _ST_DistanceUnCached(ply,pt) AS distance_uncached
> FROM data;
>
> WITH data AS (
> SELECT
> 'SRID=4326;POLYGON((-40.0 52.0, -67.0 -29.0, 102.0 -6.0, -40.0 52.0))'::geography AS ply,
> 'SRID=4326;POINT(4 11)'::geography AS pt
> )
> SELECT
> ST_Distance(ply,pt) AS distance,
> ST_DWithin(ply,pt,0) AS dwithin_0,
> ST_DWithin(ply,pt,1) AS dwithin_1,
> _ST_DistanceTree(ply,pt) AS distance_tree,
> _ST_DistanceUnCached(ply,pt) AS distance_uncached
> FROM data;
>
> One thing that leaks into this “fix” is an assumption of CCW exterior ring orientation for polygons that fail the usual external point computation phase (generally, large ones). So orientation is starting to rear its head, but only for special cases.
>
> Lots of thinking left to do,
>
> P
>
>
>
> > On Aug 12, 2019, at 9:55 AM, Paul Ramsey <pramsey at cleverelephant.ca> wrote:
> >
> > 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 <mailto: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 <mailto:postgis-users at lists.osgeo.org>
> >> 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
More information about the postgis-users
mailing list