<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">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.<div class=""><br class=""></div><div class=""> WITH data AS (<br class=""> SELECT <br class=""> '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,<br class=""> 'SRID=4326;POINT(0 0)'::geography AS pt<br class=""> )<br class=""> SELECT <br class=""> ST_Distance(ply,pt) AS distance,<br class=""> ST_DWithin(ply,pt,0) AS dwithin_0,<br class=""> ST_DWithin(ply,pt,1) AS dwithin_1,<br class=""> _ST_DistanceTree(ply,pt) AS distance_tree,<br class=""> _ST_DistanceUnCached(ply,pt) AS distance_uncached<br class=""> FROM data;<br class=""><br class="">WITH data AS (<br class=""> SELECT <br class=""> 'SRID=4326;POLYGON((-40.0 52.0, -67.0 -29.0, 102.0 -6.0, -40.0 52.0))'::geography AS ply,<br class=""> 'SRID=4326;POINT(4 11)'::geography AS pt<br class="">)<br class="">SELECT <br class=""> ST_Distance(ply,pt) AS distance,<br class=""> ST_DWithin(ply,pt,0) AS dwithin_0,<br class=""> ST_DWithin(ply,pt,1) AS dwithin_1,<br class=""> _ST_DistanceTree(ply,pt) AS distance_tree,<br class=""> _ST_DistanceUnCached(ply,pt) AS distance_uncached<br class="">FROM data;<br class=""><br class=""><div>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.</div><div><br class=""></div><div>Lots of thinking left to do,</div><div><br class=""></div><div>P</div><div><br class=""></div><div><br class=""></div><div><br class=""><blockquote type="cite" class=""><div class="">On Aug 12, 2019, at 9:55 AM, Paul Ramsey <<a href="mailto:pramsey@cleverelephant.ca" class="">pramsey@cleverelephant.ca</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><meta http-equiv="Content-Type" content="text/html; charset=utf-8" class=""><div style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">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.<div class=""><br class=""></div><div class=""><a href="https://alastaira.wordpress.com/2012/01/27/ring-orientation-bigger-than-a-hemisphere-polygons-and-the-reorientobject-method-in-sql-server-2012/" class="">https://alastaira.wordpress.com/2012/01/27/ring-orientation-bigger-than-a-hemisphere-polygons-and-the-reorientobject-method-in-sql-server-2012/</a></div><div class=""><br class=""></div><div class="">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.</div><div class=""><br class=""></div><div class="">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. </div><div class=""><br class=""></div><div class="">There’s actually two “hemisphere problems” to deal with in any geodetic database: </div><div class=""><br class=""></div><div class="">- 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</div><div class="">- what is the correct interpretation of a ring</div><div class=""><br class=""></div><div class="">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.</div><div class=""><br class=""></div><div class="">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.</div><div class=""><br class=""></div><div class="">ATB,</div><div class=""><br class=""></div><div class="">P</div><div class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On Aug 9, 2019, at 9:53 AM, Christian Pschierer <<a href="mailto:christian.pschierer@gmx.net" class="">christian.pschierer@gmx.net</a>> wrote:</div><br class="Apple-interchange-newline"><div class="">
<meta charset="UTF-8" class="">
<div class="">
<div class="">
Hi Darafei,
</div>
<div class="">
<br class="">
</div>
<div class="">
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.
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
But the problem is a different one:
<br class="">
</div>
<div class="">
SELECT ST_Distance(ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452,
<div class="">
130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105,
</div>
<div class="">
172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745,
</div>-36.5625 40.9798980696201))')::geography , ST_Point(0, 0)) / 1000
</div>
<div class="">
return 0 in PostGis 2.5.2. (as expected)
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
However: ST_Intersect or ST_DWithin on the same geometries return false. Distance 0 and Intersects=false does not make sense.
<br class="">
</div>
<div class="">
The culprit is the bounding-box check in the latter two function ( OPERATOR(public.&&) ).
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
Test Case:
<br class="">
</div>
<div class="">
Create a sample dataset with one point per 1x1 degree:
<br class="">
</div>
<div class="">
CREATE MATERIALIZED VIEW public.tmp_points_1x1 AS (
<br class="">SELECT row_number() over() AS eid, ST_Translate(point, j, i) AS geom
<br class="">FROM
<br class="">generate_series(-89, 89) AS i,
<br class="">generate_series(-180, 179) AS j,
<br class="">(SELECT ('POINT(0 0)')::geometry AS point) AS b )
</div>
<div class="">
<br class="">
</div>
<div class="">
Now select from this dataset:
<br class="">
</div>
<div class="">
SELECT eid,ST_SetSRID(geom, 4326) FROM tmp_points_1x1
<br class="">WHERE ST_Distance(geom::geography, ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452,
<br class="">130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105,
<br class="">172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745,
<br class="">-36.5625 40.9798980696201))')::geography ) <= 0.0
</div>
<div class="">
<br class="">
</div>
<div class="">
returns the expected result:
<br class="">
</div>
<div class="">
<span id="cid:b77ea88d4a8040e2acba71977eef0c3c@Open-Xchange" class=""><image.png></span>
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
ST_Intersects leaves a hole around the point 0, 0!
<br class="">
</div>
<div class="">
SELECT eid,ST_SetSRID(geom::geography, 4326) FROM tmp_points_1x1
<br class="">WHERE ST_Intersects(geom::geography, ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452,
<br class="">130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105,
<br class="">172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745,
<br class="">-36.5625 40.9798980696201))')::geography )
</div>
<div class="">
<span id="cid:f4d33ea06e63493bbaed4f5c326815e9@Open-Xchange" class=""><image.png></span>
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
&& leaves 2 holes and select a lot of additional points outside of the expected bounding box.
<br class="">
</div>
<div class="">
SELECT eid,ST_SetSRID(geom, 4326) FROM tmp_points_1x1
<br class="">WHERE geom::geography OPERATOR(public.&&) ST_GeomFromText('POLYGON((-36.5625 40.9798980696201, -22.5 -7.71099165543322,6.328125 -44.0875850282452,
<br class="">130.078125 -49.3823727870096,170.546875 -47.9899216674142,170.546875 69.162557908105,
<br class="">172.96875 75.320025232208,68.203125 77.9156689863258,4.21875 71.3007929163745,
<br class="">-36.5625 40.9798980696201))')::geography
</div>
<div class="">
<span id="cid:aff01ee403f64389986eb83bd95998a5@Open-Xchange" class=""><image.png></span>
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
Greetings
<br class="">
</div>
<div class="">
Christian
<br class="">
</div>
<div class="">
<br class="">
</div>
<blockquote type="cite" class="">
<div class="">
On 09 August 2019 at 14:56 "Darafei \"Komяpa\" Praliaskouski" <
<a href="mailto:me@komzpa.net" class="">me@komzpa.net</a>> wrote:
</div>
<div class="">
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
Hi Christian,
</div>
<div class="">
<br class="">
</div>
<div class="">
To understand the query please try plotting your polygon on a globe and
</div>
<div class="">
then imagine axis aligned box that will contain it (one side parallel to
</div>
<div class="">
equator, other parallel to 0 and third to 90th meridian). You will see that
</div>
<div class="">
your "large" geometry is really 40 degrees wide spot around antimeridian.
</div>
<div class="">
To visualize it better using common planar GIS tools, feed it into
</div>
<div class="">
geography ST_Segmentize - it will produce the actual 'curved' geometry used
</div>
<div class="">
in calculation. To see closer to behavior you expect, make sure there are
</div>
<div class="">
no points connected by longer, not shorter, part of Great Circle, as in
</div>
<div class="">
your example.
</div>
<div class="">
<br class="">
</div>
<div class="">
Hope this helps.
</div>
<div class="">
<br class="">
</div>
<div class="">
On Thu, Aug 8, 2019 at 2:46 AM Christian Pschierer <
</div>
<div class="">
<a href="mailto:christian.pschierer@gmx.net" class="">christian.pschierer@gmx.net</a>> wrote:
</div>
<div class="">
<br class="">
</div>
<blockquote type="cite" class="">
<div class="">
Hi,
</div>
<div class="">
<br class="">
</div>
<div class="">
we found some unexpected results when doing spatial queries on very
</div>
<div class="">
large geography polygons. For example
</div>
<div class="">
<br class="">
</div>
<div class="">
SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(160,
</div>
<div class="">
60),ST_Point(-160,-60)), 4326)::geography,ST_SetSRID(ST_Point(0, 0),
</div>
<div class="">
4326)::geography)/1000
</div>
<div class="">
<br class="">
</div>
<div class="">
returns 13130km instead of 0 as the point 0,0 should be inside this polygon.
</div>
<div class="">
<br class="">
</div>
<div class="">
Queries on smaller search areas, or geometry queries return 0.
</div>
<div class="">
<br class="">
</div>
<div class="">
SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(60,
</div>
<div class="">
60),ST_Point(-60,-60)), 4326)::geography,ST_SetSRID(ST_Point(0,
</div>
<div class="">
0),4326)::geography)/1000
</div>
<div class="">
<br class="">
</div>
<div class="">
SELECT ST_Distance(ST_SetSRID( ST_MakeBox2D(ST_Point(160,
</div>
<div class="">
60),ST_Point(-160,-60)), 4326),ST_SetSRID(ST_Point(0, 0), 4326))/1000
</div>
<div class="">
<br class="">
</div>
<div class="">
It seems like PostGis switches inside/outside on geographies if they
</div>
<div class="">
exceed a certain size. Is this correct? Is there a way to control this
</div>
<div class="">
behaviour?
</div>
<div class="">
<br class="">
</div>
<div class="">
Greetings
</div>
<div class="">
Christian
</div>
<div class="">
<br class="">
</div>
<div class="">
_______________________________________________
</div>
<div class="">
postgis-users mailing list
</div>
<div class="">
<a href="mailto:postgis-users@lists.osgeo.org" class="">postgis-users@lists.osgeo.org</a>
<br class="">
</div>
<div class="">
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noopener" target="_blank" class="">https://lists.osgeo.org/mailman/listinfo/postgis-users</a>
<br class="">
</div>
</blockquote>
<div class="">
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
<br class="">
</div>
<div class="">
--
</div>
<div class="">
Darafei Praliaskouski
</div>
<div class="">
Support me:
<a href="http://patreon.com/komzpa" rel="noopener" target="_blank" class="">http://patreon.com/komzpa</a>
<br class="">
</div>
<div class="">
_______________________________________________
</div>
<div class="">
postgis-users mailing list
</div>
<div class="">
<a href="mailto:postgis-users@lists.osgeo.org" class="">postgis-users@lists.osgeo.org</a>
<br class="">
</div>
<div class="">
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noopener" target="_blank" class="">https://lists.osgeo.org/mailman/listinfo/postgis-users</a>
<br class="">
</div>
</blockquote>
</div>
_______________________________________________<br class="">postgis-users mailing list<br class=""><a href="mailto:postgis-users@lists.osgeo.org" class="">postgis-users@lists.osgeo.org</a><br class=""><a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" class="">https://lists.osgeo.org/mailman/listinfo/postgis-users</a></div></blockquote></div><br class=""></div></div></div></blockquote></div><br class=""></div></body></html>