[geos-devel] GEOS NG Regression Review (GDAL feedback)

[Even Rouault]

Hi Paul,

I gave a try to GEOS git head built with -DDISABLE_OVERLAYNG=OFF, and ran it against 
GDAL autotest suite.

The needed changes are in https://github.com/OSGeo/gdal/pull/2942

Most of them are boring, due to polygon vertices not being returned in the same order, and 
fixed by using spatial equality instead of plain WKT comparison.

One non-trivial change I had to do was related to a convoluted way of how OGR tries to split 
line or polygon geometries crossing the antimeridian, when reprojecting from a projected 
CRS crossing the antimedian (e.g. UTM zone 60), to geographic coordinates. The resulting 
geometry needs to be a multilinestring or multipolygon with one (or several) part close to 
longitude +180 and another one close to longitude -180. Typically needed when reprojecting 
geometries for GeoJSON RFC7946 compliance.

1) For each segment of the input geometry that crosses the antimeridian, OGR determines 
by dichotomy the coordinates, in the projected CRS, of the intersection with the 
antimeridian. This leads to a list of points P_antimeridian[], that is sorted by increasing y.
2) Previously, OGR built a polygon, at the left of the antimeridian, whose right edge are 
points whose x is at P_antimeridian[].x - epsilon, and another one, at the right of the 
antimeridian, whose left edge are points whose x is at P_antimeridian[].x + epsilon. Then it 
creates a multipolygon with those 2 polygons.
3) Then it intersects the input geometry with those 2 multipolygons, which leads to a multi-
geometry not crossing the antimeridian.
4) And finally it reprojects that geometry to geographic coordinates.

The intersection of the input geometry with this multipolygon of 2 parts can be shown with:

from osgeo import ogr

# input geometry crossing the antimeridian (UTM 60N)
geom = ogr.CreateGeometryFromWkt('LINESTRING(832864.275023695 0,835092.849076364 
0)')

# multipolygon with one part left to the antimeridian, one part right
geom2 = ogr.CreateGeometryFromWkt('MULTIPOLYGON (((832864.275023695 
0.0,833978.556808034 -0.000110682755987,833978.556808034 0.0,833978.556808034 
0.000110682755987,832864.275023695 0.0,832864.275023695 0.0)),((835092.849076364 
0.0,833978.557030887 -0.000110682755987,833978.557030887 0.0,833978.557030887 
0.000110682755987,835092.849076364 0.0,835092.849076364 0.0)))')

# intersection
print(geom.Intersection(geom2))

With OverlayNG, the following leads to a LINESTRING EMPTY, whereas with GEOS 
3.8something, it leads to the expected result of
MULTILINESTRING ((832864.275023695 0.0,833978.556808034 0.0),(833978.557030887 
0.0,835092.849076364 0.0))

I'm not sure if the OverlayNG result is a feature or a bug.

Anyway, I've changed the GDAL code to compute the difference between the input 
geometry and a polygon with one edge with points whose x is at P_antimeridian[].x - epsilon 
and another one with points whose x is at P_antimeridian[].x + epsilon. This is more 
straightforward, and works with any GEOS version.

(by the way, I don't think the GDAL autotest suite exercices much GEOS Z handling)

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20200917/02c47f09/attachment.html>