Skip to content

Commit 268799e

Browse files
rouaultrouault
committed
PostGIS: use full geometry intersection for SetSpatialFilter(), and not only bounding box
Fixes OSGeo#13807
1 parent 81c7c68 commit 268799e

File tree

3 files changed

+133
-53
lines changed

3 files changed

+133
-53
lines changed

‎autotest/ogr/ogr_pg.py‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6376,3 +6376,50 @@ def test_ogr_pg_field_truncation(pg_ds):
63766376
lyr.ResetReading()
63776377
f = lyr.GetNextFeature()
63786378
assert f["field"] == b"abcd\xc3\xa9".decode("UTF-8")
6379+
6380+
6381+
###############################################################################
6382+
# Test real geometry intersection in spatial filter
6383+
6384+
6385+
@gdaltest.enable_exceptions()
6386+
@pytest.mark.parametrize("GEOM_TYPE", ["geometry", "geography"])
6387+
def test_ogr_pg_geometry_intersection_spatial_filter(pg_ds, use_postgis, GEOM_TYPE):
6388+
6389+
if use_postgis:
6390+
srs = osr.SpatialReference(epsg=4326)
6391+
lyr = pg_ds.CreateLayer("test", srs, options=["GEOM_TYPE=" + GEOM_TYPE])
6392+
else:
6393+
lyr = pg_ds.CreateLayer("test")
6394+
f = ogr.Feature(lyr.GetLayerDefn())
6395+
f.SetGeometry(ogr.CreateGeometryFromWkt("POINT (0.5 0.5)"))
6396+
lyr.CreateFeature(f)
6397+
f = ogr.Feature(lyr.GetLayerDefn())
6398+
f.SetGeometry(ogr.CreateGeometryFromWkt("POINT (0.9 0.5)"))
6399+
lyr.CreateFeature(f)
6400+
6401+
lyr.ResetReading()
6402+
6403+
# Pacman style polygon: almost a square except it doesn't include 0.9, 0.5
6404+
lyr.SetSpatialFilter(
6405+
ogr.CreateGeometryFromWkt(
6406+
"POLYGON ((0 0,0 1,1 1,1 0.6,0.8 0.6,0.8 0.4,1 0.4,1 0,0 0))"
6407+
)
6408+
)
6409+
f = lyr.GetNextFeature()
6410+
assert f.GetGeometryRef().ExportToWkt() == "POINT (0.5 0.5)"
6411+
f = lyr.GetNextFeature()
6412+
assert f is None
6413+
6414+
lyr.SetSpatialFilter(None)
6415+
6416+
with pg_ds.ExecuteSQL("SELECT * FROM test") as sql_lyr:
6417+
sql_lyr.SetSpatialFilter(
6418+
ogr.CreateGeometryFromWkt(
6419+
"POLYGON ((0 0,0 1,1 1,1 0.6,0.8 0.6,0.8 0.4,1 0.4,1 0,0 0))"
6420+
)
6421+
)
6422+
f = sql_lyr.GetNextFeature()
6423+
assert f.GetGeometryRef().ExportToWkt() == "POINT (0.5 0.5)"
6424+
f = sql_lyr.GetNextFeature()
6425+
assert f is None

‎ogr/ogrsf_frmts/pg/ogrpgresultlayer.cpp‎

Lines changed: 46 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -270,8 +270,8 @@ OGRFeature *OGRPGResultLayer::GetNextFeature()
270270
return nullptr;
271271

