From 9d52f689f047fecebd20a554e328a7e55603e151 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 a non-rectangle clip geometry, but doesn't intersect it (3.10.0 regression) Fixes #11652 Fixes #10341 This fixes commit 388d3ae7e47cba64e1050023bed43c4b04fe4277 (PR #10341) "ogr2ogr: speed-up -clipsrc/-clipdst by avoiding GEOS when possible", which only worked if the clipping geometry was a rectangle. --- apps/ogr2ogr_lib.cpp | 75 ++++++++++++++++---------- autotest/utilities/test_ogr2ogr_lib.py | 48 +++++++++++++++++ 2 files changed, 96 insertions(+), 27 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 698e77fd99e7..d4e27e8bd884 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -581,12 +581,14 @@ class LayerTranslator std::unique_ptr m_poClipSrcReprojectedToSrcSRS{}; const OGRSpatialReference *m_poClipSrcReprojectedToSrcSRS_SRS = nullptr; OGREnvelope m_oClipSrcEnv{}; + bool m_bClipSrcIsRectangle = false; OGRGeometry *m_poClipDstOri = nullptr; bool m_bWarnedClipDstSRS = false; std::unique_ptr m_poClipDstReprojectedToDstSRS{}; const OGRSpatialReference *m_poClipDstReprojectedToDstSRS_SRS = nullptr; OGREnvelope m_oClipDstEnv{}; + bool m_bClipDstIsRectangle = false; bool m_bExplodeCollections = false; bool m_bNativeData = false; @@ -600,10 +602,15 @@ class LayerTranslator const GDALVectorTranslateOptions *psOptions); private: - std::pair - GetDstClipGeom(const OGRSpatialReference *poGeomSRS); - std::pair - GetSrcClipGeom(const OGRSpatialReference *poGeomSRS); + struct ClipGeomDesc + { + const OGRGeometry *poGeom = nullptr; + const OGREnvelope *poEnv = nullptr; + bool bGeomIsRectangle = false; + }; + + ClipGeomDesc GetDstClipGeom(const OGRSpatialReference *poGeomSRS); + ClipGeomDesc GetSrcClipGeom(const OGRSpatialReference *poGeomSRS); }; static OGRLayer *GetLayerAndOverwriteIfNecessary(GDALDataset *poDstDS, @@ -6317,16 +6324,17 @@ bool LayerTranslator::Translate( if (poStolenGeometry->IsEmpty()) goto end_loop; - const auto [poClipGeom, poClipGeomEnvelope] = + const auto clipGeomDesc = GetSrcClipGeom(poStolenGeometry->getSpatialReference()); - if (poClipGeom && poClipGeomEnvelope) + if (clipGeomDesc.poGeom && clipGeomDesc.poEnv) { OGREnvelope oEnv; poStolenGeometry->getEnvelope(&oEnv); - if (!poClipGeomEnvelope->Contains(oEnv) && - !(poClipGeomEnvelope->Intersects(oEnv) && - poClipGeom->Intersects(poStolenGeometry.get()))) + if (!clipGeomDesc.poEnv->Contains(oEnv) && + !(clipGeomDesc.poEnv->Intersects(oEnv) && + clipGeomDesc.poGeom->Intersects( + poStolenGeometry.get()))) { goto end_loop; } @@ -6568,22 +6576,23 @@ bool LayerTranslator::Translate( if (poDstGeometry->IsEmpty()) goto end_loop; - const auto [poClipGeom, poClipGeomEnvelope] = + const auto clipGeomDesc = GetSrcClipGeom(poDstGeometry->getSpatialReference()); - if (!(poClipGeom && poClipGeomEnvelope)) + if (!(clipGeomDesc.poGeom && clipGeomDesc.poEnv)) goto end_loop; OGREnvelope oDstEnv; poDstGeometry->getEnvelope(&oDstEnv); - if (!poClipGeomEnvelope->Contains(oDstEnv)) + if (!(clipGeomDesc.bGeomIsRectangle && + clipGeomDesc.poEnv->Contains(oDstEnv))) { std::unique_ptr poClipped; - if (poClipGeomEnvelope->Intersects(oDstEnv)) + if (clipGeomDesc.poEnv->Intersects(oDstEnv)) { - poClipped.reset( - poClipGeom->Intersection(poDstGeometry.get())); + poClipped.reset(clipGeomDesc.poGeom->Intersection( + poDstGeometry.get())); } if (poClipped == nullptr || poClipped->IsEmpty()) { @@ -6750,9 +6759,9 @@ bool LayerTranslator::Translate( if (poDstGeometry->IsEmpty()) goto end_loop; - auto [poClipGeom, poClipGeomEnvelope] = GetDstClipGeom( + const auto clipGeomDesc = GetDstClipGeom( poDstGeometry->getSpatialReference()); - if (!poClipGeom || !poClipGeomEnvelope) + if (!clipGeomDesc.poGeom || !clipGeomDesc.poEnv) { goto end_loop; } @@ -6760,13 +6769,15 @@ bool LayerTranslator::Translate( OGREnvelope oDstEnv; poDstGeometry->getEnvelope(&oDstEnv); - if (!poClipGeomEnvelope->Contains(oDstEnv)) + if (!(clipGeomDesc.bGeomIsRectangle && + clipGeomDesc.poEnv->Contains(oDstEnv))) { std::unique_ptr poClipped; - if (poClipGeomEnvelope->Intersects(oDstEnv)) + if (clipGeomDesc.poEnv->Intersects(oDstEnv)) { - poClipped.reset(poClipGeom->Intersection( - poDstGeometry.get())); + poClipped.reset( + clipGeomDesc.poGeom->Intersection( + poDstGeometry.get())); } if (poClipped == nullptr || poClipped->IsEmpty()) @@ -6975,7 +6986,7 @@ bool LayerTranslator::Translate( * expressed. * @return the destination clip geometry and its envelope, or (nullptr, nullptr) */ -std::pair +LayerTranslator::ClipGeomDesc LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS) { if (m_poClipDstReprojectedToDstSRS_SRS != poGeomSRS) @@ -6988,7 +6999,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS) if (m_poClipDstReprojectedToDstSRS->transformTo(poGeomSRS) != OGRERR_NONE) { - return std::make_pair(nullptr, nullptr); + return ClipGeomDesc(); } m_poClipDstReprojectedToDstSRS_SRS = poGeomSRS; } @@ -7014,8 +7025,13 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS) if (poGeom && !m_oClipDstEnv.IsInit()) { poGeom->getEnvelope(&m_oClipDstEnv); + m_bClipDstIsRectangle = poGeom->IsRectangle(); } - return std::make_pair(poGeom, poGeom ? &m_oClipDstEnv : nullptr); + ClipGeomDesc ret; + ret.poGeom = poGeom; + ret.poEnv = poGeom ? &m_oClipDstEnv : nullptr; + ret.bGeomIsRectangle = m_bClipDstIsRectangle; + return ret; } /************************************************************************/ @@ -7028,7 +7044,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS) * expressed. * @return the source clip geometry and its envelope, or (nullptr, nullptr) */ -std::pair +LayerTranslator::ClipGeomDesc LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS) { if (m_poClipSrcReprojectedToSrcSRS_SRS != poGeomSRS) @@ -7041,7 +7057,7 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS) if (m_poClipSrcReprojectedToSrcSRS->transformTo(poGeomSRS) != OGRERR_NONE) { - return std::make_pair(nullptr, nullptr); + return ClipGeomDesc(); } m_poClipSrcReprojectedToSrcSRS_SRS = poGeomSRS; } @@ -7066,8 +7082,13 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS) if (poGeom && !m_oClipSrcEnv.IsInit()) { poGeom->getEnvelope(&m_oClipSrcEnv); + m_bClipSrcIsRectangle = poGeom->IsRectangle(); } - return std::make_pair(poGeom, poGeom ? &m_oClipSrcEnv : nullptr); + ClipGeomDesc ret; + ret.poGeom = poGeom; + ret.poEnv = poGeom ? &m_oClipSrcEnv : nullptr; + ret.bGeomIsRectangle = m_bClipDstIsRectangle; + return ret; } /************************************************************************/ diff --git a/autotest/utilities/test_ogr2ogr_lib.py b/autotest/utilities/test_ogr2ogr_lib.py index 9903cc794863..c1c97119a0ad 100755 --- a/autotest/utilities/test_ogr2ogr_lib.py +++ b/autotest/utilities/test_ogr2ogr_lib.py @@ -1116,6 +1116,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 @@ -1162,6 +1173,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) @@ -1390,6 +1414,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 @@ -1436,6 +1471,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)