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

autotest/gdrivers/gpkg.py

@@ -4216,6 +4216,118 @@ def test_gpkg_byte_nodata_value(band_count):
     gdal.Unlink(filename)
 
 
+###############################################################################
+# Test gdal_get_layer_pixel_value() function
+
+
+def test_gpkg_sql_gdal_get_layer_pixel_value():
+
+    filename = "/vsimem/test_ogr_gpkg_sql_gdal_get_layer_pixel_value.gpkg"
+    src_ds = gdal.Open("data/byte.tif")
+    gdal.GetDriverByName("GPKG").CreateCopy(
+        filename, src_ds, options=["RASTER_TABLE=byte"]
+    )
+    src_ds = gdal.Open("data/float32.tif")
+    gdal.GetDriverByName("GPKG").CreateCopy(
+        filename, src_ds, options=["APPEND_SUBDATASET=YES"]
+    )
+
+    ds = gdal.OpenEx(filename)
+    sql_lyr = ds.ExecuteSQL(
+        "select gdal_get_layer_pixel_value('byte', 1, 'georef', 440780, 3751080)"
+    )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 156
+
+    sql_lyr = ds.ExecuteSQL(
+        "select gdal_get_layer_pixel_value('byte', 1, 'pixel', 1, 4)"
+    )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 156
+
+    sql_lyr = ds.ExecuteSQL(
+        "select gdal_get_layer_pixel_value('float32', 1, 'pixel', 0, 1)"
+    )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 115.0
+
+    # Invalid column
+    sql_lyr = ds.ExecuteSQL(
+        "select gdal_get_layer_pixel_value('byte', 1, 'pixel', -1, 0)"
+    )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] is None
+
+    # NULL as 1st arg
+    with gdaltest.error_handler():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_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():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 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():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 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():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 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():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 1, 'pixel', 0, NULL)"
+        )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    # Invalid band number
+    with gdaltest.error_handler():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 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():
+        sql_lyr = ds.ExecuteSQL(
+            "select gdal_get_layer_pixel_value('byte', 1, 'invalid', 0, 0)"
+        )
+        f = sql_lyr.GetNextFeature()
+        ds.ReleaseResultSet(sql_lyr)
+        assert f[0] is None
+
+    gdal.Unlink(filename)
+
+
 ###############################################################################
 #
 

autotest/ogr/ogr_gpkg.py

@@ -8055,3 +8055,127 @@ 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', 440780, 3751080)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 156
+
+    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, 4)"
+        )
+    f = sql_lyr.GetNextFeature()
+    ds.ReleaseResultSet(sql_lyr)
+    assert f[0] == 156
+
+    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/raster/gpkg.rst

@@ -610,6 +610,33 @@ Examples
 
       gdalinfo my.gpkg -oo TABLE=a_table
 
+.. _raster.gpkg.raster:
+
+Raster SQL functions
+~~~~~~~~~~~~~~~~~~~~
+
+The raster SQL functions mentioned at :ref:`sql_sqlite_dialect_raster_functions`
+are also available.
+
+The ``gdal_get_layer_pixel_value()`` function (added in GDAL 3.7), variant of the
+generic ``gdal_get_pixel_value()``, can be used to extract the value of a pixel
+in a raster layer of the current dataset.
+
+It takes 5 arguments:
+
+* a string with the layer/table 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_layer_pixel_value('my_raster_table', 1, 'georef', 440720, 3751320)
+    SELECT gdal_get_layer_pixel_value('my_raster_table', 1, 'pixel', 0, 0)
+
 See Also
 --------
 

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:`raster.gpkg.raster`
+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,11 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource, public GDALG
         bool                m_bIsGeometryTypeAggregateInterrupted = false;
         std::string         m_osGeometryTypeAggregateResult{};
 
+        // Used by GDALGeoPackageDataset::GetRasterLayerDataset()
+        std::map< std::string, std::unique_ptr<GDALDataset> > m_oCachedRasterDS{};
+
+        void                CloseDB();
+
         CPL_DISALLOW_COPY_ASSIGN(GDALGeoPackageDataset)
 
     public:
@@ -387,6 +393,8 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource, public GDALG
         void                SetGeometryTypeAggregateResult(const std::string& s) { m_osGeometryTypeAggregateResult = s; }
         const std::string&  GetGeometryTypeAggregateResult() const { return m_osGeometryTypeAggregateResult; }
 
+        GDALDataset*        GetRasterLayerDataset(const char* pszLayerName);
+
     protected:
 
         virtual CPLErr IRasterIO( GDALRWFlag, int, int, int, int,