BUG: fix write of kml lon/lat transpose (#421)

CHANGES.md

@@ -57,6 +57,8 @@
     `write` and `write_dataframe` (#384).
 -   Fixed bug preventing reading from bytes or file-like in `read_arrow` /
     `open_arrow` (#407).
+-   Fixed bug transposing longitude and latitude when writing files with
+    coordinate transformation from EPSG:4326 (#421).
 
 ### Packaging
 

pyogrio/_io.pyx

@@ -2134,6 +2134,10 @@ cdef create_ogr_dataset_layer(
         if crs is not None:
             try:
                 ogr_crs = create_crs(crs)
+                # force geographic CRS to use lon, lat order and ignore axis order specified by CRS, in order
+                # to correctly write KML and GeoJSON coordinates in correct order
+                OSRSetAxisMappingStrategy(ogr_crs, OAMS_TRADITIONAL_GIS_ORDER)
+
 
             except Exception as exc:
                 if dataset_options != NULL:
@@ -2735,4 +2739,4 @@ cdef create_fields_from_arrow_schema(
                 f"Error while creating field from Arrow for field {i} with name "
                 f"'{get_string(child.name)}' and type {get_string(child.format)}"
                 f"{gdal_msg}."
-            )
+            )

pyogrio/_ogr.pxd

@@ -196,7 +196,10 @@ cdef extern from "ogr_srs_api.h":
     const char*             OSRGetAuthorityName(OGRSpatialReferenceH srs, const char *key)
     const char*             OSRGetAuthorityCode(OGRSpatialReferenceH srs, const char *key)
     OGRErr                  OSRImportFromEPSG(OGRSpatialReferenceH srs, int code)
-
+    ctypedef enum OSRAxisMappingStrategy:
+        OAMS_TRADITIONAL_GIS_ORDER
+
+    void                    OSRSetAxisMappingStrategy(OGRSpatialReferenceH hSRS, OSRAxisMappingStrategy)
     int                     OSRSetFromUserInput(OGRSpatialReferenceH srs, const char *pszDef)
     void                    OSRSetPROJSearchPaths(const char *const *paths)
     OGRSpatialReferenceH    OSRNewSpatialReference(const char *wkt)

pyogrio/tests/test_geopandas_io.py

@@ -6,7 +6,7 @@
 import numpy as np
 import pytest
 
-from pyogrio import list_layers, read_info, __gdal_version__
+from pyogrio import list_layers, list_drivers, read_info, __gdal_version__
 from pyogrio.errors import DataLayerError, DataSourceError, FeatureError, GeometryError
 from pyogrio.geopandas import read_dataframe, write_dataframe, PANDAS_GE_20
 from pyogrio.raw import (
@@ -2067,3 +2067,77 @@ def test_non_utf8_encoding_shapefile_sql(tmp_path, use_arrow):
     )
     assert actual.columns[0] == mandarin
     assert actual[mandarin].values[0] == mandarin
+
+
+@pytest.mark.requires_arrow_write_api
+def test_write_kml_file_coordinate_order(tmp_path, use_arrow):
+    # confirm KML coordinates are written in lon, lat order even if CRS axis specifies otherwise
+    points = [Point(10, 20), Point(30, 40), Point(50, 60)]
+    gdf = gp.GeoDataFrame(geometry=points, crs="EPSG:4326")
+    output_path = tmp_path / "test.kml"
+    write_dataframe(
+        gdf, output_path, layer="tmp_layer", driver="KML", use_arrow=use_arrow
+    )
+
+    gdf_in = read_dataframe(output_path, use_arrow=use_arrow)
+
+    assert np.array_equal(gdf_in.geometry.values, points)
+
+    if "LIBKML" in list_drivers():
+        # test appending to the existing file only if LIBKML is available
+        # as it appears to fall back on LIBKML driver when appending.
+        points_append = [Point(70, 80), Point(90, 100), Point(110, 120)]
+        gdf_append = gp.GeoDataFrame(geometry=points_append, crs="EPSG:4326")
+
+        write_dataframe(
+            gdf_append,
+            output_path,
+            layer="tmp_layer",
+            driver="KML",
+            use_arrow=use_arrow,
+            append=True,
+        )
+        # force_2d used to only compare xy geometry as z-dimension is undesirably
+        # introduced when the kml file is over-written.
+        gdf_in_appended = read_dataframe(
+            output_path, use_arrow=use_arrow, force_2d=True
+        )
+
+        assert np.array_equal(gdf_in_appended.geometry.values, points + points_append)
+
+
+@pytest.mark.requires_arrow_write_api
+def test_write_geojson_rfc7946_coordinates(tmp_path, use_arrow):
+    points = [Point(10, 20), Point(30, 40), Point(50, 60)]
+    gdf = gp.GeoDataFrame(geometry=points, crs="EPSG:4326")
+    output_path = tmp_path / "test.geojson"
+    write_dataframe(
+        gdf,
+        output_path,
+        layer="tmp_layer",
+        driver="GeoJSON",
+        RFC7946=True,
+        use_arrow=use_arrow,
+    )
+
+    gdf_in = read_dataframe(output_path, use_arrow=use_arrow)
+
+    assert np.array_equal(gdf_in.geometry.values, points)
+
+    # test appending to the existing file
+
+    points_append = [Point(70, 80), Point(90, 100), Point(110, 120)]
+    gdf_append = gp.GeoDataFrame(geometry=points_append, crs="EPSG:4326")
+
+    write_dataframe(
+        gdf_append,
+        output_path,
+        layer="tmp_layer",
+        driver="GeoJSON",
+        RFC7946=True,
+        use_arrow=use_arrow,
+        append=True,
+    )
+
+    gdf_in_appended = read_dataframe(output_path, use_arrow=use_arrow)
+    assert np.array_equal(gdf_in_appended.geometry.values, points + points_append)