Merge pull request #8626 from rouault/fix_8625

GPKG: make GetFeatureCount() do full geometry intersection and not just bounding box (fixes #8625)

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

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);

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<OGRGeoPackageTableLayer *>(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<const GByte *>(sqlite3_value_blob(argv[0]));
+    auto poGeom = std::unique_ptr<OGRGeometry>(
+        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()                           */
 /************************************************************************/

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)
     {