SQLite/GPKG: Add ST_Length(geom, use_ellipsoid)

autotest/ogr/ogr_gpkg.py

@@ -10541,6 +10541,88 @@ def test_ogr_gpkg_ST_Area_on_ellipsoid(tmp_vsimem):
             assert f[0] is None
 
 
+###############################################################################
+# Test ST_Length(geom)
+
+
+def test_ogr_gpkg_ST_Length(tmp_vsimem):
+
+    tmpfilename = tmp_vsimem / "test_ogr_sql_ST_Length.gpkg"
+
+    ds = ogr.GetDriverByName("GPKG").CreateDataSource(tmpfilename)
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4258)
+    lyr = ds.CreateLayer("my_layer", srs=srs)
+    geom_colname = lyr.GetGeometryColumn()
+    feat = ogr.Feature(lyr.GetLayerDefn())
+    feat.SetGeometryDirectly(
+        ogr.CreateGeometryFromWkt("POLYGON((2 49,3 49,3 48,2 49))")
+    )
+    lyr.CreateFeature(feat)
+    feat = None
+
+    with ds.ExecuteSQL(f"SELECT ST_Length({geom_colname}) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] == pytest.approx(3.414213562373095)
+
+    with ds.ExecuteSQL("SELECT ST_Length(null) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] is None
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL("SELECT ST_Length(X'FF') FROM my_layer") as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] is None
+
+
+###############################################################################
+# Test ST_Length(geom, use_ellipsoid=True)
+
+
+def test_ogr_gpkg_ST_Length_on_ellipsoid(tmp_vsimem):
+
+    tmpfilename = tmp_vsimem / "test_ogr_sql_ST_Length_on_ellipsoid.gpkg"
+
+    ds = ogr.GetDriverByName("GPKG").CreateDataSource(tmpfilename)
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4258)
+    lyr = ds.CreateLayer("my_layer", srs=srs)
+    geom_colname = lyr.GetGeometryColumn()
+    feat = ogr.Feature(lyr.GetLayerDefn())
+    feat.SetGeometryDirectly(
+        ogr.CreateGeometryFromWkt("POLYGON((2 49,3 49,3 48,2 49))")
+    )
+    lyr.CreateFeature(feat)
+    feat = None
+
+    with ds.ExecuteSQL(f"SELECT ST_Length({geom_colname}, 1) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] == pytest.approx(317885.7863996293)
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL(
+            f"SELECT ST_Length({geom_colname}, 0) FROM my_layer"
+        ) as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] == pytest.approx(317885.7863996293)
+
+    with ds.ExecuteSQL("SELECT ST_Length(null, 1) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] is None
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL("SELECT ST_Length(X'FF', 1) FROM my_layer") as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] is None
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL(
+            f"SELECT ST_Length(SetSRID({geom_colname}, -10), 0) FROM my_layer"
+        ) as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] is None
+
+
 ###############################################################################
 # Test LAUNDER=YES layer creation option
 

autotest/ogr/ogr_sqlite.py

@@ -4094,6 +4094,49 @@ def test_ogr_sql_ST_Area_on_ellipsoid(tmp_vsimem, require_spatialite):
         assert f[0] is None
 
 
+###############################################################################
+# Test ST_Length(geom, use_ellipsoid=True)
+
+
+def test_ogr_sql_ST_Length_on_ellipsoid(tmp_vsimem):
+
+    tmpfilename = tmp_vsimem / "test_ogr_sql_ST_Length_on_ellipsoid.db"
+
+    ds = ogr.GetDriverByName("SQLite").CreateDataSource(
+        tmpfilename, options=["SPATIALITE=YES"]
+    )
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4258)
+    lyr = ds.CreateLayer("my_layer", srs=srs)
+    geom_colname = lyr.GetGeometryColumn()
+    feat = ogr.Feature(lyr.GetLayerDefn())
+    feat.SetGeometryDirectly(
+        ogr.CreateGeometryFromWkt("LINESTRING(2 49,3 49,3 48,2 49)")
+    )
+    lyr.CreateFeature(feat)
+    feat = None
+
+    with ds.ExecuteSQL(f"SELECT ST_Length({geom_colname}, 1) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] == pytest.approx(317885.7863996293)
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL(
+            f"SELECT ST_Length({geom_colname}, 0) FROM my_layer"
+        ) as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] == pytest.approx(317885.7863996293)
+
+    with ds.ExecuteSQL("SELECT ST_Length(null, 1) FROM my_layer") as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f[0] is None
+
+    with gdal.quiet_errors():
+        with ds.ExecuteSQL("SELECT ST_Length(X'FF', 1) FROM my_layer") as sql_lyr:
+            f = sql_lyr.GetNextFeature()
+            assert f[0] is None
+
+
 def test_ogr_sqlite_stddev():
     """Test STDDEV_POP() and STDDEV_SAMP"""
 

