OGR VRT: accept SrcRegion value to be any geometry type as well as Se…

…tSpatialFilter() (#11537) Fixes #11518 Co-authored-by: Dan Baston <dbaston@gmail.com>
OSGeo · Dec 26, 2024 · ada5d1f · ada5d1f
1 parent af66409
commit ada5d1f
Show file tree

Hide file tree

Showing 5 changed files with 69 additions and 48 deletions.
diff --git a/autotest/ogr/data/vrt/invalid.vrt b/autotest/ogr/data/vrt/invalid.vrt
@@ -70,10 +70,6 @@
         <SrcDataSource relativeToVRT="1">../poly.shp</SrcDataSource>
         <SrcRegion>foo</SrcRegion>
     </OGRVRTLayer>
-    <OGRVRTLayer name="poly">
-        <SrcDataSource relativeToVRT="1">../poly.shp</SrcDataSource>
-        <SrcRegion>POINT (0 1)</SrcRegion>
-    </OGRVRTLayer>
     <OGRVRTLayer name="poly">
         <SrcDataSource relativeToVRT="1">../poly.shp</SrcDataSource>
         <GeometryField name="foo" field="bar"/>

diff --git 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(
+        """<OGRVRTDataSource>
+            <OGRVRTLayer name="poly">
+                <SrcDataSource>data/poly.shp</SrcDataSource>
+                <SrcRegion>MULTIPOLYGON(((478315 4762880,478315 4765610,481645 4765610,481645 4762880,478315 4762880)))</SrcRegion>
+            </OGRVRTLayer>
+    </OGRVRTDataSource>"""
+    )
+    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
@@ -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
@@ -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
@@ -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<OGRGeometry>(
                             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<int>(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<OGRGeometry> 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)
     {