VeniVidiVici
/
GDAL2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

OGR VRT: accept SrcRegion value to be any geometry type as well as SetSpatialFilter() (#11537) 

Fixes #11518

Co-authored-by: Dan Baston <dbaston@gmail.com>
This commit is contained in:
Even Rouault 2024-12-26 16:31:48 +01:00 committed by GitHub
parent af66409096
commit ada5d1f344
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
5 changed files with 69 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
autotest/ogr/data/vrt/invalid.vrt
View File

@ -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"/>

32
autotest/ogr/ogr_vrt.py
View File

@ -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

4
doc/source/drivers/vector/vrt.rst
View File

@ -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.

2
ogr/ogrsf_frmts/vrt/ogr_vrt.h
View File

@ -125,6 +125,8 @@ class OGRVRTLayer final : public OGRLayer
bool bError;
bool m_bEmptyResultSet = false;
bool ParseGeometryField(CPLXMLNode *psNode, CPLXMLNode *psNodeParent,
OGRVRTGeomFieldProps *poProps);

75
ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp
View File

@ -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)
bool bDoIntersection = true;
if (m_bFilterIsEnvelope)
{
CPLError(CE_Failure, CPLE_AppDefined,
"Spatial filter should be polygon when a SrcRegion is "
"defined. Ignoring it");
poSpatialGeom = poSrcRegion;
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)
{
poNewSpatialGeom = poSrcRegion;
bDoIntersection = false;
}
}
else
if (bDoIntersection)
{
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)
{
poSpatialGeom = m_poFilterGeom->Intersection(poSrcRegion);
bToDelete = true;
}
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)
{