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
_______________________________________________ geos-devel mailing list geos-devel@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/geos-devel