ogr2ogr: fix -clipsrc/-clipdst when a input geometry is within the en…

…velope of the clip geometry, but doesn't intersect it (3.10.0 regression) Fixes #10341 This essentially reverts the gist of commit 388d3ae (PR #10341) "ogr2ogr: speed-up -clipsrc/-clipdst by avoiding GEOS when possible", which based on a totally wrong assumption.
OSGeo · Jan 14, 2025 · 0f3c39b · 0f3c39b
1 parent c7e2b4f
commit 0f3c39b
Show file tree

Hide file tree

Showing 2 changed files with 101 additions and 60 deletions.
diff --git 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<OGRGeometry> poClipped;
+                    if (poClipGeomEnvelope->Intersects(oDstEnv))
                     {
-                        std::unique_ptr<OGRGeometry> 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<OGRGeometry> poClipped;
+                        if (poClipGeomEnvelope->Intersects(oDstEnv))
                         {
-                            std::unique_ptr<OGRGeometry> 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
@@ -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)