correct WKT when exporting from a XY location

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.

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:

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

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"}