Running an ST_Intersect with PostGIS Geography with large extents returns incorrect results.

[just7460]

SYNOPSIS: When trying to find which features from a table intersect with a large provided extent, only a portion of the expected features are returned. The provided extent is significantly larger than the extent holding the features.

More specifically, the table contains 4 features which are all located in North America. The SQL statements all attempt to find the features which intersect with a shape that extends to roughly the entire globe. Attempting to find the features that intersect with a polygon created by the 4 most extreme points will return 0 results. Because I assumed that PostGIS Geography uses Great Circles, I also included some SQL which attempts to find the intersecting features with a densified polygon. This seems to find some of the features, but are still missing some.

POSTGRESQL SPECIFIC INFORMATION: Version: "PostgreSQL 9.6.3, compiled by Visual C++ build 1800, 64-bit"

Special Installation requirements: N/A

Special startup parameters: N/A

OPERATING SYSTEM INFORMATION: OS: Windows 8.1 Enterprise Processor: Intel® Xeon® CPU E5-1620 v3 @ 3.50GHz System type: 64-bit Operating System, x64-based processor

STEPS TO REPRODUCE: Run the following SQL, or run the attached .sql file via psql.

===== DROP TABLE IF EXISTS geog_intersect;

CREATE TABLE geog_intersect (

oid serial NOT NULL, code VARCHAR(3), geog Geography (Geometry, 4326)

);

INSERT INTO geog_intersect (geog)

VALUES (ST_GeogFromText('SRID=4326; Point (-110.45 57.54)'));

INSERT INTO geog_intersect (geog)

VALUES (ST_GeogFromText('SRID=4326; Point ( -98.57 54.63)'));

INSERT INTO geog_intersect (geog)

VALUES (ST_GeogFromText('SRID=4326; Point ( -89.84 52.69)'));

INSERT INTO geog_intersect (geog)

VALUES (ST_GeogFromText('SRID=4326; Point (-101.24 62.14)'));

create index spa_idx on geog_intersect USING GIST (geog);

select oid, ST_AsText(geog) AS geog from geog_intersect

where (geog && ST_GeogFromText('SRID=4326; POLYGON ((

-179.99 -89.99, 179.99 -89.99, 179.99 89.99, -179.99 89.99, -179.99 -89.99))') = 't');

select oid, ST_AsText(geog) AS geog from geog_intersect

where (geog && ST_GeogFromText('SRID=4326; POLYGON ((

-179.99 -89.99, -159.99 -89.99, -139.99 -89.99, -119.99 -89.99, -99.99 -89.99,

-80.00 -89.99, -60.00 -89.99, -40.00 -89.99, -20.00 -89.99, -0.00 -89.99, 20.00 -89.99,

40.00 -89.99, 60.00 -89.99, 80.00 -89.99, 99.99 -89.99, 119.99 -89.99, 139.99 -89.99,

179.99 -89.99, 179.99 -45.00, 179.99 0.00, 179.99 45.00, 179.99 89.99, 159.99 89.99, 139.99 89.99, 119.99 89.99, 99.99 89.99, 80.00 89.99,

60.00 89.99, 40.00 89.99, 20.00 89.99, 0.00 89.99, -20.00 89.99, -40.00 89.99,

-60.00 89.99, -80.00 89.99, -99.99 89.99, -119.99 89.99, -139.99 89.99,

-179.99 89.99, -179.99 45.00, -179.99 0.00, -179.99 -45.00, -179.99 -89.99))') = 't');

select oid, ST_AsText(geog) as geog from geog_intersect where ST_Intersects(geog, ST_GeogFromText('SRID=4326; POLYGON (( -179.99 -89.99, 179.99 -89.99, 179.99 89.99, -179.99 89.99, -179.99 -89.99))')) = 't';

select oid, ST_AsText(geog) as geog from geog_intersect where ST_Intersects(geog, ST_GeogFromText('SRID=4326; POLYGON (( -179.99 -89.99, -159.99 -89.99, -139.99 -89.99, -119.99 -89.99, -99.99 -89.99,

-80.00 -89.99, -60.00 -89.99, -40.00 -89.99, -20.00 -89.99, -0.00 -89.99, 20.00 -89.99,

40.00 -89.99, 60.00 -89.99, 80.00 -89.99, 99.99 -89.99, 119.99 -89.99, 139.99 -89.99,

179.99 -89.99, 179.99 89.99, 159.99 89.99, 139.99 89.99, 119.99 89.99, 99.99 89.99, 80.00 89.99,

60.00 89.99, 40.00 89.99, 20.00 89.99, 0.00 89.99, -20.00 89.99, -40.00 89.99,

-60.00 89.99, -80.00 89.99, -99.99 89.99, -119.99 89.99, -139.99 89.99,

-179.99 89.99, -179.99 -89.99))')) = 't'; ====

[just7460]

.sql file containing the reproducing SQL.

[robe]

I didn't look at your examples, but geography can only work with a half-hemisphere and from your description it sounds like you are trying to do more than that.

​http://postgis.net/docs/manual-2.3/using_postgis_dbmanagement.html#idm1345

As such I'm marking this as wontfix. If you have a more specific example I can follow and are sure you are not falling into the geography flattened illusion, feel free to reopen.

[just7460]

I understand that it has been some time since we discussed this issue, but can you provide with an update, or a location to make suggestions, on support for features larger than a half-hemisphere?

This is a large problem for us as our customers expect to be able to see their data at 'full extent' (global scale). SQL Server originally had issues with features larger than half a hemisphere, but they've resolved this issue since SQL Server 2008.

Looking forward to your response!

[robe]

You could put in another ticket as a feature request, to handle more than half-hemisphere for geography.

pramsey would have to respond how much work that would be, as I presume we would have to also put restrictions on polygon orientation. As I recall I think that's how SQL Server handle's it by assuming polygon draw direction.

That said, Things get a done a lot faster if people are willing to fund the effort. You should express your willingness to fund the work in your ticket.

Thanks, Regina

[just7460]

Thanks for your response! I'll submit another ticket for the feature request.