Skip to content

Commit 77eb999

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

File tree

3 files changed

+143
-53
lines changed

3 files changed

+143
-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: 52 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,10 @@
1414

1515
#include "cpl_conv.h"
1616
#include "ogr_pg.h"
17+
#include "ogr_p.h"
18+
19+
#include <algorithm>
20+
#include <limits>
1721

1822
#define PQexec this_is_an_error
1923

@@ -270,8 +274,8 @@ OGRFeature *OGRPGResultLayer::GetNextFeature()
270274
return nullptr;
271275

272276
if ((m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
273-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
274-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY ||
277+
(poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
278+
poGeomFieldDefn->nSRSId > 0) ||
275279
FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter))) &&
276280
(m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
277281
return poFeature;
@@ -297,35 +301,55 @@ OGRErr OGRPGResultLayer::ISetSpatialFilter(int iGeomField,
297301
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
298302
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
299303
{
300-
if (m_poFilterGeom != nullptr)
304+
poGeomFieldDefn->GetSpatialRef(); // make sure nSRSId is resolved
305+
if (m_poFilterGeom != nullptr && poDS->sPostGISVersion.nMajor >= 0)
301306
{
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)
307+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY &&
308+
poGeomFieldDefn->nSRSId > 0)
309+
{
310+
char *pszHexEWKB = OGRGeometryToHexEWKB(
311+
m_poFilterGeom, poGeomFieldDefn->nSRSId,
312+
poDS->sPostGISVersion.nMajor,
313+
poDS->sPostGISVersion.nMinor);
314+
osWHERE.Printf(
315+
"WHERE ST_Intersects(%s, '%s'::GEOMETRY) ",
316+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef())
317+
.c_str(),
318+
pszHexEWKB);
319+
CPLFree(pszHexEWKB);
320+
}
321+
else
308322
{
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;
323+
OGREnvelope sEnvelope;
324+
325+
m_poFilterGeom->getEnvelope(&sEnvelope);
326+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY)
327+
{
328+
if (sEnvelope.MinX < -180.0)
329+
sEnvelope.MinX = -180.0;
330+
if (sEnvelope.MinY < -90.0)
331+
sEnvelope.MinY = -90.0;
332+
if (sEnvelope.MaxX > 180.0)
333+
sEnvelope.MaxX = 180.0;
334+
if (sEnvelope.MaxY > 90.0)
335+
sEnvelope.MaxY = 90.0;
336+
}
337+
// Avoid +/- infinity
338+
sEnvelope.MinX = std::max(
339+
sEnvelope.MinX, -std::numeric_limits<double>::max());
340+
sEnvelope.MinY = std::max(
341+
sEnvelope.MinY, -std::numeric_limits<double>::max());
342+
sEnvelope.MaxX = std::min(
343+
sEnvelope.MaxX, std::numeric_limits<double>::max());
344+
sEnvelope.MaxY = std::min(
345+
sEnvelope.MaxY, std::numeric_limits<double>::max());
346+
osWHERE.Printf(
347+
"WHERE %s && ST_MakeEnvelope(%.17g,%.17g,%.17g,%.17g) ",
348+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef())
349+
.c_str(),
350+
sEnvelope.MinX, sEnvelope.MinY, sEnvelope.MaxX,
351+
sEnvelope.MaxY);
317352
}
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);
329353
}
330354
else
331355
{

‎ogr/ogrsf_frmts/pg/ogrpgtablelayer.cpp‎

Lines changed: 44 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,10 @@
1818
#include "cpl_error.h"
1919
#include "ogr_p.h"
2020

21+
#include <algorithm>
2122
#include <chrono>
2223
#include <condition_variable>
24+
#include <limits>
2325
#include <mutex>
2426
#include <thread>
2527

@@ -1057,30 +1059,48 @@ void OGRPGTableLayer::BuildWhere()
10571059
(poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
10581060
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY))
10591061
{
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);
1062+
poGeomFieldDefn->GetSpatialRef(); // make sure nSRSId is resolved
1063+
if (poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY)
1064+
{
1065+
char *pszHexEWKB = OGRGeometryToHexEWKB(
1066+
m_poFilterGeom, poGeomFieldDefn->nSRSId,
1067+
poDS->sPostGISVersion.nMajor, poDS->sPostGISVersion.nMinor);
1068+
osWHERE.Printf(
1069+
"WHERE ST_Intersects(%s, '%s'::GEOMETRY) ",
1070+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
1071+
pszHexEWKB);
1072+
CPLFree(pszHexEWKB);
1073+
}
1074+
else
1075+
{
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+
// Avoid +/- infinity
1091+
sEnvelope.MinX =
1092+
std::max(sEnvelope.MinX, -std::numeric_limits<double>::max());
1093+
sEnvelope.MinY =
1094+
std::max(sEnvelope.MinY, -std::numeric_limits<double>::max());
1095+
sEnvelope.MaxX =
1096+
std::min(sEnvelope.MaxX, std::numeric_limits<double>::max());
1097+
sEnvelope.MaxY =
1098+
std::min(sEnvelope.MaxY, std::numeric_limits<double>::max());
1099+
osWHERE.Printf(
1100+
"WHERE %s && ST_MakeEnvelope(%.17g,%.17g,%.17g,%.17g) ",
1101+
OGRPGEscapeColumnName(poGeomFieldDefn->GetNameRef()).c_str(),
1102+
sEnvelope.MinX, sEnvelope.MinY, sEnvelope.MaxX, sEnvelope.MaxY);
1103+
}
10841104
}
10851105

10861106
if (!osQuery.empty())
@@ -1174,7 +1194,6 @@ OGRFeature *OGRPGTableLayer::GetNextFeature()
11741194
* request */
11751195
if (m_poFilterGeom == nullptr || poGeomFieldDefn == nullptr ||
11761196
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOMETRY ||
1177-
poGeomFieldDefn->ePostgisType == GEOM_TYPE_GEOGRAPHY ||
11781197
FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter)))
11791198
{
11801199
if (iFIDAsRegularColumnIndex >= 0)

0 commit comments

Comments
 (0)