How to retrieve the projection?

[nshiab]

Hi!

St_transform converts geometry from one projection to another. But how do we know what the original projection is? Is there a function to retrieve it?

Thank you very much for your great work! :)

[ericemc3]

Just to add some specifics to that important question, could duckdb-spatial provide a function to read geoparquet files geo metadata, such as CRS, bounding box, encoding (wkb or arrow), geometry_types ?

For instance R/geoarrow can display a metadata attribute:

nc = open_dataset("myfile-geo.parquet")
> nc$metadata
$geo
[1] "{\"version\":\"1.0.0-beta.1\",\"primary_column\":\"geometry\",\"columns\":{\"geometry\":{\"encoding\":\"WKB\",\"CRS\":\"4326\",\"bbox\":[-62.850000000000001,-21.387426000000001,55.829439000000001,51.082326000000002],\"orientation\":\"counterclockwise\",\"geometry_types\":[\"Point\"]}}}"

[Maxxen]

This will be possible in the next DuckDB release, see duckdb/duckdb#9126.

[Maxxen]

It's not possible right now, but we're planning to add something like a st_read_metadata to read metadata like projection.

[Maxxen]

Seems like that is the arrow schema as a IPC blob in base64. I don't think its particularly interesting since DuckDB can give you the schema of the file by just calling DESCRIBE FROM st_read('...')

[ericemc3]

Good catch !
with the help of on online base64 decoder, It looks like there are inside interesting geographic informations:

I am not sure that an st_read can output that. And actually, st_read don't seem to read geoparquet files.

But how to decode such a blob inside DuckDB?
This:

SELECT from_base64(decode(value)) AS v
FROM parquet_kv_metadata('myfile-geo.parquet') ;

outputs:
Error: Conversion Error: Could not decode string "{"version":"1.0.0-beta.1","primary_column":"geometry","columns":{"geometry":{"encoding":"WKB","crs":{"$schema":"https:\/\/proj.org\/schemas\/v0.7\/projjson.schema.json","type":"ProjectedCRS","name":"RGF93 v1 \/ Lambert-93","base_crs":{"name":"RGF93 v1","datum":{"type":"GeodeticReferenceFrame","name":"Reseau Geodesique Francais 1993 v1","ellipsoid":{"name":"GRS 1980","semi_major_axis":6378137,"inverse_flattening":298.257222101}},"coordinate_system":{"subtype":"ellipsoidal","axis":[{"name":"Geodetic latitude","abbreviation":"Lat","direction":"north","unit":"degree"},{"name":"Geodetic longitude","abbreviation":"Lon","direction":"east","unit":"degree"}]},"id":{"authority":"EPSG","code":4171}},"conversion":{"name":"Lambert-93","method":{"name":"Lambert Conic Conformal (2SP)","id":{"authority":"EPSG","code":9802}},"parameters":[{"name":"Latitude of false origin","value":46.5,"unit":"degree","id":{"authority":"EPSG","code":8821}},{"name":"Longitude of false origin","value":3,"unit":"degree","id":{"authority":"EPSG","code":8822}},{"name":"Latitude of 1st standard parallel","value":49,"unit":"degree","id":{"authority":"EPSG","code":8823}},{"name":"Latitude of 2nd standard parallel","value":44,"unit":"degree","id":{"authority":"EPSG","code":8824}},{"name":"Easting at false origin","value":700000,"unit":"metre","id":{"authority":"EPSG","code":8826}},{"name":"Northing at false origin","value":6600000,"unit":"metre","id":{"authority":"EPSG","code":8827}}]},"coordinate_system":{"subtype":"Cartesian","axis":[{"name":"Easting","abbreviation":"X","direction":"east","unit":"metre"},{"name":"Northing","abbreviation":"Y","direction":"north","unit":"metre"}]},"scope":"Engineering survey, topographic mapping.","area":"France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica).","bbox":{"south_latitude":41.15,"west_longitude":-9.86,"north_latitude":51.56,"east_longitude":10.38},"id":{"authority":"EPSG","code":2154}},"bbox":[-13111921.893685382,471170.27157718595,14452860.825010845,20156810.236127403],"orientation":"counterclockwise","geometry_types":["Point"]}}}" as base64: length must be a multiple of 4

[Maxxen]

Yes, st_read does not read geoparquet. This discussion is not really related to the original - but if you could share an example geoparquet file I could maybe try to help you.

[ericemc3]

That is very kind of you, thank you.
Here it is :

SELECT from_base64(decode(value)) AS v
FROM parquet_kv_metadata('https://static.data.gouv.fr/resources/irve-formt-geoparquet-20231222/20240109-200540/irve-20231222-qgis-geo.parquet') ;

This is about electric recharging stations in France

[kylebarron]

The geoparquet metadata is usually stored both within the parquet metadata and as part of the arrow schema. So the arrow schema is purely redundant unless you want to recreate arrow types

[Maxxen]

I just merged #227 so this should now be possible!

[nshiab]

Awesome! Thanks a lot, @Maxxen !