From d4bef0cba87d234c00b1b910169689af4d8157c7 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Thu, 30 Oct 2025 17:49:17 -0600 Subject: [PATCH 1/9] fix_moving attrs in dataset creation --- src/xarray_grass/xarray_grass.py | 50 +++++++++++------- tests/test_xarray_grass.py | 87 ++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+), 17 deletions(-) diff --git a/src/xarray_grass/xarray_grass.py b/src/xarray_grass/xarray_grass.py index 3615b6e..043cfa7 100644 --- a/src/xarray_grass/xarray_grass.py +++ b/src/xarray_grass/xarray_grass.py @@ -252,12 +252,14 @@ def open_grass_maps( # Open all given maps and identify non-existent data not_found = {config["not_found_key"]: [] for config in map_processing_configs} data_array_list = [] + raw_coords_list = [] for config in map_processing_configs: for map_name in config["input_list"]: if not config["existence_check_method"](map_name): not_found[config["not_found_key"]].append(map_name) continue data_array = config["open_function"](map_name, gi) + raw_coords_list.append(data_array.coords) data_array_list.append(data_array) if raise_on_not_found and any(not_found.values()): raise ValueError(f"Objects not found: {not_found}") @@ -265,15 +267,25 @@ def open_grass_maps( if session is not None: session.__exit__(None, None, None) - dataset = xr.merge(data_array_list) - dataset.attrs["crs_wkt"] = gi.get_crs_wkt_str() - dataset.attrs["Conventions"] = "CF-1.13-draft" - # dataset.attrs["title"] = "" - # dataset.attrs["history"] = "" - # dataset.attrs["source"] = "" - # dataset.attrs["references"] = "" - # dataset.attrs["institution"] = "" - # dataset.attrs["comment"] = "" + # TODO: time coordinates must be merged + coords_dict = {} + for i, coords in enumerate(raw_coords_list): + for k, v in coords.items(): + coords_dict[k] = v + + data_array_dict = {da.name: da for da in data_array_list} + + attrs = { + "crs_wkt": gi.get_crs_wkt_str(), + "Conventions": "CF-1.13-draft", + # "title": "", + # "history": "", + # "source": "", + # "references": "", + # "institution": "", + # "comment": "", + } + dataset = xr.Dataset(data_vars=data_array_dict, coords=coords_dict, attrs=attrs) return dataset @@ -370,7 +382,11 @@ def open_grass_strds(strds_name: str, grass_i: GrassInterface) -> xr.DataArray: # Add CF attributes r_infos = grass_i.get_raster_info(map_list[0].id) da_with_attrs = set_cf_coordinates( - da_concat, gi=grass_i, is_3d=False, time_dim=start_time_dim, time_unit=time_unit + da_concat, + gi=grass_i, + is_3d=False, + time_dims=[start_time_dim, end_time_dim], + time_unit=time_unit, ) da_with_attrs.attrs["long_name"] = strds_infos.title da_with_attrs.attrs["source"] = ",".join([r_infos["source1"], r_infos["source2"]]) @@ -423,7 +439,7 @@ def open_grass_str3ds(str3ds_name: str, grass_i: GrassInterface) -> xr.DataArray gi=grass_i, is_3d=True, z_unit=r3_infos["vertical_units"], - time_dim=start_time_dim, + time_dims=[start_time_dim, end_time_dim], time_unit=time_unit, ) da_with_attrs.attrs["long_name"] = strds_infos.title @@ -440,7 +456,7 @@ def set_cf_coordinates( gi: GrassInterface, is_3d: bool, z_unit: str = "", - time_dim: str = "", + time_dims: Optional[list[str, str]] = None, time_unit: str = "", ): """Set coordinate attributes according to CF conventions""" @@ -454,11 +470,11 @@ def set_cf_coordinates( else: y_coord = "y" x_coord = "x" - if time_dim: - da[time_dim].attrs["axis"] = "T" - da[time_dim].attrs["standard_name"] = "time" - if time_unit: - da[time_dim].attrs["units"] = time_unit + if time_dims is not None: + for time_dim in time_dims: + da[time_dim].attrs["axis"] = "T" + da[time_dim].attrs["standard_name"] = "time" + da[time_dim].attrs["units"] = time_unit da[x_coord].attrs["axis"] = "X" da[y_coord].attrs["axis"] = "Y" if gi.is_latlon(): diff --git a/tests/test_xarray_grass.py b/tests/test_xarray_grass.py index 6f42a08..c8ca1a1 100644 --- a/tests/test_xarray_grass.py +++ b/tests/test_xarray_grass.py @@ -252,4 +252,91 @@ def test_drop_variables(self, grass_i, temp_gisdb) -> None: assert len(test_dataset.y_3d) == region.rows3 assert len(test_dataset.x) == region.cols assert len(test_dataset.y) == region.rows + + def test_attributes_separation(self, grass_i, temp_gisdb) -> None: + """Test that DataArray attributes don't leak to Dataset level.""" + mapset_path = os.path.join( + str(temp_gisdb.gisdb), str(temp_gisdb.project), str(temp_gisdb.mapset) + ) + region = grass_i.get_region() + + # Test with multiple types of data to ensure attributes are handled correctly + test_dataset = xr.open_dataset( + mapset_path, + raster=ACTUAL_RASTER_MAP, + raster_3d=ACTUAL_RASTER3D_MAP, + strds=ACTUAL_STRDS, + ) + + # Dataset-level attributes: only these should be present + expected_dataset_attrs = {"crs_wkt", "Conventions"} + actual_dataset_attrs = set(test_dataset.attrs.keys()) + assert actual_dataset_attrs == expected_dataset_attrs, ( + f"Dataset has unexpected attributes. " + f"Expected: {expected_dataset_attrs}, " + f"Got: {actual_dataset_attrs}" + ) + + # Verify Dataset attrs have correct values + assert "Conventions" in test_dataset.attrs + assert test_dataset.attrs["Conventions"] == "CF-1.13-draft" + assert "crs_wkt" in test_dataset.attrs + assert isinstance(test_dataset.attrs["crs_wkt"], str) + + # DataArray-level attributes that should NOT appear at Dataset level + dataarray_only_attrs = {"long_name", "source", "units", "comment"} + + # Check that DataArray attributes don't leak to Dataset level + for attr in dataarray_only_attrs: + assert attr not in test_dataset.attrs, ( + f"DataArray attribute '{attr}' should not appear at Dataset level" + ) + + # Verify each DataArray has its own attributes + for var_name in test_dataset.data_vars: + data_array = test_dataset[var_name] + + # Each DataArray should have at least some of these attributes + # (depending on the data source, some might be empty strings) + for attr in ["long_name", "source", "units"]: + assert attr in data_array.attrs, ( + f"DataArray '{var_name}' missing expected attribute '{attr}'" + ) + + # Verify coordinate attributes are set correctly (these are set by set_cf_coordinates) + # Check x coordinate + assert "axis" in test_dataset.x.attrs + assert test_dataset.x.attrs["axis"] == "X" + + # Check y coordinate + assert "axis" in test_dataset.y.attrs + assert test_dataset.y.attrs["axis"] == "Y" + + # Check z coordinate for 3D data + if "z" in test_dataset.coords: + assert "axis" in test_dataset.z.attrs + assert test_dataset.z.attrs["axis"] == "Z" + assert "positive" in test_dataset.z.attrs + assert test_dataset.z.attrs["positive"] == "up" + assert len(test_dataset.z) == region.depths + + # Check time coordinate attributes (from STRDS) + time_coords = [coord for coord in test_dataset.coords if "time" in coord] + assert len(time_coords) > 0, "Expected at least one time coordinate from STRDS" + + for time_coord in time_coords: + # Time coordinates should have axis="T" + assert "axis" in test_dataset[time_coord].attrs + assert test_dataset[time_coord].attrs["axis"] == "T" + + # Time coordinates should have standard_name="time" + assert "standard_name" in test_dataset[time_coord].attrs + assert test_dataset[time_coord].attrs["standard_name"] == "time" + + # For relative time coordinates, units should be present + if time_coord.startswith("start_time_") or time_coord.startswith( + "end_time_" + ): + # This is a relative time coordinate, should have units + assert "units" in test_dataset[time_coord].attrs assert len(test_dataset.z) == region.depths From 800ef85f3e3ed61dfdcf6b8455f722306df658c8 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Thu, 30 Oct 2025 17:58:59 -0600 Subject: [PATCH 2/9] add cleanup to test_region_fixture --- tests/conftest.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/conftest.py b/tests/conftest.py index 81010dc..33eb09f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -110,7 +110,15 @@ def grass_session_fixture(temp_gisdb: GrassConfig): @pytest.fixture(scope="function") def grass_test_region(grass_session_fixture): """A test region fixture to make sure the tests are run in a predictable region""" + # Save current region + gs.run_command("g.region", save="test_backup_region") + # Set test region gs.run_command("g.region", flags="o", b=0, t=1500, res3=500, res=500) + yield + # Restore original region + gs.run_command("g.region", region="test_backup_region") + # Clean up the saved region + gs.run_command("g.remove", type="region", name="test_backup_region", flags="f") @pytest.fixture(scope="class") From dd78808642e628d4bb93060d6d968f64806e132c Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 11:20:23 -0600 Subject: [PATCH 3/9] ci: use apt-get instead of apt --- .github/workflows/tests.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 1ae3933..b42dbf3 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -17,8 +17,8 @@ jobs: steps: - name: Install dependencies run: | - apt update - apt install -y ca-certificates + apt-get update + apt-get install -y ca-certificates - uses: actions/checkout@v5 From 81397fc3fa8aab0c969272a86c43b63e72f70c1e Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 11:22:16 -0600 Subject: [PATCH 4/9] each stds has its own time dimension --- src/xarray_grass/xarray_grass.py | 39 ++++----- tests/conftest.py | 135 ++++++++++++++++++++++++++++++- tests/test_grass_interface.py | 8 +- tests/test_xarray_grass.py | 26 +++--- 4 files changed, 176 insertions(+), 32 deletions(-) diff --git a/src/xarray_grass/xarray_grass.py b/src/xarray_grass/xarray_grass.py index 043cfa7..3b148f7 100644 --- a/src/xarray_grass/xarray_grass.py +++ b/src/xarray_grass/xarray_grass.py @@ -14,6 +14,7 @@ """ import os +from datetime import datetime, timezone from pathlib import Path from typing import Iterable, Optional @@ -22,6 +23,7 @@ from xarray.backends import BackendArray import xarray as xr import grass_session # noqa: F401 +import xarray_grass from xarray_grass.grass_interface import GrassInterface @@ -267,9 +269,8 @@ def open_grass_maps( if session is not None: session.__exit__(None, None, None) - # TODO: time coordinates must be merged coords_dict = {} - for i, coords in enumerate(raw_coords_list): + for coords in raw_coords_list: for k, v in coords.items(): coords_dict[k] = v @@ -279,8 +280,8 @@ def open_grass_maps( "crs_wkt": gi.get_crs_wkt_str(), "Conventions": "CF-1.13-draft", # "title": "", - # "history": "", - # "source": "", + "history": f"{datetime.now(timezone.utc)}: Created with xarray-grass version {xarray_grass.__version__}", + "source": f"GRASS database. project: <{Path(project).name}>, mapset:<{Path(mapset).name}>", # "references": "", # "institution": "", # "comment": "", @@ -349,21 +350,21 @@ def open_grass_strds(strds_name: str, grass_i: GrassInterface) -> xr.DataArray: """must be called from within a grass session TODO: lazy loading """ + strds_id = grass_i.get_id_from_name(strds_name) + strds_name = grass_i.get_name_from_id(strds_id) x_coords, y_coords, _ = get_coordinates(grass_i, raster_3d=False).values() - strds_infos = grass_i.get_stds_infos(strds_name, stds_type="strds") + strds_infos = grass_i.get_stds_infos(strds_id, stds_type="strds") if strds_infos.temporal_type == "absolute": - start_time_dim = "start_time" - end_time_dim = "end_time" time_unit = None else: time_unit = strds_infos.time_unit - start_time_dim = f"start_time_{time_unit}" - end_time_dim = f"end_time_{time_unit}" + start_time_dim = f"start_time_{strds_name}" + end_time_dim = f"end_time_{strds_name}" dims = [start_time_dim, "y", "x"] coordinates = dict.fromkeys(dims) coordinates["x"] = x_coords coordinates["y"] = y_coords - map_list = grass_i.list_maps_in_strds(strds_name) + map_list = grass_i.list_maps_in_strds(strds_id) array_list = [] for map_data in map_list: coordinates[start_time_dim] = [map_data.start_time] @@ -375,7 +376,7 @@ def open_grass_strds(strds_name: str, grass_i: GrassInterface) -> xr.DataArray: ndarray, coords=coordinates, dims=dims, - name=grass_i.get_name_from_id(strds_name), + name=strds_name, ) array_list.append(data_array) da_concat = xr.concat(array_list, dim=start_time_dim) @@ -400,22 +401,22 @@ def open_grass_strds(strds_name: str, grass_i: GrassInterface) -> xr.DataArray: def open_grass_str3ds(str3ds_name: str, grass_i: GrassInterface) -> xr.DataArray: """Open a series of 3D raster maps. TODO: Figure out what to do when the z value of the maps is time.""" + str3ds_id = grass_i.get_id_from_name(str3ds_name) + str3ds_name = grass_i.get_name_from_id(str3ds_id) x_coords, y_coords, z_coords = get_coordinates(grass_i, raster_3d=True).values() - strds_infos = grass_i.get_stds_infos(str3ds_name, stds_type="str3ds") + strds_infos = grass_i.get_stds_infos(str3ds_id, stds_type="str3ds") if strds_infos.temporal_type == "absolute": - start_time_dim = "start_time" - end_time_dim = "end_time" time_unit = None else: time_unit = strds_infos.time_unit - start_time_dim = f"start_time_{time_unit}" - end_time_dim = f"end_time_{time_unit}" + start_time_dim = f"start_time_{str3ds_name}" + end_time_dim = f"end_time_{str3ds_name}" dims = [start_time_dim, "z", "y_3d", "x_3d"] coordinates = dict.fromkeys(dims) coordinates["x_3d"] = x_coords coordinates["y_3d"] = y_coords coordinates["z"] = z_coords - map_list = grass_i.list_maps_in_str3ds(str3ds_name) + map_list = grass_i.list_maps_in_str3ds(str3ds_id) array_list = [] for map_data in map_list: coordinates[start_time_dim] = [map_data.start_time] @@ -427,7 +428,7 @@ def open_grass_str3ds(str3ds_name: str, grass_i: GrassInterface) -> xr.DataArray ndarray, coords=coordinates, dims=dims, - name=grass_i.get_name_from_id(str3ds_name), + name=str3ds_name, ) array_list.append(data_array) @@ -462,7 +463,7 @@ def set_cf_coordinates( """Set coordinate attributes according to CF conventions""" spatial_unit = gi.get_spatial_units() if is_3d: - da["z"].attrs["positive"] = "up" + # da["z"].attrs["positive"] = "up" # Not defined by grass da["z"].attrs["axis"] = "Z" da["z"].attrs["units"] = z_unit y_coord = "y_3d" diff --git a/tests/conftest.py b/tests/conftest.py index 33eb09f..ceb19af 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -103,6 +103,39 @@ def grass_session_fixture(temp_gisdb: GrassConfig): gen_str3ds(temporal_type="relative") gen_str3ds(temporal_type="absolute") + # Add custom str3ds with specific absolute times + gen_str3ds_custom_absolute( + start_times=[ + "2014-10-01 00:10:10", + "2015-10-01 01:11:59", + "2017-06-01 18:13:56", + ], + end_times=[ + "2015-09-01 05:11:23", + "2017-06-01 18:13:56", + "2018-02-01 13:10:45.67", + ], + str3ds_name="test_str3ds_custom_absolute", + ) + + # Add strds with relative times in seconds + gen_strds_relative( + start_times=[0, 1000, 1500], + end_times=[1000, 1500, 3600], + time_unit="seconds", + strds_name="test_strds_relative_seconds", + seed_add=200, + ) + + # Add strds with relative times in days + gen_strds_relative( + start_times=[30, 76, 78], + end_times=[76, 78, 90], + time_unit="days", + strds_name="test_strds_relative_days", + seed_add=300, + ) + yield session session.close() @@ -151,7 +184,6 @@ def gen_str3ds( semantic="mean", ) # create MapDataset objects list - str3ds_length = 10 map_dts_lst = [] for i, map_time in enumerate(str3ds_times): # Generate random map. Given seed for reproducibility @@ -187,6 +219,107 @@ def gen_str3ds( return stds +def gen_str3ds_custom_absolute( + start_times: list, end_times: list, str3ds_name: str = "test_str3ds_custom_absolute" +) -> tgis.Raster3DDataset: + """Generate a synthetic str3ds with custom absolute start and end times.""" + grass_i = GrassInterface() + + # Parse datetime strings to datetime objects + start_dts = [datetime.fromisoformat(t.replace(" ", "T")) for t in start_times] + end_dts = [datetime.fromisoformat(t.replace(" ", "T")) for t in end_times] + + # Create the str3ds + stds_id = grass_i.get_id_from_name(str3ds_name) + stds = tgis.open_new_stds( + name=stds_id, + type="str3ds", + temporaltype="absolute", + title="", + descr="", + semantic="mean", + ) + + # Create MapDataset objects list + map_dts_lst = [] + for i, (start_time, end_time) in enumerate(zip(start_dts, end_dts)): + # Generate random map. Given seed for reproducibility + formatted_date = start_time.strftime("%Y%m%d_%H%M%S") + map_name = f"test3d_abs_{formatted_date}" + gs.raster3d.mapcalc3d(exp=f"{map_name}=rand(10,100)", seed=i + 100) + + # Create MapDataset + map_id = grass_i.get_id_from_name(map_name) + map_dts = tgis.Raster3DDataset(map_id) + # Load spatial data from map + map_dts.load() + # Set time + map_dts.set_absolute_time(start_time=start_time, end_time=end_time) + # Populate the list + map_dts_lst.append(map_dts) + + # Finally register the maps + tgis.register.register_map_object_list( + type="raster3d", + map_list=map_dts_lst, + output_stds=stds, + delete_empty=False, + unit="", + ) + return stds + + +def gen_strds_relative( + start_times: list, + end_times: list, + time_unit: str, + strds_name: str, + seed_add: int = 200, +) -> tgis.RasterDataset: + """Generate a synthetic strds with relative times in seconds.""" + grass_i = GrassInterface() + + # Create the strds + stds_id = grass_i.get_id_from_name(strds_name) + stds = tgis.open_new_stds( + name=stds_id, + type="strds", + temporaltype="relative", + title="", + descr="", + semantic="mean", + ) + + # Create MapDataset objects list + map_dts_lst = [] + for i, (start_time, end_time) in enumerate(zip(start_times, end_times)): + # Generate random map + map_name = f"test_raster_sec_{start_time}" + gs.mapcalc(exp=f"{map_name}=rand(0,100)", seed=i + seed_add) + + # Create MapDataset + map_id = grass_i.get_id_from_name(map_name) + map_dts = tgis.RasterDataset(map_id) + # Load spatial data from map + map_dts.load() + # Set time + map_dts.set_relative_time( + start_time=start_time, end_time=end_time, unit=time_unit + ) + # Populate the list + map_dts_lst.append(map_dts) + + # Finally register the maps + tgis.register.register_map_object_list( + type="raster", + map_list=map_dts_lst, + output_stds=stds, + delete_empty=False, + unit=time_unit, + ) + return stds + + def create_sample_dataarray( dims_spec: dict, shape: tuple, diff --git a/tests/test_grass_interface.py b/tests/test_grass_interface.py index 87fe537..677a4d2 100644 --- a/tests/test_grass_interface.py +++ b/tests/test_grass_interface.py @@ -196,7 +196,13 @@ def test_has_mask(self): def test_list_strds(self): strds_list = GrassInterface.list_strds() - assert strds_list == [ACTUAL_STRDS] + assert ACTUAL_STRDS in strds_list + assert len(strds_list) == 3 # modis + 2 synthetic ones + + def test_list_str3ds(self): + str3ds_list = GrassInterface.list_str3ds() + assert GrassInterface.get_id_from_name(RELATIVE_STR3DS) in str3ds_list + assert len(str3ds_list) == 3 # 3 synthetic str3ds def test_get_stds_infos(self, grass_i): strds_infos = grass_i.get_stds_infos(ACTUAL_STRDS, stds_type="strds") diff --git a/tests/test_xarray_grass.py b/tests/test_xarray_grass.py index c8ca1a1..22deda0 100644 --- a/tests/test_xarray_grass.py +++ b/tests/test_xarray_grass.py @@ -149,8 +149,8 @@ def test_load_multiple_rasters(self, grass_i, temp_gisdb) -> None: ) region = grass_i.get_region() assert isinstance(test_dataset, xr.Dataset) - # z, y_3d, x_3d, y, x, absolute and relative time - assert len(test_dataset.dims) == 7 + # z, y_3d, x_3d, y, x, and one time dimension for each stds + assert len(test_dataset.dims) == 8 assert len(test_dataset) == 7 assert len(test_dataset.x_3d) == region.cols3 assert len(test_dataset.y_3d) == region.rows3 @@ -213,7 +213,10 @@ def test_load_whole_mapset(self, grass_i, temp_gisdb) -> None: + list_str3ds_name ) assert isinstance(whole_mapset, xr.Dataset) - assert len(whole_mapset.dims) == 7 + dim_num = len( + list_strds_name + list_str3ds_name + ["y", "x", "y_3d", "x_3d", "z"] + ) + assert len(whole_mapset.dims) == dim_num assert len(whole_mapset) == len(all_variables) assert all(var in whole_mapset for var in all_variables) assert len(whole_mapset.x_3d) == region.cols3 @@ -241,13 +244,16 @@ def test_drop_variables(self, grass_i, temp_gisdb) -> None: raster_3d=[ACTUAL_RASTER3D_MAP, ACTUAL_RASTER3D_MAP2], str3ds=[RELATIVE_STR3DS, ABSOLUTE_STR3DS], strds=ACTUAL_STRDS, - drop_variables=[ACTUAL_RASTER_MAP], + drop_variables=[ABSOLUTE_STR3DS, ACTUAL_RASTER_MAP], ) region = grass_i.get_region() assert isinstance(test_dataset, xr.Dataset) - # z, y_3d, x_3d, y, x, absolute and relative time - assert len(test_dataset.dims) == 7 - assert len(test_dataset) == 6 # 7 - 1 dropped variable + # Each stds has its own time dimension + num_dims = len( + ["z", "y_3d", "x_3d", "y", "x"] + [ACTUAL_STRDS, RELATIVE_STR3DS] + ) + assert len(test_dataset.dims) == num_dims + assert len(test_dataset) == 7 - 2 # dropped vars assert len(test_dataset.x_3d) == region.cols3 assert len(test_dataset.y_3d) == region.rows3 assert len(test_dataset.x) == region.cols @@ -269,7 +275,7 @@ def test_attributes_separation(self, grass_i, temp_gisdb) -> None: ) # Dataset-level attributes: only these should be present - expected_dataset_attrs = {"crs_wkt", "Conventions"} + expected_dataset_attrs = {"crs_wkt", "Conventions", "source", "history"} actual_dataset_attrs = set(test_dataset.attrs.keys()) assert actual_dataset_attrs == expected_dataset_attrs, ( f"Dataset has unexpected attributes. " @@ -284,7 +290,7 @@ def test_attributes_separation(self, grass_i, temp_gisdb) -> None: assert isinstance(test_dataset.attrs["crs_wkt"], str) # DataArray-level attributes that should NOT appear at Dataset level - dataarray_only_attrs = {"long_name", "source", "units", "comment"} + dataarray_only_attrs = {"long_name", "units", "comment"} # Check that DataArray attributes don't leak to Dataset level for attr in dataarray_only_attrs: @@ -316,8 +322,6 @@ def test_attributes_separation(self, grass_i, temp_gisdb) -> None: if "z" in test_dataset.coords: assert "axis" in test_dataset.z.attrs assert test_dataset.z.attrs["axis"] == "Z" - assert "positive" in test_dataset.z.attrs - assert test_dataset.z.attrs["positive"] == "up" assert len(test_dataset.z) == region.depths # Check time coordinate attributes (from STRDS) From 05f038a06a53228169c231385ee2708be12dc52a Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 12:55:33 -0600 Subject: [PATCH 5/9] update readme --- README.md | 52 +++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 43 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 39efb52..1a5f3f1 100644 --- a/README.md +++ b/README.md @@ -21,23 +21,27 @@ You need to install GRASS independently. >>> import xarray as xr >>> test_ds = xr.open_dataset("/home/lc/grassdata/nc_spm_08_grass7/PERMANENT/", raster=["boundary_county_500m", "elevation"]) >>> test_ds - Size: 244kB -Dimensions: (y: 150, x: 135) + Size: 253kB +Dimensions: (y: 150, x: 140) Coordinates: * y (y) float32 600B 2.2e+05 2.2e+05 ... 2.207e+05 - * x (x) float32 540B 6.383e+05 6.383e+05 ... 6.39e+05 + * x (x) float32 560B 6.383e+05 6.383e+05 ... 6.39e+05 Data variables: - boundary_county_500m (y, x) float64 162kB ... - elevation (y, x) float32 81kB ... + boundary_county_500m (y, x) float64 168kB ... + elevation (y, x) float32 84kB ... Attributes: - crs_wkt: PROJCRS["NAD83(HARN) / North Carolina",BASEGEOGCRS["NAD83(HARN... + crs_wkt: PROJCRS["NAD83(HARN) / North Carolina",BASEGEOGCRS["NAD83(H... + Conventions: CF-1.13-draft + history: 2025-10-31 18:22:16.644873+00:00: Created with xarray-grass... + source: GRASS database. project: , mapset:>> test_ds = xr.open_dataset("/home/lc/grassdata/nc_spm_08_grass7/PERMANENT/", raster=["boundary_county_500m", "elevation"], strds="LST_Day_monthly@modis_lst") +>>> test_ds + Size: 2MB +Dimensions: (y: 150, x: 140, start_time_LST_Day_monthly: 24) +Coordinates: + * y (y) float32 600B 2.2e+05 2.2e+05 ... 2.207e+05 + * x (x) float32 560B 6.383e+05 ... 6.39e+05 + * start_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... + end_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... +Data variables: + boundary_county_500m (y, x) float64 168kB ... + elevation (y, x) float32 84kB ... + LST_Day_monthly (start_time_LST_Day_monthly, y, x) int32 2MB ... +Attributes: + crs_wkt: PROJCRS["NAD83(HARN) / North Carolina",BASEGEOGCRS["NAD83(H... + Conventions: CF-1.13-draft + history: 2025-10-31 17:40:15.661101+00:00: Created with xarray-grass... + source: GRASS database. project: , mapset: Date: Fri, 31 Oct 2025 19:38:58 -0600 Subject: [PATCH 6/9] to_grass(): more flexible dims --- src/xarray_grass/to_grass.py | 409 ++++++++++++++++----------- tests/test_coord_utils.py | 4 +- tests/test_tograss.py | 16 +- tests/test_tograss_error_handling.py | 26 +- 4 files changed, 277 insertions(+), 178 deletions(-) diff --git a/src/xarray_grass/to_grass.py b/src/xarray_grass/to_grass.py index e7ea4ba..06f6742 100644 --- a/src/xarray_grass/to_grass.py +++ b/src/xarray_grass/to_grass.py @@ -15,7 +15,7 @@ import os from pathlib import Path -from typing import Mapping +from typing import Mapping, Optional from pyproj import CRS import xarray as xr @@ -28,22 +28,10 @@ from xarray_grass.coord_utils import get_region_from_xarray -# Default dimension names -default_dims = { - "start_time": "start_time", - "end_time": "end_time", - "x": "x", - "y": "y", - "x_3d": "x_3d", - "y_3d": "y_3d", - "z": "z", -} - - def to_grass( dataset: xr.Dataset | xr.DataArray, mapset: str | Path, - dims: Mapping[str, str] = None, + dims: Optional[Mapping[str, Mapping[str, str]]] = None, create: bool = False, ) -> None: """Convert an xarray.Dataset or xarray.DataArray to GRASS GIS maps. @@ -62,11 +50,12 @@ def to_grass( will be converted. mapset : str | Path Path to the target GRASS mapset. - dims : Mapping[str, str], optional + dims : Mapping[str, Mapping[str, str]], optional A mapping from standard dimension names ('start_time', 'end_time', 'x', 'y', 'x_3d', 'y_3d', 'z',) - to the actual dimension names in the dataset. For example, if your 3D dataset - east-west coordinate is named 'lon', you would pass `dims={'x_3d': 'lon'}`. + to the actual dimension names in the dataset. For example, + if your 3D variable named "pressure" has an east-west coordinate named 'lon', + you would pass `dims={'pressure': {'x_3d': 'lon'}}`. Defaults to None, which implies standard dimension names are used. create : bool, optional If True (default), the mapset will be created if it does not exist. @@ -91,6 +80,11 @@ def to_grass( project_path = mapset_path.parent gisdb = project_path.parent + if not isinstance(dataset, xr.Dataset | xr.DataArray): + raise TypeError( + f"'dataset must be either an Xarray DataArray or Dataset, not {type(dataset)}" + ) + if create: # Not until GRASS 8.5 raise NotImplementedError("'create' not yet available.") @@ -122,17 +116,14 @@ def to_grass( # Skip until grass 8.5 pass - # set the dimensions dict - if dims is not None: - if not isinstance(dims, Mapping): - raise TypeError("dims parameter must be a mapping (e.g., a dictionary).") - # Start with a copy of defaults, then update with valid user-provided dims - processed_dims = default_dims.copy() - for k, v in dims.items(): - processed_dims[k] = v - dims = processed_dims - else: - dims = default_dims.copy() + try: + input_var_names = [var_name for var_name, _ in dataset.data_vars.items()] + except AttributeError: # a dataarray + input_var_names = [dataset.name] + + # set dimensions Mapping + dim_formatter = DimensionsFormatter(input_var_names, dims) + dim_dataset = dim_formatter.get_formatted_dims() # Check if we're already in a GRASS session session = None @@ -147,161 +138,239 @@ def to_grass( else: # We're in an existing session, check if it matches the requested path gi = GrassInterface() - gisenv = gi.get_gisenv() - current_gisdb = gisenv["GISDBASE"] - current_location = gisenv["LOCATION_NAME"] - accessible_mapsets = gi.get_accessible_mapsets() + check_grass_session(gi, mapset_path) - requested_path = Path(gisdb) / Path(project_name) - current_path = Path(current_gisdb) / Path(current_location) - - if requested_path != current_path or str(mapset) not in accessible_mapsets: - raise ValueError( - f"Cannot access {mapset_path} " - f"from current GRASS session in project {current_path}. " - f"Accessible mapsets: {accessible_mapsets}." - ) try: - xarray_to_grass(dataset, gi, dims) + xarray_to_grass = XarrayToGrass(dataset, gi, dim_dataset) + xarray_to_grass.to_grass() finally: if session is not None: session.__exit__(None, None, None) -def xarray_to_grass( - dataset: xr.Dataset | xr.DataArray, - gi: GrassInterface, - dims: Mapping[str, str] = None, -) -> None: - """Convert an xarray Dataset or DataArray to GRASS maps. - This function validates the CRS and pass the individual DataArrays to the - `datarray_to_grass` function""" - grass_crs = CRS(gi.get_crs_wkt_str()) - dataset_crs = CRS(dataset.attrs["crs_wkt"]) - # TODO: reproj if not same crs - # TODO: handle no CRS for xy locations - if grass_crs != dataset_crs: - raise ValueError( - f"CRS mismatch: GRASS project CRS is {grass_crs}, " - f"but dataset CRS is {dataset_crs}." - ) - try: - for var_name, data in dataset.data_vars.items(): - datarray_to_grass(data, gi, dims) - except AttributeError: - datarray_to_grass(dataset, gi, dims) +def check_grass_session(grass_interface, mapset_path): + mapset = mapset_path.stem + project_name = mapset_path.parent.stem + project_path = mapset_path.parent + gisdb = project_path.parent + gisenv = grass_interface.get_gisenv() + current_gisdb = gisenv["GISDBASE"] + current_location = gisenv["LOCATION_NAME"] + accessible_mapsets = grass_interface.get_accessible_mapsets() -def datarray_to_grass( - data: xr.DataArray, - gi: GrassInterface, - dims: Mapping[str, str] = None, -) -> None: - """Convert an xarray DataArray to GRASS maps. + requested_path = Path(gisdb) / Path(project_name) + current_path = Path(current_gisdb) / Path(current_location) - Uses standardized (x, y) dimension naming internally. For datasets with - latitude/longitude dimensions, provide explicit mapping via dims parameter. - """ - if len(data.dims) > 4 or len(data.dims) < 2: + if requested_path != current_path or str(mapset) not in accessible_mapsets: raise ValueError( - f"Only DataArray with 2 to 4 dimensions are supported. " - f"Found {len(data.dims)} dimension(s)." + f"Cannot access {mapset_path} " + f"from current GRASS session in project {current_path}. " + f"Accessible mapsets: {accessible_mapsets}." ) - # Check for 2D spatial dimensions - is_spatial_2d = dims["x"] in data.dims and dims["y"] in data.dims - - # Check for 3D spatial dimensions - is_spatial_3d = ( - dims["x_3d"] in data.dims - and dims["y_3d"] in data.dims - and dims["z"] in data.dims - ) - - # Check for time dimension - has_time = dims["start_time"] in data.dims - - # Note: 'end_time' is also a potential temporal dimension but GRASS STRDS/STR3DS - # are typically defined by a start time. - # For simplicity 'start_time' is the primary indicator here. - - # Determine dataset type based on number of dimensions and identified dimension types - is_raster = len(data.dims) == 2 and is_spatial_2d - is_raster_3d = len(data.dims) == 3 and is_spatial_3d - is_strds = len(data.dims) == 3 and has_time and is_spatial_2d - is_str3ds = len(data.dims) == 4 and has_time and is_spatial_3d - - # Set temp region - current_region = gi.get_region() - temp_region = get_region_from_xarray(data, dims) - gi.set_region(temp_region) - # TODO: reshape to match user dims - try: - if is_raster: - gi.write_raster_map(data, data.name) - elif is_strds: - write_stds(data, gi, dims) - elif is_raster_3d: - gi.write_raster3d_map(data, data.name) - elif is_str3ds: - write_stds(data, gi, dims) - else: + +class DimensionsFormatter: + """Populate the dimension mapping based on default values and user-provided ones""" + + # Default dimension names + default_dims = { + "start_time": "start_time", + "end_time": "end_time", + "x": "x", + "y": "y", + "x_3d": "x_3d", + "y_3d": "y_3d", + "z": "z", + } + + def __init__(self, input_var_names, input_dims): + self.input_var_names = input_var_names + self.input_dims = input_dims + self._dataset_dims = {} + + # Instantiate the dimensions with default values + for var_name in input_var_names: + self._dataset_dims[var_name] = ( + self.default_dims.copy() + ) # copy avoids shared state + + if self.input_dims is not None: + self.check_input_dims() + + def check_input_dims(self): + """Check conformity of provided dims Mapping""" + if not isinstance(self.input_dims, Mapping): + raise TypeError( + "The 'dims' parameter must be of type Mapping[str, Mapping[str, str]]." + ) + for var_name, var_dims in self.input_dims.items(): + if not isinstance(var_dims, Mapping): + raise TypeError( + "The 'dims' parameter must be of type Mapping[str, Mapping[str, str]]." + ) + if var_name not in self.input_var_names: + raise ValueError( + f"Variable {var_name} not found in the input dataset. " + f"Variables found: {self.input_var_names}" + ) + + def fill_dims(self): + """Replace the default values with those given by the user.""" + for var_name, dims in self.input_dims.items(): + for dim_key, dim_value in dims.items(): + self._dataset_dims[var_name][dim_key] = dim_value + + def get_formatted_dims(self): + if self.input_dims is not None: + self.fill_dims() + return self._dataset_dims + + +class XarrayToGrass: + def __init__( + self, + dataset: xr.Dataset | xr.DataArray, + grass_interface: GrassInterface, + dims: Mapping[str, str] = None, + ): + self.dataset = dataset + self.grass_interface = grass_interface + self.dataset_dims = dims + + def to_grass(self) -> None: + """Convert an xarray Dataset or DataArray to GRASS maps. + This function validates the CRS and pass the individual DataArrays to the + `datarray_to_grass` function""" + grass_crs = CRS(self.grass_interface.get_crs_wkt_str()) + dataset_crs = CRS(self.dataset.attrs["crs_wkt"]) + # TODO: reproj if not same crs + # TODO: handle no CRS for xy locations + if grass_crs != dataset_crs: raise ValueError( - f"DataArray '{data.name}' does not match any supported GRASS dataset type. " - f"Expected 2D, 3D, STRDS, or STR3DS." + f"CRS mismatch: GRASS project CRS is {grass_crs}, " + f"but dataset CRS is {dataset_crs}." ) - finally: - # Restore the original region - gi.set_region(current_region) - - -def write_stds(data: xr.DataArray, gi: GrassInterface, dims: Mapping): - # 1. Determine the temporal coordinate and type - time_coord = data[dims["start_time"]] - time_dtype = time_coord.dtype - if isinstance(time_dtype, np.dtypes.DateTime64DType): - temporal_type = "absolute" - elif np.issubdtype(time_dtype, np.integer): - temporal_type = "relative" - time_unit = time_coord.attrs.get("units", None) - else: - raise ValueError(f"Temporal type not supported: {time_dtype}") - # 2. Determine the semantic type - # TODO: find actual type - semantic_type = "mean" - # 2.5 determine if 2D or 3D - is_3d = False - stds_type = "strds" - if len(data.isel({dims["start_time"]: 0}).dims) == 3: - is_3d = True - stds_type = "str3ds" - - # 3. Loop through the time dim: - map_list = [] - for index, time in enumerate(time_coord): - darray = data.sel({dims["start_time"]: time}) - nd_array = darray.values - # 3.1 Write each map individually - raster_name = f"{data.name}_{temporal_type}_{index}" - if not is_3d: - gi.write_raster_map(arr=nd_array, rast_name=raster_name) - else: - gi.write_raster3d_map(arr=nd_array, rast_name=raster_name) - # 3.2 populate an iterable[tuple[str, datetime | timedelta]] - time_value = time.values.item() - if temporal_type == "absolute": - absolute_time = pd.Timestamp(time_value) - map_list.append((raster_name, absolute_time.to_pydatetime())) + try: + for var_name, data in self.dataset.data_vars.items(): + self._datarray_to_grass(data, self.dataset_dims[var_name]) + except AttributeError: # DataArray + self._datarray_to_grass(self.dataset, self.dataset_dims[self.dataset.name]) + + def _datarray_to_grass( + self, + data: xr.DataArray, + dims: Mapping[str, str], + ) -> None: + """Convert an xarray DataArray to GRASS maps. + + Uses standardized (x, y) dimension naming internally. For datasets with + latitude/longitude dimensions, provide explicit mapping via dims parameter. + """ + if len(data.dims) > 4 or len(data.dims) < 2: + raise ValueError( + f"Only DataArray with 2 to 4 dimensions are supported. " + f"Found {len(data.dims)} dimension(s)." + ) + + # Check for 2D spatial dimensions + is_spatial_2d = dims["x"] in data.dims and dims["y"] in data.dims + + # Check for 3D spatial dimensions + is_spatial_3d = ( + dims["x_3d"] in data.dims + and dims["y_3d"] in data.dims + and dims["z"] in data.dims + ) + + # Check for time dimension + has_time = dims["start_time"] in data.dims + + # Note: 'end_time' is also a potential temporal dimension but GRASS STRDS/STR3DS + # are typically defined by a start time. + # For simplicity 'start_time' is the primary indicator here. + + # Determine dataset type based on number of dimensions and identified dimension types + is_raster = len(data.dims) == 2 and is_spatial_2d + is_raster_3d = len(data.dims) == 3 and is_spatial_3d + is_strds = len(data.dims) == 3 and has_time and is_spatial_2d + is_str3ds = len(data.dims) == 4 and has_time and is_spatial_3d + + # Set temp region + current_region = self.grass_interface.get_region() + temp_region = get_region_from_xarray(data, dims) + print(temp_region) + self.grass_interface.set_region(temp_region) + # TODO: reshape to match user dims + try: + if is_raster: + self.grass_interface.write_raster_map(data, data.name) + elif is_strds: + self._write_stds(data, dims) + elif is_raster_3d: + self.grass_interface.write_raster3d_map(data, data.name) + elif is_str3ds: + self._write_stds(data, dims) + else: + raise ValueError( + f"DataArray '{data.name}' does not match any supported GRASS dataset type. " + f"Expected 2D, 3D, STRDS, or STR3DS." + ) + finally: + # Restore the original region + self.grass_interface.set_region(current_region) + + def _write_stds(self, data: xr.DataArray, dims: Mapping): + # 1. Determine the temporal coordinate and type + time_coord = data[dims["start_time"]] + time_dtype = time_coord.dtype + if isinstance(time_dtype, np.dtypes.DateTime64DType): + temporal_type = "absolute" + elif np.issubdtype(time_dtype, np.integer): + temporal_type = "relative" + time_unit = time_coord.attrs.get("units", None) else: - relative_time = pd.Timedelta(time_value, unit=time_unit) - map_list.append((raster_name, relative_time.to_pytimedelta())) - # 4. Create STDS and register the maps in it - gi.register_maps_in_stds( - stds_title="", - stds_name=data.name, - stds_desc="", - map_list=map_list, - semantic=semantic_type, - t_type=temporal_type, - stds_type=stds_type, - ) + raise ValueError(f"Temporal type not supported: {time_dtype}") + # 2. Determine the semantic type + # TODO: find actual type + semantic_type = "mean" + # 2.5 determine if 2D or 3D + is_3d = False + stds_type = "strds" + if len(data.isel({dims["start_time"]: 0}).dims) == 3: + is_3d = True + stds_type = "str3ds" + + # 3. Loop through the time dim: + map_list = [] + for index, time in enumerate(time_coord): + darray = data.sel({dims["start_time"]: time}) + nd_array = darray.values + # 3.1 Write each map individually + raster_name = f"{data.name}_{temporal_type}_{index}" + if not is_3d: + self.grass_interface.write_raster_map( + arr=nd_array, rast_name=raster_name + ) + else: + self.grass_interface.write_raster3d_map( + arr=nd_array, rast_name=raster_name + ) + # 3.2 populate an iterable[tuple[str, datetime | timedelta]] + time_value = time.values.item() + if temporal_type == "absolute": + absolute_time = pd.Timestamp(time_value) + map_list.append((raster_name, absolute_time.to_pydatetime())) + else: + relative_time = pd.Timedelta(time_value, unit=time_unit) + map_list.append((raster_name, relative_time.to_pytimedelta())) + # 4. Create STDS and register the maps in it + self.grass_interface.register_maps_in_stds( + stds_title="", + stds_name=data.name, + stds_desc="", + map_list=map_list, + semantic=semantic_type, + t_type=temporal_type, + stds_type=stds_type, + ) diff --git a/tests/test_coord_utils.py b/tests/test_coord_utils.py index 8d2e6c2..2208b77 100644 --- a/tests/test_coord_utils.py +++ b/tests/test_coord_utils.py @@ -17,9 +17,11 @@ import numpy as np import pytest -from xarray_grass.to_grass import default_dims +from xarray_grass.to_grass import DimensionsFormatter from xarray_grass.coord_utils import get_region_from_xarray +default_dims = DimensionsFormatter.default_dims + def test_get_region_from_xarray_2d_xy(): """Test get_region_from_xarray with a 2D DataArray using x, y.""" diff --git a/tests/test_tograss.py b/tests/test_tograss.py index dfe8660..c3397b6 100644 --- a/tests/test_tograss.py +++ b/tests/test_tograss.py @@ -400,7 +400,7 @@ def test_dataarray_to_str3ds_conversion( dataset=sample_da, mapset=mapset_arg, create=False, - dims={"start_time": "time"}, + dims={"test_str3ds_vol": {"start_time": "time"}}, ) try: available_str3ds = grass_i.list_str3ds() @@ -529,7 +529,10 @@ def test_dataset_conversion_mixed_types( dataset=sample_ds, mapset=mapset_arg, create=False, - dims={"start_time": "time"}, + dims={ + "strds_var": {"start_time": "time"}, + "str3ds_var": {"start_time": "time"}, + }, ) try: grass_objects = grass_i.list_grass_objects() @@ -728,7 +731,12 @@ def test_mapset_creation_false_existing_mapset( ) @pytest.mark.parametrize( - "rename_map", [None, {"y": "northing", "x": "easting"}, {"y": "custom_y"}] + "rename_map", + [ + None, + {"dims_test_raster": {"y": "northing", "x": "easting"}}, + {"dims_test_raster": {"y": "custom_y"}}, + ], ) def test_dims_mapping( self, @@ -756,7 +764,7 @@ def test_dims_mapping( ) if rename_map: - sample_da = sample_da.rename(rename_map) + sample_da = sample_da.rename(rename_map[da_name]) to_grass( dataset=sample_da, mapset=str(mapset_path), diff --git a/tests/test_tograss_error_handling.py b/tests/test_tograss_error_handling.py index ac4e535..9538eed 100644 --- a/tests/test_tograss_error_handling.py +++ b/tests/test_tograss_error_handling.py @@ -226,8 +226,8 @@ def test_invalid_dataset_type(self, temp_gisdb, grass_i: GrassInterface): invalid_datasets = [123, "a string", [1, 2, 3], {"data": np.array([1])}, None] for invalid_ds in invalid_datasets: with pytest.raises( - AttributeError, - match=r"object has no attribute 'attrs'", + TypeError, + match=r"'dataset must be either an Xarray DataArray or Dataset", ): to_grass(dataset=invalid_ds, mapset=str(mapset_path), create=False) @@ -253,7 +253,7 @@ def test_invalid_dims_parameter_type(self, temp_gisdb, grass_i: GrassInterface): for invalid_dims in invalid_dims_params: with pytest.raises( TypeError, - match=r"dims parameter must be a mapping", # More specific regex + match=r"'dims' parameter must be of type", ): to_grass( dataset=sample_da, @@ -284,3 +284,23 @@ def test_invalid_mapset_parameter_type(self, temp_gisdb, grass_i: GrassInterface match=r"(mapset parameter must be a string or a Path|Invalid mapset type|argument should be a str or an os.PathLike object)", ): to_grass(dataset=sample_da, mapset=invalid_mapset, create=True) + + def test_invalid_dims_invalid_var(self, temp_gisdb, grass_i: GrassInterface): + session_crs_wkt = grass_i.get_crs_wkt_str() + mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / temp_gisdb.mapset + sample_da = create_sample_dataarray( + dims_spec={"y": np.arange(2.0), "x": np.arange(2.0)}, + shape=(2, 2), + crs_wkt=session_crs_wkt, + name="data_for_dims_validation", + ) + with pytest.raises( + ValueError, + match=r"not found in the input dataset. Variables found", + ): + to_grass( + dataset=sample_da, + mapset=str(mapset_path), + dims={"invalid_var_name": {"x": "my_x"}}, + create=False, + ) From a7ba14e45e525bfcfbefc3487006e2be37f45da6 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 19:54:44 -0600 Subject: [PATCH 7/9] remove print statement --- src/xarray_grass/to_grass.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/xarray_grass/to_grass.py b/src/xarray_grass/to_grass.py index 06f6742..dd012e9 100644 --- a/src/xarray_grass/to_grass.py +++ b/src/xarray_grass/to_grass.py @@ -299,9 +299,8 @@ def _datarray_to_grass( # Set temp region current_region = self.grass_interface.get_region() temp_region = get_region_from_xarray(data, dims) - print(temp_region) self.grass_interface.set_region(temp_region) - # TODO: reshape to match user dims + # TODO: reshape to match userGRASS expected dims order try: if is_raster: self.grass_interface.write_raster_map(data, data.name) From 4e9543e633c1892e943f182950b3b863a2cf94a9 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 20:18:34 -0600 Subject: [PATCH 8/9] add overwrite switch to to_grass() --- src/xarray_grass/to_grass.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/xarray_grass/to_grass.py b/src/xarray_grass/to_grass.py index dd012e9..8813150 100644 --- a/src/xarray_grass/to_grass.py +++ b/src/xarray_grass/to_grass.py @@ -32,6 +32,7 @@ def to_grass( dataset: xr.Dataset | xr.DataArray, mapset: str | Path, dims: Optional[Mapping[str, Mapping[str, str]]] = None, + overwrite: bool = False, create: bool = False, ) -> None: """Convert an xarray.Dataset or xarray.DataArray to GRASS GIS maps. @@ -133,11 +134,11 @@ def to_grass( gisdb=str(gisdb), location=str(project_name), mapset=str(mapset) ) session.__enter__() - gi = GrassInterface() + gi = GrassInterface(overwrite) else: # We're in an existing session, check if it matches the requested path - gi = GrassInterface() + gi = GrassInterface(overwrite) check_grass_session(gi, mapset_path) try: From 5930daa67c81b97778ca07ffbed9ec24d67bde06 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Fri, 31 Oct 2025 20:19:02 -0600 Subject: [PATCH 9/9] add write to grass to readme --- README.md | 78 +++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 53 insertions(+), 25 deletions(-) diff --git a/README.md b/README.md index 1a5f3f1..67338df 100644 --- a/README.md +++ b/README.md @@ -55,29 +55,7 @@ You can list the accessible mapsets with `g.mapsets` from GRASS. In GRASS, the time dimension of various STDSs is not homogeneous, as it is for the spatial coordinates. To reflect this, xarray-grass will create one time dimension for each STDS loaded. -From within a grass session, it is possible to access various mapsets: - -```python ->>> test_ds = xr.open_dataset("/home/lc/grassdata/nc_spm_08_grass7/PERMANENT/", raster=["boundary_county_500m", "elevation"], strds="LST_Day_monthly@modis_lst") ->>> test_ds - Size: 2MB -Dimensions: (y: 150, x: 140, start_time_LST_Day_monthly: 24) -Coordinates: - * y (y) float32 600B 2.2e+05 2.2e+05 ... 2.207e+05 - * x (x) float32 560B 6.383e+05 ... 6.39e+05 - * start_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... - end_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... -Data variables: - boundary_county_500m (y, x) float64 168kB ... - elevation (y, x) float32 84kB ... - LST_Day_monthly (start_time_LST_Day_monthly, y, x) int32 2MB ... -Attributes: - crs_wkt: PROJCRS["NAD83(HARN) / North Carolina",BASEGEOGCRS["NAD83(H... - Conventions: CF-1.13-draft - history: 2025-10-31 17:40:15.661101+00:00: Created with xarray-grass... - source: GRASS database. project: , mapset: Size: 3MB +Dimensions: (y: 165, x: 179, start_time_LST_Day_monthly: 24) +Coordinates: + * y (y) float32 660B 2.196e+05 ... 2.212e+05 + * x (x) float32 716B 6.377e+05 ... 6.395e+05 + * start_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... + end_time_LST_Day_monthly (start_time_LST_Day_monthly) datetime64[ns] 192B ... +Data variables: + boundary_county_500m (y, x) float64 236kB ... + elevation (y, x) float32 118kB ... + LST_Day_monthly (start_time_LST_Day_monthly, y, x) int32 3MB ... +Attributes: + crs_wkt: PROJCRS["NAD83(HARN) / North Carolina",BASEGEOGCRS["NAD83(H... + Conventions: CF-1.13-draft + history: 2025-11-01 02:10:24.652838+00:00: Created with xarray-grass... + source: GRASS database. project: , mapset: