[postgis-tickets] [PostGIS] #3825: Running an ST_Intersect with PostGIS Geography with large extents returns incorrect results.

[PostGIS]

#3825: Running an ST_Intersect with PostGIS Geography with large extents returns
incorrect results.
--------------------------------+---------------------------
 Reporter:  just7460            |      Owner:  pramsey
     Type:  defect              |     Status:  new
 Priority:  medium              |  Milestone:  PostGIS 2.3.4
Component:  postgis             |    Version:  2.3.x
 Keywords:  Geography, postgis  |
--------------------------------+---------------------------
 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';
 ====

--
Ticket URL: <https://trac.osgeo.org/postgis/ticket/3825>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.