From 84e97628ed1fbc9116adc2d44627a336271d8086 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Mon, 27 Jul 2026 14:31:54 -0600 Subject: [PATCH] correct WKT when exporting from a XY location --- README.md | 10 +++++++++- src/xarray_grass/grass_interface.py | 14 +++++++++++++- tests/test_grass_interface.py | 12 ++++++++++++ tests/test_xarray_grass.py | 3 +++ 4 files changed, 37 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index fcd04aa..4322c4f 100644 --- a/README.md +++ b/README.md @@ -95,11 +95,19 @@ The attributes of the coordinates are in line with CF Conventions. ### Dataset attributes The attributes set at the dataset level are: - - `crs_wkt` from the g.proj command + - `crs_wkt` from the `g.proj` command - `Conventions`, the CF Convention version - `history`, the time of creation and version of xarray-grass - `source`, the name of the current grass project and mapset +For a GRASS XY (unprojected) location, `g.proj` does not provide a CRS and +returns the display label `XY location (unprojected)`. xarray-grass represents +such a location as a WKT2 engineering CRS named `XY location (unprojected)`, +with a Cartesian east/north coordinate system measured in metres. This keeps +the location explicitly local rather than assigning a geographic or projected +CRS, while ensuring that `crs_wkt` can be parsed by WKT consumers such as +`pyproj.CRS.from_wkt`. + ## Writing an Xarray Dataset or DataArray to GRASS Continuing the script fro above, we can now write back the STRDS to GRASS. diff --git a/src/xarray_grass/grass_interface.py b/src/xarray_grass/grass_interface.py index 7174c8b..7b8c18b 100644 --- a/src/xarray_grass/grass_interface.py +++ b/src/xarray_grass/grass_interface.py @@ -40,6 +40,15 @@ gs.core.set_raise_on_error(True) +XY_LOCATION_CRS_WKT = ( + 'ENGCRS["XY location (unprojected)",' + 'EDATUM["Unknown engineering datum"],' + "CS[Cartesian,2]," + 'AXIS["x",east,ORDER[1],LENGTHUNIT["metre",1]],' + 'AXIS["y",north,ORDER[2],LENGTHUNIT["metre",1]]]' +) + + @dataclass class GrassConfig: gisdb: str | Path @@ -244,7 +253,8 @@ def set_region(region_data: RegionData) -> None: def is_latlon() -> bool: return bool(gs.locn_is_latlong()) - def is_xy(self) -> bool: + @staticmethod + def is_xy() -> bool: """return True if the location is neither projected or latlon""" proj_code = gs.parse_command("g.region", flags="pug")["projection"] if int(proj_code) == 0: @@ -312,6 +322,8 @@ def name_is_raster_3d(self, raster3d_name: str) -> bool: @staticmethod def get_crs_wkt_str() -> str: + if GrassInterface.is_xy(): + return XY_LOCATION_CRS_WKT return gs.read_command("g.proj", flags="wf").replace("\n", "") def grass_dtype(self, dtype: str) -> str: diff --git a/tests/test_grass_interface.py b/tests/test_grass_interface.py index 677a4d2..a1d8a6d 100644 --- a/tests/test_grass_interface.py +++ b/tests/test_grass_interface.py @@ -24,6 +24,7 @@ import grass_session # noqa: F401 import grass.script as gs import grass.exceptions as gexceptions +from pyproj import CRS from xarray_grass import GrassInterface, RegionData @@ -186,6 +187,17 @@ def test_get_crs_wkt_str(self): ref_str = gs.read_command("g.proj", flags="wf") assert crs_str == ref_str.replace("\n", "") assert isinstance(crs_str, str) + assert CRS.from_wkt(crs_str).equals(CRS.from_wkt(ref_str)) + + def test_get_crs_wkt_str_for_xy_location(self, monkeypatch): + monkeypatch.setattr(GrassInterface, "is_xy", staticmethod(lambda: True)) + + crs = CRS.from_wkt(GrassInterface.get_crs_wkt_str()) + + assert crs.name == "XY location (unprojected)" + assert crs.is_engineering + assert [axis.direction for axis in crs.axis_info] == ["east", "north"] + assert all(axis.unit_name == "metre" for axis in crs.axis_info) def test_has_mask(self): assert GrassInterface.has_mask() is False diff --git a/tests/test_xarray_grass.py b/tests/test_xarray_grass.py index ddfb548..3b9caa2 100644 --- a/tests/test_xarray_grass.py +++ b/tests/test_xarray_grass.py @@ -18,6 +18,7 @@ import pytest import xarray as xr +from pyproj import CRS from xarray_grass.xarray_grass import dir_is_grass_mapset from xarray_grass.xarray_grass import dir_is_grass_project @@ -304,6 +305,8 @@ def test_attributes_separation(self, grass_i, temp_gisdb) -> None: assert test_dataset.attrs["Conventions"] == "CF-1.13-draft" assert "crs_wkt" in test_dataset.attrs assert isinstance(test_dataset.attrs["crs_wkt"], str) + exported_crs = CRS.from_wkt(test_dataset.attrs["crs_wkt"]) + assert exported_crs.equals(CRS.from_wkt(grass_i.get_crs_wkt_str())) # DataArray-level attributes that should NOT appear at Dataset level dataarray_only_attrs = {"long_name", "units", "comment"}