272272
if ((m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
273-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
274-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY ||
273+
(poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
274+
poGeomFieldDefn->nSRSId > 0) ||
275275
FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter))) &&
276276
(m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
277277
return poFeature;
@@ -297,35 +297,53 @@ OGRErr OGRPGResultLayer::ISetSpatialFilter(int iGeomField,
297297
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
298298
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
299299
{
300-
if (m_poFilterGeom != nullptr)
300+
poGeomFieldDefn->GetSpatialRef(); // make sure nSRSId is resolved
301+
if (m_poFilterGeom != nullptr && poDS->sPostGISVersion.nMajor >= 0)
301302
{
302-
char szBox3D_1[128];
303-
char szBox3D_2[128];
304-
OGREnvelope sEnvelope;
305-
306-
m_poFilterGeom->getEnvelope(&sEnvelope);
307-
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
303+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
304+
poGeomFieldDefn->nSRSId > 0)
305+
{
306+
OGRWktOptions sOptions;
307+
sOptions.variant = wkbVariantIso;
308+
const std::string osWKT =
309+
m_poFilterGeom->exportToWkt(sOptions);
310+
osWHERE.Printf(
311+
"WHERE ST_Intersects(%s, ST_SetSRID("
312+
"ST_GeomFromText('%s'), %d)) ",
313+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef())
314+
.c_str(),
315+
osWKT.c_str(), poGeomFieldDefn->nSRSId);
316+
}
317+
else
308318
{
309-
if (sEnvelope.MinX < -180.0)
310-
sEnvelope.MinX = -180.0;
311-
if (sEnvelope.MinY < -90.0)
312-
sEnvelope.MinY = -90.0;
313-
if (sEnvelope.MaxX > 180.0)
314-
sEnvelope.MaxX = 180.0;
315-
if (sEnvelope.MaxY > 90.0)
316-
sEnvelope.MaxY = 90.0;
319+
char szBox3D_1[128];
320+
char szBox3D_2[128];
321+
OGREnvelope sEnvelope;
322+
323+
m_poFilterGeom->getEnvelope(&sEnvelope);
324+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
325+
{
326+
if (sEnvelope.MinX < -180.0)
327+
sEnvelope.MinX = -180.0;
328+
if (sEnvelope.MinY < -90.0)
329+
sEnvelope.MinY = -90.0;
330+
if (sEnvelope.MaxX > 180.0)
331+
sEnvelope.MaxX = 180.0;
332+
if (sEnvelope.MaxY > 90.0)
333+
sEnvelope.MaxY = 90.0;
334+
}
335+
CPLsnprintf(szBox3D_1, sizeof(szBox3D_1), "%.17g %.17g",
336+
sEnvelope.MinX, sEnvelope.MinY);
337+
CPLsnprintf(szBox3D_2, sizeof(szBox3D_2), "%.17g %.17g",
338+
sEnvelope.MaxX, sEnvelope.MaxY);
339+
osWHERE.Printf(
340+
"WHERE %s && %s('BOX3D(%s, %s)'::box3d,%d) ",
341+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef())
342+
.c_str(),
343+
(poDS->sPostGISVersion.nMajor >= 2) ? "ST_SetSRID"
344+
: "SetSRID",
345+
szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
317346
}
318-
CPLsnprintf(szBox3D_1, sizeof(szBox3D_1), "%.17g %.17g",
319-
sEnvelope.MinX, sEnvelope.MinY);
320-
CPLsnprintf(szBox3D_2, sizeof(szBox3D_2), "%.17g %.17g",
321-
sEnvelope.MaxX, sEnvelope.MaxY);
322-
osWHERE.Printf(
323-
"WHERE %s && %s('BOX3D(%s, %s)'::box3d,%d) ",
324-
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef())
325-
.c_str(),
326-
(poDS->sPostGISVersion.nMajor >= 2) ? "ST_SetSRID"
327-
: "SetSRID",
328-
szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
329347
}
330348
else
331349
{

‎ogr/ogrsf_frmts/pg/ogrpgtablelayer.cpp‎

Lines changed: 40 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1057,30 +1057,46 @@ void OGRPGTableLayer::BuildWhere()
10571057
(poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
10581058
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY))
10591059
{
1060-
char szBox3D_1[128];
1061-
char szBox3D_2[128];
1062-
OGREnvelope sEnvelope;
1063-
1064-
m_poFilterGeom->getEnvelope(&sEnvelope);
1065-
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
1066-
{
1067-
if (sEnvelope.MinX < -180.0)
1068-
sEnvelope.MinX = -180.0;
1069-
if (sEnvelope.MinY < -90.0)
1070-
sEnvelope.MinY = -90.0;
1071-
if (sEnvelope.MaxX > 180.0)
1072-
sEnvelope.MaxX = 180.0;
1073-
if (sEnvelope.MaxY > 90.0)
1074-
sEnvelope.MaxY = 90.0;
1075-
}
1076-
CPLsnprintf(szBox3D_1, sizeof(szBox3D_1), "%.17g %.17g", sEnvelope.MinX,
1077-
sEnvelope.MinY);
1078-
CPLsnprintf(szBox3D_2, sizeof(szBox3D_2), "%.17g %.17g", sEnvelope.MaxX,
1079-
sEnvelope.MaxY);
1080-
osWHERE.Printf(
1081-
"WHERE %s && ST_SetSRID('BOX3D(%s, %s)'::box3d,%d) ",
1082-
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
1083-
szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
1060+
poGeomFieldDefn->GetSpatialRef(); // make sure nSRSId is resolved
1061+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY)
1062+
{
1063+
OGRWktOptions sOptions;
1064+
sOptions.variant = wkbVariantIso;
1065+
const std::string osWKT = m_poFilterGeom->exportToWkt(sOptions);
1066+
osWHERE.Printf(
1067+
"WHERE ST_Intersects(%s, ST_SetSRID(ST_GeomFromText('%s'), "
1068+
"%d)) ",
1069+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
1070+
osWKT.c_str(), poGeomFieldDefn->nSRSId);
1071+
}
1072+
else
1073+
{
1074+
char szBox3D_1[128];
1075+
char szBox3D_2[128];
1076+
OGREnvelope sEnvelope;
1077+
1078+
m_poFilterGeom->getEnvelope(&sEnvelope);
1079+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
1080+
{
1081+
if (sEnvelope.MinX < -180.0)
1082+
sEnvelope.MinX = -180.0;
1083+
if (sEnvelope.MinY < -90.0)
1084+
sEnvelope.MinY = -90.0;
1085+
if (sEnvelope.MaxX > 180.0)
1086+
sEnvelope.MaxX = 180.0;
1087+
if (sEnvelope.MaxY > 90.0)
1088+
sEnvelope.MaxY = 90.0;
1089+
}
1090+
CPLsnprintf(szBox3D_1, sizeof(szBox3D_1), "%.17g %.17g",
1091+
sEnvelope.MinX, sEnvelope.MinY);
1092+
CPLsnprintf(szBox3D_2, sizeof(szBox3D_2), "%.17g %.17g",
1093+
sEnvelope.MaxX, sEnvelope.MaxY);
1094+
1095+
osWHERE.Printf(
1096+
"WHERE %s && ST_SetSRID('BOX3D(%s, %s)'::box3d,%d) ",
1097+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
1098+
szBox3D_1, szBox3D_2, poGeomFieldDefn->nSRSId);
1099+
}
10841100
}
10851101

10861102
if (!osQuery.empty())
@@ -1174,7 +1190,6 @@ OGRFeature *OGRPGTableLayer::GetNextFeature()
11741190
* request */
11751191
if (m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
11761192
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
1177-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY ||
11781193
FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter)))
11791194
{
11801195
if (iFIDAsRegularColumnIndex >= 0)

0 commit comments

Comments
 (0)