[gdal-dev] UTM local correction for height and distance-to-reference-meridian ?

[Rahkonen Jukka]

Hi,

I was just testing the approach with ogrinfo (GDAL 3.11.0dev) but maybe I am doing something wrong.

This works in SpatiaLite gui (or at least the length is not the same):

Cartesian length
select ST_Length
(ST_GeomFromText
('LINESTRING ( 339775 6936992, 548057 6951815 )',3067));

208808.794961

Length along ellipsoid (accepts only Long/Lat coordinates, see https://www.gaia-gis.it/gaia-sins/spatialite-sql-5.1.0.html:

select ST_Length
(ST_Transform(
ST_GeomFromText
('LINESTRING ( 339775 6936992, 548057 6951815 )',3067),4326)
,1);

208875.050186

Now the same with ogrinfo. I simplified the SQL by  saving the test linesting into a file in EPSG:3067.

ogrinfo -dialect sqlite -sql "select ST_Length(geometry) from length" length.json
208808.794960845

Good so far. But with "use_ellipsoid" gives very unexpected result:

ogrinfo -dialect sqlite -sql "select ST_Length(geometry,1) from length" length.json

417750.100376289

Same wrong result also with this:

ogrinfo -dialect sqlite -sql "select ST_Length(ST_Transform(geometry,4326),1) from length" length.json

length.json is this:

{
"type": "FeatureCollection",
"crs": {"type": "name", "properties": { "name": "EPSG:3067" }},
"features": [
{ "type": "Feature", "properties": null, "geometry": {"type":"LineString","coordinates":[[339775,6936992],[548057,6951815]]} }
]
}

-Jukka Rahkonen-

Lähettäjä: gdal-dev <gdal-dev-bounces at lists.osgeo.org> Puolesta Even Rouault via gdal-dev
Lähetetty: tiistai 21. tammikuuta 2025 12.45
Vastaanottaja: Thomas Knudsen <knudsen.thomas at gmail.com>; Peter Bennewitz <geo25 at pab-opto.de>
Kopio: gdal-dev at lists.osgeo.org
Aihe: Re: [gdal-dev] UTM local correction for height and distance-to-reference-meridian ?



Le 21/01/2025 à 11:37, Thomas Knudsen via gdal-dev a écrit :
If you plan to transform to a different system anyway, I would definitely prefer going to geographical coordinates, and computing the length of the geodesic between the two points.

actually no need to do the explicit conversion to geographical coordinates, if the appropriate SRS is attached to your geometry, OGRGeometry::get_GeodesicLength() will do the right thing: https://gdal.org/en/stable/api/ogrgeometry_cpp.html#_CPPv4NK13OGRLineString18get_GeodesicLengthEPK19OGRSpatialReference

Python API: https://gdal.org/en/stable/api/python/vector_api.html#osgeo.ogr.Geometry.GeodesicLength

Python examples: https://github.com/OSGeo/gdal/blob/master/autotest/ogr/ogr_geom.py#L4560

--

http://www.spatialys.com<http://www.spatialys.com/>

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20250121/37519d5c/attachment-0001.htm>