doc/source/drivers/vector/gpkg.rst

@@ -123,6 +123,10 @@ Spatialite, are also available :
 -  ST_Area(geom *Geometry*, use_ellipsoid *boolean*): (GDAL >= 3.9): compute
    the area in square meters, considering the geometry on the ellipsoid
    (use_ellipsoid must be set to true/1).
+-  ST_Length(geom *Geometry*): (GDAL >= 3.10): compute the area in units of the geometry SRS.
+-  ST_Length(geom *Geometry*, use_ellipsoid *boolean*): (GDAL >= 3.10): compute
+   the area in meters, considering the geometry on the ellipsoid
+   (use_ellipsoid must be set to true/1).
 -  SetSRID(geom *Geometry*, srs_id *Integer*): overrides the geometry' SRS ID,
    without reprojection.
 -  ST_Transform(geom *Geometry*, target_srs_id *Integer*): reproject the geometry

ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp

@@ -8750,6 +8750,77 @@ static void OGRGeoPackageGeodesicArea(sqlite3_context *pContext, int argc,
         pContext, OGR_G_GeodesicArea(OGRGeometry::ToHandle(poGeom.get())));
 }
 
+/************************************************************************/
+/*                   OGRGeoPackageLengthOrGeodesicLength()              */
+/************************************************************************/
+
+static void OGRGeoPackageLengthOrGeodesicLength(sqlite3_context *pContext,
+                                                int argc, sqlite3_value **argv)
+{
+    if (sqlite3_value_type(argv[0]) != SQLITE_BLOB)
+    {
+        sqlite3_result_null(pContext);
+        return;
+    }
+    if (argc == 2 && sqlite3_value_int(argv[1]) != 1)
+    {
+        CPLError(CE_Warning, CPLE_NotSupported,
+                 "ST_Length(geom, use_ellipsoid) is only supported for "
+                 "use_ellipsoid = 1");
+    }
+
+    const int nBLOBLen = sqlite3_value_bytes(argv[0]);
+    const GByte *pabyBLOB =
+        reinterpret_cast<const GByte *>(sqlite3_value_blob(argv[0]));
+    GPkgHeader sHeader;
+    if (!OGRGeoPackageGetHeader(pContext, argc, argv, &sHeader, false, false))
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Invalid geometry");
+        sqlite3_result_blob(pContext, nullptr, 0, nullptr);
+        return;
+    }
+
+    GDALGeoPackageDataset *poDS =
+        static_cast<GDALGeoPackageDataset *>(sqlite3_user_data(pContext));
+
+    OGRSpatialReference *poSrcSRS = nullptr;
+    if (argc == 2)
+    {
+        poSrcSRS = poDS->GetSpatialRef(sHeader.iSrsId, true);
+        if (!poSrcSRS)
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "SRID set on geometry (%d) is invalid", sHeader.iSrsId);
+            sqlite3_result_blob(pContext, nullptr, 0, nullptr);
+            return;
+        }
+    }
+
+    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_blob(pContext, nullptr, 0, nullptr);
+            return;
+        }
+        poGeom.reset(poGeomSpatialite);
+    }
+
+    if (argc == 2)
+        poGeom->assignSpatialReference(poSrcSRS);
+
+    sqlite3_result_double(
+        pContext,
+        argc == 1 ? OGR_G_Length(OGRGeometry::ToHandle(poGeom.get()))
+                  : OGR_G_GeodesicLength(OGRGeometry::ToHandle(poGeom.get())));
+}
+
 /************************************************************************/
 /*                      OGRGeoPackageTransform()                        */
 /************************************************************************/
@@ -9452,6 +9523,13 @@ void GDALGeoPackageDataset::InstallSQLFunctions()
                                 OGRGeoPackageSTMakeValid, nullptr, nullptr);
     }
 
+    sqlite3_create_function(hDB, "ST_Length", 1, UTF8_INNOCUOUS, nullptr,
+                            OGRGeoPackageLengthOrGeodesicLength, nullptr,
+                            nullptr);
+    sqlite3_create_function(hDB, "ST_Length", 2, UTF8_INNOCUOUS, this,
+                            OGRGeoPackageLengthOrGeodesicLength, nullptr,
+                            nullptr);
+
     sqlite3_create_function(hDB, "ST_Area", 1, UTF8_INNOCUOUS, nullptr,
                             OGRGeoPackageSTArea, nullptr, nullptr);
     sqlite3_create_function(hDB, "ST_Area", 2, UTF8_INNOCUOUS, this,

ogr/ogrsf_frmts/sqlite/ogrsqlitesqlfunctions.cpp

@@ -671,6 +671,52 @@ static void OGR2SQLITE_ST_GeodesicArea(sqlite3_context *pContext, int argc,
     }
 }
 
