BUG: long, lat geometry swapped ! geopandas.to_file( "out.kml", driver="KML",engine='pyogrio' )

[nicholas-ys-tan]

Hello, opening an issue here that was raised in geopandas by u/phisan-chula

https://github.com/geopandas/geopandas/issues/3332

because I wanted to have a go at submitting a PR on this. I assume pyogrio is the right place to put in a fix.

This seems to be an issue that came up in GDAL 3.0 when it started using the convention set by the CRS - for EPSG:4326 - that was lat/lon - in conflict with the lon/lat convention by WKT/WKB (to my understanding).

This resulted in the KML's exported inferring a lat/lon convention when writing despite the explicit lon/lat input by the geodataframe.

This was resolved by GDAL allowing the AxisMappingStrategy to be set as OAMS_TRADITIONAL_GIS_ORDER (i.e. lon/lat)

Apologies if I got any details above wrong, I am probably operating at the limit of my geospatial knowledge at the moment.

[brendan-ward]

It looks like the GDAL KML driver tries to set the OAMS_TRADITIONAL_GIS_ORDER if the incoming CRS is not null, and then always uses that . However, it seems that is not working when we pass in that CRS (not sure why yet), maybe it is considering EPSG:4326 and EPSG:4326 (with OAMS_TRADITIONAL_GIS_ORDER applied) as equal and not applying a transformation?

Looking at this makes me think that setting OAMS_TRADITIONAL_GIS_ORDER should occur within GDAL drivers instead of on our end, since GDAL is already trying to handle this as part of the driver.

I notice that ogr2ogr preserves the correct coordinate order when converting from a dataset in EPSG:4326 coordinate order to a KML.

[nicholas-ys-tan]

That's interesting, does look like it's intended to be managed on GDAL side. Noted that was implemented in June 2019 but QGIS implemented a fix in their end in Nov 2019 here: https://github.com/qgis/QGIS/commit/bd170fa7fbb6b7edaabd1212ba4b926b7a1d6cbd

[jorisvandenbossche]

Also on the Fiona side they set https://github.com/Toblerity/Fiona/blob/991d7fceeba88f72ddda2275608bf97bfcbddf8f/fiona/ogrext.pyx#L1412-L1417 the traditional GIS order on the SRS object when creating the dataset layer to write to (calling OSRSetAxisMappingStrategy with OAMS_TRADITIONAL_GIS_ORDER)

[jorisvandenbossche]

From reading https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order, it sounds like this should be set on our side. What I don't fully understand, though, is why this never came up before if that is actually an issue. Or the reason might be that for most file formats no coordinate transformation happens while writing, and so we don't notice that?

Testing that hypothesis, another example where a coordinate transformation can happen is when writing to GeoJSON and asking for the RFC7946 version of geojson (this is not the default, but when enabled in will force the result to use WGS84 (OGC::CRS84)). So if we then start from a gdf with a CRS that is not lon-lat order (or x-y order, as most projected CRS use), it will trigger a transformation that might be using a wrong order (it thinks the data comes in as lat/lon, while we pass it as lon/lat):

In [28]: gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy([50, 51], [14, 15]), crs="EPSG:4326")

In [29]: gdf.to_file("/tmp/test-axis-order.json", driver="GeoJSON", engine='pyogrio')

In [30]: gdf.to_file("/tmp/test-axis-order2.json", driver="GeoJSON", engine='pyogrio', RFC7946=True)

In [31]: !cat /tmp/test-axis-order.json
{
"type": "FeatureCollection",
"name": "test-axis-order",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 50.0, 14.0 ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 51.0, 15.0 ] } }
]
}

In [32]: !cat /tmp/test-axis-order2.json
{
"type": "FeatureCollection",
"name": "test-axis-order2",
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 14.0, 50.0 ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 15.0, 51.0 ] } }
]
}

So also here we see a wrong result ..