From 0f3c39bd990ed5627e9066b50593cb7b2c49d8ef Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 15 Jan 2025 00:57:47 +0100 Subject: [PATCH] ogr2ogr: fix -clipsrc/-clipdst when a input geometry is within the envelope of the clip geometry, but doesn't intersect it (3.10.0 regression) Fixes #10341 This essentially reverts the gist of commit 388d3ae7e47cba64e1050023bed43c4b04fe4277 (PR #10341) "ogr2ogr: speed-up -clipsrc/-clipdst by avoiding GEOS when possible", which based on a totally wrong assumption. --- apps/ogr2ogr_lib.cpp | 113 ++++++++++++------------- autotest/utilities/test_ogr2ogr_lib.py | 48 +++++++++++ 2 files changed, 101 insertions(+), 60 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index ea3e52bba62e..febd87520dfa 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -6753,39 +6753,36 @@ bool LayerTranslator::Translate( OGREnvelope oDstEnv; poDstGeometry->getEnvelope(&oDstEnv); - if (!poClipGeomEnvelope->Contains(oDstEnv)) + std::unique_ptr poClipped; + if (poClipGeomEnvelope->Intersects(oDstEnv)) { - std::unique_ptr poClipped; - if (poClipGeomEnvelope->Intersects(oDstEnv)) - { - poClipped.reset( - poClipGeom->Intersection(poDstGeometry.get())); - } - if (poClipped == nullptr || poClipped->IsEmpty()) - { - goto end_loop; - } + poClipped.reset( + poClipGeom->Intersection(poDstGeometry.get())); + } - const int nDim = poDstGeometry->getDimension(); - if (poClipped->getDimension() < nDim && - wkbFlatten(poDstFDefn->GetGeomFieldDefn(iGeom) - ->GetType()) != wkbUnknown) - { - CPLDebug( - "OGR2OGR", - "Discarding feature " CPL_FRMT_GIB - " of layer %s, " - "as its intersection with -clipsrc is a %s " - "whereas the input is a %s", - nSrcFID, poSrcLayer->GetName(), - OGRToOGCGeomType(poClipped->getGeometryType()), - OGRToOGCGeomType( - poDstGeometry->getGeometryType())); - goto end_loop; - } + if (poClipped == nullptr || poClipped->IsEmpty()) + { + goto end_loop; + } - poDstGeometry = std::move(poClipped); + const int nDim = poDstGeometry->getDimension(); + if (poClipped->getDimension() < nDim && + wkbFlatten( + poDstFDefn->GetGeomFieldDefn(iGeom)->GetType()) != + wkbUnknown) + { + CPLDebug( + "OGR2OGR", + "Discarding feature " CPL_FRMT_GIB " of layer %s, " + "as its intersection with -clipsrc is a %s " + "whereas the input is a %s", + nSrcFID, poSrcLayer->GetName(), + OGRToOGCGeomType(poClipped->getGeometryType()), + OGRToOGCGeomType(poDstGeometry->getGeometryType())); + goto end_loop; } + + poDstGeometry = std::move(poClipped); } OGRCoordinateTransformation *const poCT = @@ -6936,41 +6933,37 @@ bool LayerTranslator::Translate( OGREnvelope oDstEnv; poDstGeometry->getEnvelope(&oDstEnv); - if (!poClipGeomEnvelope->Contains(oDstEnv)) + std::unique_ptr poClipped; + if (poClipGeomEnvelope->Intersects(oDstEnv)) { - std::unique_ptr poClipped; - if (poClipGeomEnvelope->Intersects(oDstEnv)) - { - poClipped.reset(poClipGeom->Intersection( - poDstGeometry.get())); - } - - if (poClipped == nullptr || poClipped->IsEmpty()) - { - goto end_loop; - } + poClipped.reset( + poClipGeom->Intersection(poDstGeometry.get())); + } - const int nDim = poDstGeometry->getDimension(); - if (poClipped->getDimension() < nDim && - wkbFlatten(poDstFDefn->GetGeomFieldDefn(iGeom) - ->GetType()) != wkbUnknown) - { - CPLDebug( - "OGR2OGR", - "Discarding feature " CPL_FRMT_GIB - " of layer %s, " - "as its intersection with -clipdst is a %s " - "whereas the input is a %s", - nSrcFID, poSrcLayer->GetName(), - OGRToOGCGeomType( - poClipped->getGeometryType()), - OGRToOGCGeomType( - poDstGeometry->getGeometryType())); - goto end_loop; - } + if (poClipped == nullptr || poClipped->IsEmpty()) + { + goto end_loop; + } - poDstGeometry = std::move(poClipped); + const int nDim = poDstGeometry->getDimension(); + if (poClipped->getDimension() < nDim && + wkbFlatten(poDstFDefn->GetGeomFieldDefn(iGeom) + ->GetType()) != wkbUnknown) + { + CPLDebug( + "OGR2OGR", + "Discarding feature " CPL_FRMT_GIB + " of layer %s, " + "as its intersection with -clipdst is a %s " + "whereas the input is a %s", + nSrcFID, poSrcLayer->GetName(), + OGRToOGCGeomType(poClipped->getGeometryType()), + OGRToOGCGeomType( + poDstGeometry->getGeometryType())); + goto end_loop; } + + poDstGeometry = std::move(poClipped); } if (psOptions->dfXYRes != diff --git a/autotest/utilities/test_ogr2ogr_lib.py b/autotest/utilities/test_ogr2ogr_lib.py index 1db27aa165ac..bb7fcd43cfc0 100755 --- a/autotest/utilities/test_ogr2ogr_lib.py +++ b/autotest/utilities/test_ogr2ogr_lib.py @@ -1115,6 +1115,17 @@ def test_ogr2ogr_lib_clipsrc_datasource(tmp_vsimem): f.SetField("filter_field", "exact_overlap_full_result") f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))")) clip_layer.CreateFeature(f) + # Clip geometry envelope contains envelope of input geometry, but does not intersect it + f = ogr.Feature(clip_layer.GetLayerDefn()) + f.SetField( + "filter_field", "clip_geometry_envelope_contains_envelope_but_no_intersect" + ) + f.SetGeometry( + ogr.CreateGeometryFromWkt( + "POLYGON ((-2 -1,-2 4,4 4,4 -1,3 -1,3 3,-1 3,-1 -1,-2 -1))" + ) + ) + clip_layer.CreateFeature(f) clip_ds = None # Test clip with 'half_overlap_line_result' using sql statement @@ -1161,6 +1172,19 @@ def test_ogr2ogr_lib_clipsrc_datasource(tmp_vsimem): dst_ds = None gdal.Unlink(dst_filename) + # Test clip with the "clip_geometry_envelope_contains_envelope_but_no_intersect" using only clipSrcWhere + dst_ds = gdal.VectorTranslate( + dst_filename, + src_filename, + format="GPKG", + clipSrc=clip_path, + clipSrcWhere="filter_field = 'clip_geometry_envelope_contains_envelope_but_no_intersect'", + ) + dst_lyr = dst_ds.GetLayer(0) + assert dst_lyr.GetFeatureCount() == 0 + dst_ds = None + gdal.Unlink(dst_filename) + # Cleanup gdal.Unlink(clip_path) @@ -1389,6 +1413,17 @@ def test_ogr2ogr_lib_clipdst_datasource(tmp_vsimem): f.SetField("filter_field", "exact_overlap_full_result") f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))")) clip_layer.CreateFeature(f) + # Clip geometry envelope contains envelope of input geometry, but does not intersect it + f = ogr.Feature(clip_layer.GetLayerDefn()) + f.SetField( + "filter_field", "clip_geometry_envelope_contains_envelope_but_no_intersect" + ) + f.SetGeometry( + ogr.CreateGeometryFromWkt( + "POLYGON ((-2 -1,-2 4,4 4,4 -1,3 -1,3 3,-1 3,-1 -1,-2 -1))" + ) + ) + clip_layer.CreateFeature(f) clip_ds = None # Test clip with 'half_overlap_line_result' using sql statement @@ -1435,6 +1470,19 @@ def test_ogr2ogr_lib_clipdst_datasource(tmp_vsimem): dst_ds = None gdal.Unlink(dst_filename) + # Test clip with the "clip_geometry_envelope_contains_envelope_but_no_intersect" using only clipSrcWhere + dst_ds = gdal.VectorTranslate( + dst_filename, + src_filename, + format="GPKG", + clipDst=clip_path, + clipDstWhere="filter_field = 'clip_geometry_envelope_contains_envelope_but_no_intersect'", + ) + dst_lyr = dst_ds.GetLayer(0) + assert dst_lyr.GetFeatureCount() == 0 + dst_ds = None + gdal.Unlink(dst_filename) + # Cleanup gdal.Unlink(clip_path)