+/************************************************************************/
+/*                     OGR2SQLITE_ST_GeodesicLength()                   */
+/************************************************************************/
+
+static void OGR2SQLITE_ST_GeodesicLength(sqlite3_context *pContext, int argc,
+                                         sqlite3_value **argv)
+{
+    if (sqlite3_value_int(argv[1]) != 1)
+    {
+        CPLError(CE_Warning, CPLE_NotSupported,
+                 "ST_Length(geom, use_ellipsoid) is only supported for "
+                 "use_ellipsoid = 1");
+    }
+
+    int nSRSId = -1;
+    auto poGeom = OGR2SQLITE_GetGeom(pContext, argc, argv, &nSRSId);
+    if (poGeom != nullptr)
+    {
+        OGRSpatialReference oSRS;
+        oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+        if (nSRSId > 0)
+        {
+            if (oSRS.importFromEPSG(nSRSId) != OGRERR_NONE)
+            {
+                sqlite3_result_null(pContext);
+                return;
+            }
+        }
+        else
+        {
+            CPLDebug("OGR_SQLITE",
+                     "Assuming EPSG:4326 for GeodesicLength() computation");
+            oSRS.importFromEPSG(4326);
+        }
+        poGeom->assignSpatialReference(&oSRS);
+        sqlite3_result_double(
+            pContext,
+            OGR_G_GeodesicLength(OGRGeometry::ToHandle(poGeom.get())));
+        poGeom->assignSpatialReference(nullptr);
+    }
+    else
+    {
+        sqlite3_result_null(pContext);
+    }
+}
+
 #ifdef MINIMAL_SPATIAL_FUNCTIONS
 
 /************************************************************************/
@@ -926,6 +972,25 @@ static void OGR2SQLITE_ST_Area(sqlite3_context *pContext, int argc,
         sqlite3_result_null(pContext);
 }
 
+/************************************************************************/
+/*                      OGR2SQLITE_ST_Length()                          */
+/************************************************************************/
+
+static void OGR2SQLITE_ST_Length(sqlite3_context *pContext, int argc,
+                                 sqlite3_value **argv)
+{
+    auto poGeom = OGR2SQLITE_GetGeom(pContext, argc, argv, nullptr);
+    if (poGeom != nullptr)
+    {
+        CPLPushErrorHandler(CPLQuietErrorHandler);
+        sqlite3_result_double(
+            pContext, OGR_G_Length(OGRGeometry::ToHandle(poGeom.get())));
+        CPLPopErrorHandler();
+    }
+    else
+        sqlite3_result_null(pContext);
+}
+
 /************************************************************************/
 /*                     OGR2SQLITE_ST_Buffer()                           */
 /************************************************************************/
@@ -1148,6 +1213,7 @@ static void *OGRSQLiteRegisterSQLFunctions(sqlite3 *hDB)
 
         REGISTER_ST_op(1, SRID);
         REGISTER_ST_op(1, Area);
+        REGISTER_ST_op(1, Length);
         REGISTER_ST_op(2, Buffer);
         REGISTER_ST_op(2, MakePoint);
         REGISTER_ST_op(3, MakePoint);
@@ -1163,6 +1229,13 @@ static void *OGRSQLiteRegisterSQLFunctions(sqlite3 *hDB)
         sqlite3_create_function(hDB, "ST_Area", 2, UTF8_INNOCUOUS, nullptr,
                                 OGR2SQLITE_ST_GeodesicArea, nullptr, nullptr);
 
+        // We add a ST_Length() method with 2 arguments even when Spatialite
+        // is there to indicate we want to use the ellipsoid version
+        sqlite3_create_function(hDB, "Length", 2, UTF8_INNOCUOUS, nullptr,
+                                OGR2SQLITE_ST_GeodesicLength, nullptr, nullptr);
+        sqlite3_create_function(hDB, "ST_Length", 2, UTF8_INNOCUOUS, nullptr,
+                                OGR2SQLITE_ST_GeodesicLength, nullptr, nullptr);
+
         static bool gbRegisterMakeValid = [bSpatialiteAvailable, hDB]()
         {
             bool bRegisterMakeValid = false;