OGR/GPKG: add a gdal_get_pixel_value() SQL function.

The ``gdal_get_pixel_value()`` function (added in GDAL 3.7) can be used to extract the value
of a pixel in a GDAL dataset. It requires the configuration option OGR_SQLITE_ALLOW_EXTERNAL_ACCESS
to be set to YES (for security reasons).

It takes 5 arguments:

* a string with the dataset name
* a band number (numbering starting at 1)
* a string being "georef" to indicate that subsequent values will be georeferenced
  coordinates, or "pixel" to indicate that subsequent values will be in column, line
  pixel space
* georeferenced X value or column number
* georeferenced Y value or line number

.. code-block::

    SELECT gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'georef', 440720, 3751320)
    SELECT gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', 0, 0)

autotest/ogr/ogr_gpkg.py

@@ -8055,3 +8055,135 @@ def test_ogr_gpkg_read_generated_column():
     ds = None
 
     gdal.Unlink(filename)
+
+
+###############################################################################
+# Test gdal_get_pixel_value() function
+
+
+def test_ogr_gpkg_sql_gdal_get_pixel_value():
+
+    filename = "/vsimem/test_ogr_gpkg_sql_gdal_get_pixel_value.gpkg"
+    ds = ogr.GetDriverByName("GPKG").CreateDataSource(filename)
+
+    with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'georef', 440720, 3751320)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 107
+
+    with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', 0, 0)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 107
+
+    with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', 0, 1)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 115
+
+    with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/float64.tif', 1, 'pixel', 0, 1)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 115.0
+
+    # Invalid column
+    with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', -1, 0)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] is None
+
+    # Missing OGR_SQLITE_ALLOW_EXTERNAL_ACCESS
+    with gdaltest.error_handler():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'georef', 440720, 3751320)"
+        )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # NULL as 1st arg
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value(NULL, 1, 'pixel', 0, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # NULL as 2nd arg
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', NULL, 'pixel', 0, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # NULL as 3rd arg
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, NULL, 0, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # NULL as 4th arg
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', NULL, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # NULL as 5th arg
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', 0, NULL)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # Invalid band number
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', 0, 'pixel', 0, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # Invalid value for 3rd argument
+    with gdaltest.error_handler():
+        with gdaltest.config_option("OGR_SQLITE_ALLOW_EXTERNAL_ACCESS", "YES"):
+            sql_lyr = ds.ExecuteSQL(
+                "select gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'invalid', 0, 0)"
+            )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    gdal.Unlink(filename)

doc/source/drivers/vector/gpkg.rst

@@ -124,6 +124,9 @@ Spatialite, are also available :
    in gpkg_spatial_ref_sys, starting with GDAL 3.2, it will be interpreted as
    a EPSG code.
 
+The raster SQL functions mentioned at :ref:`sql_sqlite_dialect_raster_functions`
+are also available.
+
 Link with Spatialite
 ~~~~~~~~~~~~~~~~~~~~
 

doc/source/user/sql_sqlite_dialect.rst

@@ -257,6 +257,31 @@ a value associate to a key from a HSTORE string, formatted like "key=>value,othe
 
     SELECT hstore_get_value('a => b, "key with space"=> "value with space"', 'key with space') --> 'value with space'
 
+.. _sql_sqlite_dialect_raster_functions:
+
+Raster related functions
+++++++++++++++++++++++++
+
+The ``gdal_get_pixel_value()`` function (added in GDAL 3.7) can be used to extract the value
+of a pixel in a GDAL dataset. It requires the configuration option OGR_SQLITE_ALLOW_EXTERNAL_ACCESS
+to be set to YES (for security reasons).
+
+It takes 5 arguments:
+
+* a string with the dataset name
+* a band number (numbering starting at 1)
+* a string being "georef" to indicate that subsequent values will be georeferenced
+  coordinates, or "pixel" to indicate that subsequent values will be in column, line
+  pixel space
+* georeferenced X value or column number
+* georeferenced Y value or line number
+
+.. code-block::
+
+    SELECT gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'georef', 440720, 3751320)
+    SELECT gdal_get_pixel_value('../gcore/data/byte.tif', 1, 'pixel', 0, 0)
+
+
 OGR geocoding functions
 +++++++++++++++++++++++
 

ogr/ogrsf_frmts/gpkg/ogr_geopackage.h

@@ -122,6 +122,7 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource, public GDALG
                                        int argc,
                                        sqlite3_value** argv);
 
