From a9b26d35b823cf1b2f59d4cea96909e1b11f4128 Mon Sep 17 00:00:00 2001 From: Nicholas YS Tan <39729063+nicholas-ys-tan@users.noreply.github.com> Date: Sat, 15 Jun 2024 02:28:24 +1000 Subject: [PATCH] BUG: fix write of kml lon/lat transpose (#421) --- CHANGES.md | 2 + pyogrio/_io.pyx | 6 ++- pyogrio/_ogr.pxd | 5 +- pyogrio/tests/test_geopandas_io.py | 76 +++++++++++++++++++++++++++++- 4 files changed, 86 insertions(+), 3 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index d0967fdd..8a63f95d 100644 --- a/CHANGES.md +++ b/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 diff --git a/pyogrio/_io.pyx b/pyogrio/_io.pyx index d4430b5d..de8de43f 100644 --- a/pyogrio/_io.pyx +++ b/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}." - ) + ) \ No newline at end of file diff --git a/pyogrio/_ogr.pxd b/pyogrio/_ogr.pxd index 9c53c5c8..9369ba71 100644 --- a/pyogrio/_ogr.pxd +++ b/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) diff --git a/pyogrio/tests/test_geopandas_io.py b/pyogrio/tests/test_geopandas_io.py index f782a300..9f788cf9 100644 --- a/pyogrio/tests/test_geopandas_io.py +++ b/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)