> On Sep 17, 2020, at 8:50 AM, Martin Davis <mtncl...@gmail.com> wrote: > > Thanks for testing, Even. > > Interestingly, in JTS the test case you give below for OverlayNG.intersection > actually returns the same result as the old code: > > MULTILINESTRING ((833978.557030887 0, 835092.849076364 0), (832864.275023695 > 0, 833978.556808034 0)) > > So this is something to look into on the GEOS side.
I have confirmed this happens in GEOS, with a little unit test, so it's nothing to do with GDAL. Looks like a morning in the debugger awaits. P > > On Thu, Sep 17, 2020 at 4:39 AM Even Rouault <even.roua...@spatialys.com> > wrote: > > > 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)) > _______________________________________________ > geos-devel mailing list > geos-devel@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/geos-devel _______________________________________________ geos-devel mailing list geos-devel@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/geos-devel