ST_Intersects polygon geography empty results

[royjacksonalertmedia]

I have a geography dataset containing various (valid) unioned polygons, including the United States (incl Alaska and Hawaii).

I use a materialized view with a spatial index to st_union the shapes for a given id, which include one or many states, countries, circles.

I published this view as a feature service with arcgis enterprise 11.1, and noticed a feature tile was returning no data. Specifically, a large chunk of alaska is missing.

I enabled logging in the database to figure out the query being used:

select objectid,shape from maps_db.public.threat_polygon_shapes_layer where (id='faee6fc1d8f2f583577f85c54da3900b') AND (public.ST_Intersects(shape,public.ST_setSRID(public.ST_GeogFromWKB('\x01030000000100000005000000386094f6ff7f66c000415442d9a05040386094f6ff7f66c0047fa5b145435540fcffffffff7f56c0047fa5b145435540fcffffffff7f56c000415442d9a05040386094f6ff7f66c000415442d9a05040'),'4326')) = 't') ORDER BY objectid ASC

The threat polygon here (alaska + US) is st_valid. The srids are the same. The extent polygon is touching the northern end of the globe.

I can view the two shapes together in pgadmin, and they overlap. I can also cast as geometry and union them together.

If I cast the shapes as geometry in the ST_Intersects query, the query returns expected results.

select objectid,shape from maps_db.public.threat_polygon_shapes_layer where (id='faee6fc1d8f2f583577f85c54da3900b') AND (public.ST_Intersects(shape::geometry,public.ST_setSRID(public.ST_GeogFromWKB('\x01030000000100000005000000386094f6ff7f66c000415442d9a05040386094f6ff7f66c0047fa5b145435540fcffffffff7f56c0047fa5b145435540fcffffffff7f56c000415442d9a05040386094f6ff7f66c000415442d9a05040'),'4326')::geometry) = 't') ORDER BY objectid ASC

The binary representation of the US + Alaska polygon is found in an attached file.

PostgreSQL 14.7 on aarch64-unknown-linux-gnu, compiled by aarch64-unknown-linux-gnu-gcc (GCC) 9.5.0, 64-bit

PostGIS version details:

POSTGIS="3.3.2 0" [EXTENSION] PGSQL="140" GEOS="3.10.3-CAPI-1.16.1" PROJ="9.1.0" LIBXML="2.9.9" LIBJSON="0.12.99" LIBPROTOBUF="1.3.0" WAGYU="0.5.0 (Internal)"

Same results with: POSTGIS="3.3.3 2355e8e" [EXTENSION] PGSQL="140" GEOS="3.9.0-CAPI-1.16.2" PROJ="7.2.1" LIBXML="2.9.10" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"

[robe]

Unfortunately PostGIS has no function to check validity in geography space. I checked in SQL Server which does have an IsValid check for geography and it deemed your alaskan/US boundary to be invalid. However running STMakeValid on it did make it intersect. So I fear this is just a limitation of how we assume shortest length is line in the geography plumbing.

@pramsey you see anything that can be done about this or is it just a hopeless case?

I checked the usual suspect of different distances:

SELECT ST_Distance(shape::geography, geog), _ST_DistanceTree(shape::geography, geog)
, _ST_DistanceUncached(shape::geography, geog) --, shape, geog 
FROM threat_polygon_shapes_layer_valid, public.ST_setSRID(public.ST_GeogFromWKB('\x01030000000100000005000000386094f6ff7f66c000415442d9a05040386094f6ff7f66c0047fa5b145435540fcffffffff7f56c0047fa5b145435540fcffffffff7f56c000415442d9a05040386094f6ff7f66c000415442d9a05040'),'4326') AS f(geog);

 st_distance   | _st_distancetree | _st_distanceuncached
----------------+------------------+----------------------
 32968.11078731 |   32968.11078731 |    32968.11078730672
(1 row)

and it really thinks they don't intersect regardless if I use the validated one or original invalid. I get the same answer using the original geometry as well as the one I had converted to valid in SQL Server using

select objectid,shape.MakeValid().STDistance( geography::STGeomFromWKB(0x01030000000100000005000000386094f6ff7f66c000415442d9a05040386094f6ff7f66c0047fa5b145435540fcffffffff7f56c0047fa5b145435540fcffffffff7f56c000415442d9a05040386094f6ff7f66c000415442d9a05040,'4326') ) from threat_polygon_shapes_layer where (id='faee6fc1d8f2f583577f85c54da3900b') AND shape.MakeValid().STIntersects(geography::STGeomFromWKB(0x01030000000100000005000000386094f6ff7f66c000415442d9a05040386094f6ff7f66c0047fa5b145435540fcffffffff7f56c0047fa5b145435540fcffffffff7f56c000415442d9a05040386094f6ff7f66c000415442d9a05040,'4326') ) = 1 ORDER BY objectid ASC;

SQL Server says they should intersect and distance is 0.

[robe]

geography_made_valid with SQL Server .MakeValid() function

[pramsey]

boxes_arent_boxes.png

[pramsey]

Remember that all lines in geography are great circles, and that the further north you get the more a "horizontal" line in equirectangular projection will in fact arc northward, until in the limit it is going over the pole. This constantly trips people up applying "rectangular" boxes from flat mapping UIs to geography, and I do not see any way to "fix" it that does not do damage to peoples understanding in some other way. We just aren't use to reasoning about spherical geometry in general.

[pramsey]

Oh, if you want to play this game yourself, use the geography version of ST_Segmentize. I used this:

create table uhhuh as 
select st_segmentize(
  'POLYGON((
    -179.99999550799998 66.51326044300004,
    -179.99999550799998 85.05112878000006,
    -89.99999999999994 85.05112878000006,
    -89.99999999999994 66.51326044300004,
    -179.99999550799998 66.51326044300004))'::geography, 
  10000) as geom;

[robe]

I also figured out why SQL Server thinks they intersect. It is seeing the box as covering everything but the box, presumably because it is using the orientation winding of the vertices where as PostGIS is just going by shortest length.

So this is what PostGIS sees: (using the query that Paul showed above)

This is what SQL Server sees: