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

[Szymon Haas]

 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20190813/ea5336e1/attachment.html>