+    void               *m_pSQLFunctionData = nullptr;
     GUInt32             m_nApplicationId = GPKG_APPLICATION_ID;
     GUInt32             m_nUserVersion = GPKG_1_2_VERSION;
     OGRGeoPackageTableLayer** m_papoLayers = nullptr;
@@ -272,6 +273,8 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource, public GDALG
         bool                m_bIsGeometryTypeAggregateInterrupted = false;
         std::string         m_osGeometryTypeAggregateResult{};
 
+        virtual void        CloseDB() override;
+
         CPL_DISALLOW_COPY_ASSIGN(GDALGeoPackageDataset)
 
     public:

ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp

@@ -45,6 +45,8 @@
 #include <limits>
 #include <sstream>
 
+#define COMPILATION_ALLOWED
+#include "ogrsqlitesqlfunctionscommon.cpp"
 
 // Keep in sync prototype of those 2 functions between gdalopeninfo.cpp,
 // ogrsqlitedatasource.cpp and ogrgeopackagedatasource.cpp
@@ -185,7 +187,6 @@ static std::unique_ptr<TilingSchemeDefinition> GetTilingScheme(const char* pszNa
     return poTilingScheme;
 }
 
-
 static const char* pszCREATE_GPKG_GEOMETRY_COLUMNS =
     "CREATE TABLE gpkg_geometry_columns ("
     "table_name TEXT NOT NULL,"
@@ -259,6 +260,13 @@ OGRErr GDALGeoPackageDataset::SetApplicationAndUserVersionId()
 #endif
 }
 
+void GDALGeoPackageDataset::CloseDB()
+{
+    OGRSQLiteUnregisterSQLFunctions(m_pSQLFunctionData);
+    m_pSQLFunctionData = nullptr;
+    OGRSQLiteBaseDataSource::CloseDB();
+}
+
 bool GDALGeoPackageDataset::ReOpenDB()
 {
     CPLAssert( hDB != nullptr );
@@ -999,6 +1007,8 @@ GDALGeoPackageDataset::~GDALGeoPackageDataset()
         if( poSRS )
             poSRS->Release();
     }
+
+    CloseDB();
 }
 
 /************************************************************************/
@@ -8203,6 +8213,8 @@ void GDALGeoPackageDataset::InstallSQLFunctions()
                                 SQLITE_UTF8 | SQLITE_DETERMINISTIC, nullptr,
                                 GPKG_GDAL_HasColorTable, nullptr, nullptr);
     }
+
+    m_pSQLFunctionData = OGRSQLiteRegisterSQLFunctionsCommon(hDB);
 }
 
 /************************************************************************/

ogr/ogrsf_frmts/sqlite/ogrsqlitebase.h

@@ -141,7 +141,7 @@ class OGRSQLiteBaseDataSource CPL_NON_FINAL: public GDALPamDataset
     bool                SetCacheSize();
     void                LoadExtensions();
 
-    void                CloseDB();
+    virtual void        CloseDB();
 
     std::map<CPLString, OGREnvelope> oMapSQLEnvelope{};
 

ogr/ogrsf_frmts/sqlite/ogrsqliteregexp.h

@@ -30,7 +30,7 @@
 #ifndef OGR_SQLITE_REGEXP_INCLUDED
 #define OGR_SQLITE_REGEXP_INCLUDED
 
-#include "ogr_sqlite.h"
+#include "sqlite3.h"
 
 static void* OGRSQLiteRegisterRegExpFunction(sqlite3* hDB);
 static void OGRSQLiteFreeRegExpCache(void* hRegExpCache);

ogr/ogrsf_frmts/sqlite/ogrsqlitesqlfunctions.cpp

@@ -35,7 +35,6 @@
 #include "ogrsqlitesqlfunctions.h"
 #include "ogr_geocoding.h"
 
-#include "ogrsqliteregexp.cpp" /* yes the .cpp file, to make it work on Windows with load_extension('gdalXX.dll') */
 #include "ogr_swq.h"
 
 #include <limits>
@@ -49,90 +48,11 @@
 #define MINIMAL_SPATIAL_FUNCTIONS
 #endif
 
-class OGRSQLiteExtensionData
-{
-#ifdef DEBUG
-    void* pDummy; /* to track memory leaks */
-#endif
-    std::map< std::pair<int,int>, OGRCoordinateTransformation*> oCachedTransformsMap;
-
-    void* hRegExpCache;
-
-    OGRGeocodingSessionH hGeocodingSession;
-
-  public:
-    explicit                     OGRSQLiteExtensionData(sqlite3* hDB);
-                                ~OGRSQLiteExtensionData();
-
-    OGRCoordinateTransformation* GetTransform(int nSrcSRSId, int nDstSRSId);
-
-    OGRGeocodingSessionH         GetGeocodingSession() { return hGeocodingSession; }
-    void                         SetGeocodingSession(OGRGeocodingSessionH hGeocodingSessionIn) { hGeocodingSession = hGeocodingSessionIn; }
-
-    void                         SetRegExpCache(void* hRegExpCacheIn) { hRegExpCache = hRegExpCacheIn; }
-};
-
-/************************************************************************/
-/*                     OGRSQLiteExtensionData()                         */
-/************************************************************************/
-
-OGRSQLiteExtensionData::OGRSQLiteExtensionData(CPL_UNUSED sqlite3* hDB) :
-#ifdef DEBUG
-    pDummy(CPLMalloc(1)),
-#endif
-    hRegExpCache(nullptr),
-    hGeocodingSession(nullptr)
-{}
-
-/************************************************************************/
-/*                       ~OGRSQLiteExtensionData()                      */
-/************************************************************************/
-
-OGRSQLiteExtensionData::~OGRSQLiteExtensionData()
-{
-#ifdef DEBUG
-    CPLFree(pDummy);
-#endif
-
-    std::map< std::pair<int,int>, OGRCoordinateTransformation*>::iterator oIter =
-        oCachedTransformsMap.begin();
-    for(; oIter != oCachedTransformsMap.end(); ++oIter)
-        delete oIter->second;
-
-    OGRSQLiteFreeRegExpCache(hRegExpCache);
-
-    OGRGeocodeDestroySession(hGeocodingSession);
-}
+#define DEFINE_OGRSQLiteExtensionData_GetTransform
+#include "ogrsqlitesqlfunctionscommon.cpp"
 
 /************************************************************************/
-/*                          GetTransform()                              */
-/************************************************************************/
-
-OGRCoordinateTransformation* OGRSQLiteExtensionData::GetTransform(int nSrcSRSId,
-                                                                  int nDstSRSId)
-{
-    std::map< std::pair<int,int>, OGRCoordinateTransformation*>::iterator oIter =
-        oCachedTransformsMap.find(std::pair<int,int>(nSrcSRSId, nDstSRSId));
-    if( oIter == oCachedTransformsMap.end() )
-    {
-        OGRCoordinateTransformation* poCT = nullptr;
-        OGRSpatialReference oSrcSRS, oDstSRS;
-        oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
-        oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
-        if (oSrcSRS.importFromEPSG(nSrcSRSId) == OGRERR_NONE &&
-            oDstSRS.importFromEPSG(nDstSRSId) == OGRERR_NONE )
-        {
-            poCT = OGRCreateCoordinateTransformation( &oSrcSRS, &oDstSRS );
-        }
-        oCachedTransformsMap[std::pair<int,int>(nSrcSRSId, nDstSRSId)] = poCT;
-        return poCT;
-    }
-    else
-        return oIter->second;
-}
-
-/************************************************************************/
-/*                        OGR2SQLITE_ogr_version()                     */
+/*                        OGR2SQLITE_ogr_version()                      */
 /************************************************************************/
 
 static
@@ -1108,7 +1028,7 @@ void OGRSQLITE_hstore_get_value(sqlite3_context* pContext,
 static
 void* OGRSQLiteRegisterSQLFunctions(sqlite3* hDB)
 {
-    OGRSQLiteExtensionData* pData = new OGRSQLiteExtensionData(hDB);
+    OGRSQLiteExtensionData* pData = OGRSQLiteRegisterSQLFunctionsCommon(hDB);
 
     sqlite3_create_function(hDB, "ogr_version", 0,
                             SQLITE_UTF8 | SQLITE_DETERMINISTIC, nullptr,
@@ -1272,18 +1192,5 @@ void* OGRSQLiteRegisterSQLFunctions(sqlite3* hDB)
         }
     }
 
-    pData->SetRegExpCache(OGRSQLiteRegisterRegExpFunction(hDB));
-
     return pData;
 }
-
-/************************************************************************/
-/*                   OGRSQLiteUnregisterSQLFunctions()                  */
-/************************************************************************/
-
-static
-void OGRSQLiteUnregisterSQLFunctions(void* hHandle)
-{
-    OGRSQLiteExtensionData* pData = (OGRSQLiteExtensionData* )hHandle;
-    delete pData;
-}