From faa41c20f1e3a215443e8b761dcd9cfa3f5924f8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Oct 2023 17:46:03 +0200 Subject: [PATCH] GPKG: make GetFeatureCount() do full geometry intersection and not just bounding box (fixes #8625) --- autotest/ogr/ogr_gpkg.py | 27 +++++++++ ogr/ogrsf_frmts/gpkg/ogr_geopackage.h | 8 +++ .../gpkg/ogrgeopackagedatasource.cpp | 56 +++++++++++++++++++ .../gpkg/ogrgeopackagetablelayer.cpp | 39 ++++++++++++- 4 files changed, 127 insertions(+), 3 deletions(-) diff --git a/autotest/ogr/ogr_gpkg.py b/autotest/ogr/ogr_gpkg.py index 181d4c000320..22ab0732f05b 100755 --- a/autotest/ogr/ogr_gpkg.py +++ b/autotest/ogr/ogr_gpkg.py @@ -9488,3 +9488,30 @@ def test_ogr_gpkg_sql_first_geom_null(): sql_lyr.GetNextFeature() f = sql_lyr.GetNextFeature() assert f.GetGeometryRef().GetGeometryType() == ogr.wkbPolygon + + +############################################################################### + + +@pytest.mark.require_geos +def test_ogr_gpkg_sql_exact_spatial_filter_for_feature_count(tmp_vsimem): + + filename = str( + tmp_vsimem / "test_ogr_gpkg_sql_exact_spatial_filter_for_feature_count.gpkg" + ) + ds = gdal.GetDriverByName("GPKG").Create(filename, 0, 0, 0, gdal.GDT_Unknown) + lyr = ds.CreateLayer( + "test", + ) + f = ogr.Feature(lyr.GetLayerDefn()) + f.SetGeometry(ogr.CreateGeometryFromWkt("LINESTRING (0 0,1 1)")) + lyr.CreateFeature(f) + f = ogr.Feature(lyr.GetLayerDefn()) + f.SetGeometry(ogr.CreateGeometryFromWkt("LINESTRING (0.1 0.25,0.15 0.25)")) + lyr.CreateFeature(f) + ds = None + + ds = ogr.Open(filename) + lyr = ds.GetLayer(0) + lyr.SetSpatialFilterRect(0.1, 0.2, 0.15, 0.3) + assert lyr.GetFeatureCount() == 1 diff --git a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h index 7a68ef6b1c40..f9cecfeb656d 100644 --- a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h +++ b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h @@ -117,6 +117,9 @@ struct OGRGPKGTableLayerFillArrowArray const OGRLayer *poLayerForFilterGeom = nullptr; }; +void OGR_GPKG_Intersects_Spatial_Filter(sqlite3_context *pContext, int argc, + sqlite3_value **argv); + /************************************************************************/ /* GDALGeoPackageDataset */ /************************************************************************/ @@ -764,6 +767,11 @@ class OGRGeoPackageTableLayer final : public OGRGeoPackageLayer void GetNextArrowArrayAsynchronousWorker(); void CancelAsyncNextArrowArray(); + protected: + friend void OGR_GPKG_Intersects_Spatial_Filter(sqlite3_context *pContext, + int /*argc*/, + sqlite3_value **argv); + public: OGRGeoPackageTableLayer(GDALGeoPackageDataset *poDS, const char *pszTableName); diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp index 22f10423a13f..7cbf4fe3af7a 100644 --- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp +++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp @@ -7706,6 +7706,7 @@ static bool OGRGeoPackageGetHeader(sqlite3_context *pContext, int /*argc*/, &(psHeader->MaxY)) == OGRERR_NONE) { psHeader->bEmpty = bEmpty; + psHeader->bExtentHasXY = !bEmpty; if (!(bEmpty && bNeedExtent)) return true; } @@ -7739,6 +7740,61 @@ static bool OGRGeoPackageGetHeader(sqlite3_context *pContext, int /*argc*/, return true; } +/************************************************************************/ +/* OGR_GPKG_Intersects_Spatial_Filter() */ +/************************************************************************/ + +void OGR_GPKG_Intersects_Spatial_Filter(sqlite3_context *pContext, int argc, + sqlite3_value **argv) +{ + if (sqlite3_value_type(argv[0]) != SQLITE_BLOB) + { + sqlite3_result_int(pContext, 0); + return; + } + + auto poLayer = + static_cast(sqlite3_user_data(pContext)); + + GPkgHeader sHeader; + if (poLayer->m_bFilterIsEnvelope && + OGRGeoPackageGetHeader(pContext, argc, argv, &sHeader, false) && + sHeader.bExtentHasXY) + { + OGREnvelope sEnvelope; + sEnvelope.MinX = sHeader.MinX; + sEnvelope.MinY = sHeader.MinY; + sEnvelope.MaxX = sHeader.MaxX; + sEnvelope.MaxY = sHeader.MaxY; + if (poLayer->m_sFilterEnvelope.Contains(sEnvelope)) + { + sqlite3_result_int(pContext, 1); + return; + } + } + + const int nBLOBLen = sqlite3_value_bytes(argv[0]); + const GByte *pabyBLOB = + reinterpret_cast(sqlite3_value_blob(argv[0])); + auto poGeom = std::unique_ptr( + GPkgGeometryToOGR(pabyBLOB, nBLOBLen, nullptr)); + if (poGeom == nullptr) + { + // Try also spatialite geometry blobs + OGRGeometry *poGeomSpatialite = nullptr; + if (OGRSQLiteImportSpatiaLiteGeometry(pabyBLOB, nBLOBLen, + &poGeomSpatialite) != OGRERR_NONE) + { + CPLError(CE_Failure, CPLE_AppDefined, "Invalid geometry"); + sqlite3_result_int(pContext, 0); + return; + } + poGeom.reset(poGeomSpatialite); + } + + sqlite3_result_int(pContext, poLayer->FilterGeometry(poGeom.get())); +} + /************************************************************************/ /* OGRGeoPackageSTMinX() */ /************************************************************************/ diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp index 1d3d8654e7ed..df105e5a0728 100644 --- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp +++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp @@ -3858,9 +3858,6 @@ GIntBig OGRGeoPackageTableLayer::GetFeatureCount(int /*bForce*/) } #endif - if (m_poFilterGeom != nullptr && !m_bFilterIsEnvelope) - return OGRGeoPackageLayer::GetFeatureCount(); - if (m_bDeferredCreation && RunDeferredCreationIfNecessary() != OGRERR_NONE) return 0; @@ -3869,6 +3866,7 @@ GIntBig OGRGeoPackageTableLayer::GetFeatureCount(int /*bForce*/) /* Ignore bForce, because we always do a full count on the database */ OGRErr err; CPLString soSQL; + bool bUnregisterSQLFunction = false; if (m_bIsTable && m_poFilterGeom != nullptr && m_pszAttrQueryString == nullptr && HasSpatialIndex()) { @@ -3885,6 +3883,34 @@ GIntBig OGRGeoPackageTableLayer::GetFeatureCount(int /*bForce*/) SQLEscapeName(m_osRTreeName).c_str(), sEnvelope.MinX - 1e-11, sEnvelope.MaxX + 1e-11, sEnvelope.MinY - 1e-11, sEnvelope.MaxY + 1e-11); + + if (OGRGeometryFactory::haveGEOS() && + !(m_bFilterIsEnvelope && + wkbFlatten( + m_poFeatureDefn->GetGeomFieldDefn(m_iGeomFieldFilter) + ->GetType()) == wkbPoint)) + { + bUnregisterSQLFunction = true; + sqlite3_create_function( + m_poDS->hDB, "OGR_GPKG_Intersects_Spatial_Filter", 1, + SQLITE_UTF8, this, OGR_GPKG_Intersects_Spatial_Filter, + nullptr, nullptr); + const char *pszC = + m_poFeatureDefn->GetGeomFieldDefn(m_iGeomFieldFilter) + ->GetNameRef(); + soSQL.Printf("SELECT COUNT(*) FROM \"%s\" m " + "JOIN \"%s\" r " + "ON m.\"%s\" = r.id WHERE " + "r.maxx >= %.12f AND r.minx <= %.12f AND " + "r.maxy >= %.12f AND r.miny <= %.12f AND " + "OGR_GPKG_Intersects_Spatial_Filter(m.\"%s\")", + SQLEscapeName(m_pszTableName).c_str(), + SQLEscapeName(m_osRTreeName).c_str(), + SQLEscapeName(m_osFIDForRTree).c_str(), + sEnvelope.MinX - 1e-11, sEnvelope.MaxX + 1e-11, + sEnvelope.MinY - 1e-11, sEnvelope.MaxY + 1e-11, + SQLEscapeName(pszC).c_str()); + } } } @@ -3903,6 +3929,13 @@ GIntBig OGRGeoPackageTableLayer::GetFeatureCount(int /*bForce*/) GIntBig iFeatureCount = SQLGetInteger64(m_poDS->GetDB(), soSQL.c_str(), &err); + if (bUnregisterSQLFunction) + { + sqlite3_create_function(m_poDS->hDB, + "OGR_GPKG_Intersects_Spatial_Filter", 1, + SQLITE_UTF8, this, nullptr, nullptr, nullptr); + } + /* Generic implementation uses -1 for error condition, so we will too */ if (err == OGRERR_NONE) {