[postgis-users] Results of ST_DWithin missing data around 0, 0 for large request areas

[Paul Ramsey]

You may find that this problem is recently addressed (how does this bug exist for 10 years and get reported twice in two weeks?) in https://trac.osgeo.org/postgis/ticket/4480 <https://trac.osgeo.org/postgis/ticket/4480> 

While in general the rule “things that are larger than a hemisphere are going to be problematic” does hold, the reported case was in fact addressable and fixable, if not in general then at least for most common cases.

In general:

- try to break your large things into equivalent sets of smaller things
- remember that no edge will span more than 180 degrees, as we do interpret an edge as the shortest distance between a coordinate pair

ATB,

P

> On Aug 13, 2019, at 9:55 AM, Szymon Haas <haasel77 at gmail.com> wrote:
> 
> Dear PostGIS users,
> 
> We are using PostGIS as an engine to retrieve data in GeoJSON format.
> As the data is worldwide and requests concern large geographic areas we did a choice of using "geography" data type.
> 
> The main functionality of our system is retrieving data in some "places" (points, lines, areas), including buffers. 
> Thus, after several tries, we decided to use ST_DWithin function.
> 
> Recently, during the tests, we found an issue with that function.
> Briefly - for very large areas (from the observations BBOX wider than 180 degrees EW) results of the ST_DWithin are missing data in a circular area 
> around the 0,0 point (the prime meridian and equator cross-section).
> 
> To reproduce the issue one can create a table/view with a grid of points (lines or polygons would have the same effect) around the globe:
> 
> CREATE MATERIALIZED VIEW tmp_points_1x1 AS (
> SELECT row_number() over() AS eid, ST_Translate(point, j, i)::geography AS geog
>     FROM 
>       generate_series(-89, 89) AS i,
>       generate_series(-180, 179) AS j,
>       (SELECT ('POINT(0 0)')::geometry AS point) AS b )
> 
> adding these indexes will speedup test queries:
> CREATE INDEX tmp_points_1x1_geog
>     ON tmp_points_1x1 USING gist(geog)
>     TABLESPACE pg_default;
> 
> CREATE INDEX tmp_points_1x1_eid
>     ON tmp_points_1x1 USING btree(eid)
>     TABLESPACE pg_default;
>     
> The below test query returns in Json format data for requested geojson area (here it's 190 degrees wide):
> SELECT row_to_json(fc)
> FROM (
>     SELECT 'FeatureCollection' AS type, array_to_json(array_agg(f)) AS features
>     FROM (
>         SELECT 'Feature' AS type, 
>         ST_AsGeoJSON(lg.geog)::json AS geometry, 
>         row_to_json(lp) AS properties
>         FROM tmp_points_1x1 AS lg 
>         INNER JOIN (SELECT eid 
>                     FROM tmp_points_1x1 
>                     WHERE st_dwithin(geog, ST_GEOMFROMGEOJSON('{"type":"Polygon","coordinates":[[[-20,-50],[75,-50],[170,-50],[170,50],[75,50],[-20,50],[-20,-50]]]}' ),0)
>                     ) AS lp 
>         ON lg.eid = lp.eid  
>         ) AS f 
>     ) AS fc;
> 
> On the visualization of the output you can see the missing data:
>  above180wide.png <https://drive.google.com/file/d/1wgkpLi4xyCviN3AzE8hAnRZUFVUbH3ST/view?usp=drive_web>
> If the request area is narrowed to <180 degrees wide, e.g. 170 degrees: {"type":"Polygon","coordinates":[[[-20,-50],[20,-50],[20,50],[-20,50],[-20,-50]]]}
> the problem of missing data doesn't appear any more:
> 
>  below180wide.png <https://drive.google.com/file/d/1FuXOmESE6N4_LWAjqAxg0WbOAjwJLeEe/view?usp=drive_web>Is it a known issue?
> Are there chances for that to be resolved?
> Or maybe it's not a bug (e.g. request geometries' BBOXes should not exceed 180 degrees?)
> 
> Thank you for any explanation,
> Eliasz Haas
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190813/2c22f630/attachment.html>