From ada5d1f34413ce046d337fefb3f6b43a1742eda2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 26 Dec 2024 16:31:48 +0100 Subject: [PATCH] OGR VRT: accept SrcRegion value to be any geometry type as well as SetSpatialFilter() (#11537) Fixes #11518 Co-authored-by: Dan Baston --- autotest/ogr/data/vrt/invalid.vrt | 4 -- autotest/ogr/ogr_vrt.py | 32 ++++++++++++ doc/source/drivers/vector/vrt.rst | 4 +- ogr/ogrsf_frmts/vrt/ogr_vrt.h | 2 + ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp | 75 +++++++++++++---------------- 5 files changed, 69 insertions(+), 48 deletions(-) diff --git a/autotest/ogr/data/vrt/invalid.vrt b/autotest/ogr/data/vrt/invalid.vrt index 38863a8dafad..3485ef97fa48 100644 --- a/autotest/ogr/data/vrt/invalid.vrt +++ b/autotest/ogr/data/vrt/invalid.vrt @@ -70,10 +70,6 @@ ../poly.shp foo - - ../poly.shp - POINT (0 1) - ../poly.shp diff --git a/autotest/ogr/ogr_vrt.py b/autotest/ogr/ogr_vrt.py index 8c46606b026a..0a2c80630d73 100755 --- a/autotest/ogr/ogr_vrt.py +++ b/autotest/ogr/ogr_vrt.py @@ -3386,3 +3386,35 @@ def test_ogr_vrt_warped_arrow(tmp_vsimem): assert ds.GetLayer(0).GetExtent() == pytest.approx( (166021.443080541, 500000, 0.0, 331593.179548329) ) + + +############################################################################### +# Test SetSpatialFilter() on a point and SrcRegion not being a polygon + + +@pytest.mark.require_geos +def test_ogr_vrt_srcregion_and_point_filter(): + + src_ds = ogr.Open( + """ + + data/poly.shp + MULTIPOLYGON(((478315 4762880,478315 4765610,481645 4765610,481645 4762880,478315 4762880))) + + """ + ) + lyr = src_ds.GetLayer(0) + + lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt("POINT(479751 4764703)")) + lyr.ResetReading() + assert lyr.GetNextFeature() is not None + assert lyr.GetNextFeature() is None + assert lyr.GetFeatureCount() == 1 + + lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt("POINT(-479751 -4764703)")) + lyr.ResetReading() + assert lyr.GetNextFeature() is None + assert lyr.GetFeatureCount() == 0 + + lyr.SetSpatialFilter(None) + assert lyr.GetFeatureCount() == 10 diff --git a/doc/source/drivers/vector/vrt.rst b/doc/source/drivers/vector/vrt.rst index a961bd2371ca..2d4022a33da6 100644 --- a/doc/source/drivers/vector/vrt.rst +++ b/doc/source/drivers/vector/vrt.rst @@ -247,10 +247,10 @@ layer name, and may have the following subelements: - **SrcRegion** (optional) : This element is used to - define an initial spatial filter for the source features. This spatial + define an initial spatial filter for the source features, such that only features intersecting ``SrcRegion`` will be returned in the VRT layer. This spatial filter will be combined with any spatial filter explicitly set on the VRT layer with the SetSpatialFilter() method. The value of the element - must be a valid WKT string defining a polygon. An optional **clip** + must be a valid WKT string defining a geometry in the spatial reference system of the source layer. An optional **clip** attribute can be set to "TRUE" to clip the geometries to the source region, otherwise the source geometries are not modified. diff --git a/ogr/ogrsf_frmts/vrt/ogr_vrt.h b/ogr/ogrsf_frmts/vrt/ogr_vrt.h index 84fc27dca65f..fbeb44394a6a 100644 --- a/ogr/ogrsf_frmts/vrt/ogr_vrt.h +++ b/ogr/ogrsf_frmts/vrt/ogr_vrt.h @@ -125,6 +125,8 @@ class OGRVRTLayer final : public OGRLayer bool bError; + bool m_bEmptyResultSet = false; + bool ParseGeometryField(CPLXMLNode *psNode, CPLXMLNode *psNodeParent, OGRVRTGeomFieldProps *poProps); diff --git a/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp b/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp index 98dcfca75fec..8890f93af726 100644 --- a/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp +++ b/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp @@ -468,13 +468,10 @@ bool OGRVRTLayer::ParseGeometryField(CPLXMLNode *psNode, { OGRGeometryFactory::createFromWkt(pszSrcRegion, nullptr, &poProps->poSrcRegion); - if (poProps->poSrcRegion == nullptr || - wkbFlatten(poProps->poSrcRegion->getGeometryType()) != wkbPolygon) + if (poProps->poSrcRegion == nullptr) { CPLError(CE_Warning, CPLE_AppDefined, - "Ignoring SrcRegion. It must be a valid WKT polygon"); - delete poProps->poSrcRegion; - poProps->poSrcRegion = nullptr; + "Ignoring SrcRegion. It must be a valid WKT geometry"); } poProps->bSrcClip = @@ -1211,9 +1208,9 @@ bool OGRVRTLayer::ResetSourceReading() } else { - OGRGeometry *poIntersection = + auto poIntersection = std::unique_ptr( apoGeomFieldProps[i]->poSrcRegion->Intersection( - m_poFilterGeom); + m_poFilterGeom)); if (poIntersection && !poIntersection->IsEmpty()) { poIntersection->getEnvelope(&sEnvelope); @@ -1225,7 +1222,6 @@ bool OGRVRTLayer::ResetSourceReading() sEnvelope.MinY = 0; sEnvelope.MaxY = 0; } - delete poIntersection; } } else @@ -1318,63 +1314,54 @@ bool OGRVRTLayer::ResetSourceReading() CPLFree(pszFilter); + m_bEmptyResultSet = false; + // Clear spatial filter (to be safe) for non direct geometries // and reset reading. if (m_iGeomFieldFilter < static_cast(apoGeomFieldProps.size()) && apoGeomFieldProps[m_iGeomFieldFilter]->eGeometryStyle == VGS_Direct && apoGeomFieldProps[m_iGeomFieldFilter]->iGeomField >= 0) { - OGRGeometry *poSpatialGeom = nullptr; + OGRGeometry *poNewSpatialGeom = nullptr; OGRGeometry *poSrcRegion = apoGeomFieldProps[m_iGeomFieldFilter]->poSrcRegion; - bool bToDelete = false; + std::unique_ptr poIntersection; if (poSrcRegion == nullptr) { - poSpatialGeom = m_poFilterGeom; + poNewSpatialGeom = m_poFilterGeom; } else if (m_poFilterGeom == nullptr) { - poSpatialGeom = poSrcRegion; + poNewSpatialGeom = poSrcRegion; } else { - if (wkbFlatten(m_poFilterGeom->getGeometryType()) != wkbPolygon) - { - CPLError(CE_Failure, CPLE_AppDefined, - "Spatial filter should be polygon when a SrcRegion is " - "defined. Ignoring it"); - poSpatialGeom = poSrcRegion; - } - else + bool bDoIntersection = true; + if (m_bFilterIsEnvelope) { - bool bDoIntersection = true; - if (m_bFilterIsEnvelope) - { - OGREnvelope sEnvelope; - m_poFilterGeom->getEnvelope(&sEnvelope); - if (std::isinf(sEnvelope.MinX) && - std::isinf(sEnvelope.MinY) && - std::isinf(sEnvelope.MaxX) && - std::isinf(sEnvelope.MaxY) && sEnvelope.MinX < 0 && - sEnvelope.MinY < 0 && sEnvelope.MaxX > 0 && - sEnvelope.MaxY > 0) - { - poSpatialGeom = poSrcRegion; - bDoIntersection = false; - } - } - if (bDoIntersection) + OGREnvelope sEnvelope; + m_poFilterGeom->getEnvelope(&sEnvelope); + if (std::isinf(sEnvelope.MinX) && std::isinf(sEnvelope.MinY) && + std::isinf(sEnvelope.MaxX) && std::isinf(sEnvelope.MaxY) && + sEnvelope.MinX < 0 && sEnvelope.MinY < 0 && + sEnvelope.MaxX > 0 && sEnvelope.MaxY > 0) { - poSpatialGeom = m_poFilterGeom->Intersection(poSrcRegion); - bToDelete = true; + poNewSpatialGeom = poSrcRegion; + bDoIntersection = false; } } + if (bDoIntersection) + { + poIntersection.reset(m_poFilterGeom->Intersection(poSrcRegion)); + poNewSpatialGeom = poIntersection.get(); + if (!poIntersection) + m_bEmptyResultSet = true; + } } poSrcLayer->SetSpatialFilter( - apoGeomFieldProps[m_iGeomFieldFilter]->iGeomField, poSpatialGeom); - if (bToDelete) - delete poSpatialGeom; + apoGeomFieldProps[m_iGeomFieldFilter]->iGeomField, + poNewSpatialGeom); } else { @@ -1393,6 +1380,8 @@ bool OGRVRTLayer::ResetSourceReading() OGRFeature *OGRVRTLayer::GetNextFeature() { + if (m_bEmptyResultSet) + return nullptr; if (!bHasFullInitialized) FullInitialize(); if (!poSrcLayer || poDS->GetRecursionDetected()) @@ -2218,6 +2207,8 @@ OGRErr OGRVRTLayer::GetExtent(int iGeomField, OGREnvelope *psExtent, int bForce) GIntBig OGRVRTLayer::GetFeatureCount(int bForce) { + if (m_bEmptyResultSet) + return 0; if (nFeatureCount >= 0 && m_poFilterGeom == nullptr && m_poAttrQuery == nullptr) {