diff --git a/README.md b/README.md index 347c098..fcd04aa 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,7 @@ A [GRASS](https://grass.osgeo.org/) backend for [Xarray](https://xarray.dev/). Explore all your GRASS rasters with Xarray. +Import zarr or NetCDF into your GRASS database. ## Installation @@ -18,29 +19,56 @@ You need to install GRASS independently. ## Loading GRASS data as an Xarray Dataset ```python ->>> 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: 253kB -Dimensions: (y: 150, x: 140) +import xarray as xr +import grass_session + +import grass.script as gscript + + +with grass_session.Session( + gisdb="/home/laurent/data/grassdata/", + location="nc_spm_08_grass7_test", + mapset="PERMANENT", +): + # Make modis_lst mapset accessible + gscript.run_command("g.mapsets", mapset="modis_lst", operation="add") + test_ds = xr.open_dataset( + "", # No need to pass a path, the information is taken from the active grass session + engine="xarray_grass", # If no path is given, then the engine must be specified + raster=["boundary_county_500m", "elevation"], + strds="LST_Day_monthly@modis_lst", + ) + + print(test_ds) +``` + +``` +Search path not modified + Size: 6MB +Dimensions: (y: 165, x: 179, 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.383e+05 ... 6.39e+05 + * 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 168kB ... - elevation (y, x) float32 84kB ... + boundary_county_500m (y, x) float64 236kB ... + elevation (y, x) float32 118kB ... + LST_Day_monthly (start_time_LST_Day_monthly, y, x) int64 6MB ... Attributes: 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:=0.5", "numpy>=2.2.5", "pyproj>=3.7.1", "pandas>=2.2.3", @@ -28,6 +27,7 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", "Topic :: Scientific/Engineering", ] @@ -39,6 +39,7 @@ Issues = "https://github.com/lrntct/xarray-grass/issues" [dependency-groups] dev = [ + "grass-session>=0.5", "pre-commit>=4.2.0", "pytest>=8.3.5", "pytest-random-order>=1.1.1", diff --git a/src/xarray_grass/grass_interface.py b/src/xarray_grass/grass_interface.py index be57a94..f66f8d3 100644 --- a/src/xarray_grass/grass_interface.py +++ b/src/xarray_grass/grass_interface.py @@ -25,7 +25,6 @@ # Needed to import grass modules -import grass_session # noqa: F401 import grass.script as gs from grass.script import array as garray import grass.pygrass.utils as gutils @@ -98,6 +97,10 @@ def __init__(self, overwrite: bool = False): if "GISRC" not in os.environ: raise RuntimeError("GRASS session not set.") self.overwrite = overwrite + if self.overwrite: + os.environ["GRASS_OVERWRITE"] = "1" + else: + os.environ["GRASS_OVERWRITE"] = "0" tgis.init() @staticmethod @@ -505,3 +508,33 @@ def register_maps_in_stds( unit=t_unit, ) return self + + def get_coordinates(self, raster_3d: bool) -> dict[str : np.ndarray]: + """return np.ndarray of coordinates from the GRASS region.""" + current_region = self.get_region() + lim_e = current_region.e + lim_w = current_region.w + lim_n = current_region.n + lim_s = current_region.s + lim_t = current_region.t + lim_b = current_region.b + dz = current_region.tbres + if raster_3d: + dx = current_region.ewres3 + dy = current_region.nsres3 + else: + dx = current_region.ewres + dy = current_region.nsres + # GRASS limits are at the edge of the region. + # In the exported arrays, coordinates are at the center of the cell + # Stop not changed to include it in the range + start_w = lim_w + dx / 2 + stop_e = lim_e + start_s = lim_s + dy / 2 + stop_n = lim_n + start_b = lim_b + dz / 2 + stop_t = lim_t + x_coords = np.arange(start=start_w, stop=stop_e, step=dx, dtype=np.float32) + y_coords = np.arange(start=start_s, stop=stop_n, step=dy, dtype=np.float32) + z_coords = np.arange(start=start_b, stop=stop_t, step=dz, dtype=np.float32) + return {"x": x_coords, "y": y_coords, "z": z_coords} diff --git a/src/xarray_grass/to_grass.py b/src/xarray_grass/to_grass.py index 4a2b39a..96e2c96 100644 --- a/src/xarray_grass/to_grass.py +++ b/src/xarray_grass/to_grass.py @@ -13,35 +13,30 @@ GNU General Public License for more details. """ +from __future__ import annotations import os -from pathlib import Path -from typing import Mapping, Optional +from typing import TYPE_CHECKING, Mapping, Optional from pyproj import CRS import xarray as xr import numpy as np import pandas as pd -import grass_session # noqa: F401 -from xarray_grass.grass_interface import GrassInterface -from xarray_grass.xarray_grass import dir_is_grass_mapset, dir_is_grass_project from xarray_grass.coord_utils import get_region_from_xarray +if TYPE_CHECKING: + from xarray_grass.grass_interface import GrassInterface + 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. - This function handles the setup of the GRASS environment and session - management. It can create a new mapset if specified and not already - existing. It then calls the appropriate internal functions to perform - the conversion of the xarray object's data variables into GRASS raster, - raster 3D, STRDS, or STR3DS object. + Perform the conversion of the xarray object's data variables into GRASS + raster, raster 3D, STRDS, or STR3DS object. Parameters @@ -49,8 +44,6 @@ def to_grass( dataset : xr.Dataset | xr.DataArray The xarray object to convert. If a Dataset, each data variable will be converted. - mapset : str | Path - Path to the target GRASS mapset. dims : Mapping[str, Mapping[str, str]], optional A mapping from standard dimension names ('start_time', 'end_time', 'x', 'y', 'x_3d', 'y_3d', 'z',) @@ -66,109 +59,32 @@ def to_grass( Returns ------- None - - Raises - ------ - ValueError - If the provided `mapset` path is invalid, not a GRASS mapset, - or if its parent directory is not a valid GRASS project when - `create` is False or the mapset doesn't exist. - If the target mapset is not accessible from an existing GRASS session. """ - mapset_path = Path(mapset) - mapset = mapset_path.stem - project_name = mapset_path.parent.stem - 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.") - - if mapset_path.is_file(): - raise ValueError(f"Mapset path '{mapset_path}' is a file, not a directory.") - - # Prioritize check for create=False with a non-existent mapset - if not mapset_path.is_dir() and not create: - raise ValueError( - f"Mapset path '{mapset_path}' is not a valid directory and create is False." - ) - - if mapset_path.is_dir() and not dir_is_grass_mapset(mapset_path): - raise ValueError( - f"Path '{mapset_path}' exists but is not a valid GRASS mapset." - ) - - if not mapset_path.is_dir() and create and not dir_is_grass_project(project_path): - raise ValueError( - f"Mapset '{mapset_path}' not found and its parent directory " - f"'{project_path}' is not a valid GRASS project." + if "GISRC" not in os.environ: + raise RuntimeError( + "GRASS session not detected. " + "Please setup a GRASS session before trying to access GRASS data." ) - if not mapset_path.is_dir() and dir_is_grass_project(project_path) and create: - # gs.run_command( - # "g.mapset", mapset=mapset_path.name, project=project_name, flags="c" - # ) - # Skip until grass 8.5 - pass + from xarray_grass.grass_interface import GrassInterface - try: + if isinstance(dataset, xr.Dataset): input_var_names = [var_name for var_name, _ in dataset.data_vars.items()] - except AttributeError: # a dataarray + elif isinstance(dataset, xr.DataArray): input_var_names = [dataset.name] + else: + raise TypeError( + f"'dataset must be either an Xarray DataArray or Dataset, not {type(dataset)}" + ) # 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 - if "GISRC" not in os.environ: - # No existing session, create a new one - session = grass_session.Session( - gisdb=str(gisdb), location=str(project_name), mapset=str(mapset) - ) - session.__enter__() - gi = GrassInterface(overwrite) - - else: - # We're in an existing session, check if it matches the requested path - gi = GrassInterface(overwrite) - check_grass_session(gi, mapset_path) - - try: - 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 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() - - 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}." - ) + # Write to grass + gi = GrassInterface(overwrite) + xarray_to_grass = XarrayToGrass(dataset, gi, dim_dataset) + xarray_to_grass.to_grass() class DimensionsFormatter: @@ -367,6 +283,26 @@ def _write_stds(self, data: xr.DataArray, dims: Mapping): stds_type = "str3ds" arr_type = "raster3d" + # Check if exists + if "strds" == stds_type: + if ( + not self.grass_interface.overwrite + and self.grass_interface.name_is_strds(data.name) + ): + raise RuntimeError( + f"STRDS {data.name} already exists and will not be overwritten." + ) + elif "str3ds" == stds_type: + if ( + not self.grass_interface.overwrite + and self.grass_interface.name_is_str3ds(data.name) + ): + raise RuntimeError( + f"STR3DS {data.name} already exists and will not be overwritten." + ) + else: + raise ValueError(f"Unknown STDS type '{stds_type}'.") + # 3. Loop through the time dim: map_list = [] for index, time in enumerate(time_coord): diff --git a/src/xarray_grass/xarray_grass.py b/src/xarray_grass/xarray_grass.py index 9dbd7f5..44c500f 100644 --- a/src/xarray_grass/xarray_grass.py +++ b/src/xarray_grass/xarray_grass.py @@ -17,10 +17,8 @@ from pathlib import Path from typing import Iterable, Optional -import numpy as np from xarray.backends import BackendEntrypoint import xarray as xr -import grass_session # noqa: F401 import xarray_grass from xarray_grass.grass_interface import GrassInterface @@ -53,8 +51,23 @@ def open_dataset( drop_variables: Iterable[str], ) -> xr.Dataset: """Open GRASS project or mapset as an xarray.Dataset. + Requires an active GRASS session. TODO: add support for whole project. """ + if "GISRC" not in os.environ: + raise RuntimeError( + "GRASS session not detected. " + "Please setup a GRASS session before trying to access GRASS data." + ) + + if filename_or_obj: + dirpath = Path(filename_or_obj) + if not dir_is_grass_mapset(dirpath): + raise ValueError(f"{filename_or_obj} is not a GRASS mapset") + self.check_accessible_mapset(filename_or_obj) + + self.grass_interface = GrassInterface() + open_func_params = dict( raster_list=raster, raster_3d_list=raster_3d, @@ -62,38 +75,7 @@ def open_dataset( str3ds_list=str3ds, ) if not any([raster, raster_3d, strds, str3ds]): - # Load the whole mapset. - # If a map is in a STDS, do not load it as a single map. - gi = GrassInterface() - grass_objects = gi.list_grass_objects() - # strds - rasters_in_strds = [] - for strds_name in grass_objects["strds"]: - maps_in_strds = gi.list_maps_in_strds(strds_name) - rasters_in_strds.extend([map_data.id for map_data in maps_in_strds]) - if open_func_params["strds_list"] is None: - open_func_params["strds_list"] = [strds_name] - else: - open_func_params["strds_list"].append(strds_name) - raster3ds_in_str3ds = [] - # str3ds - for str3ds_name in grass_objects["str3ds"]: - maps_in_str3ds = gi.list_maps_in_str3ds(str3ds_name) - raster3ds_in_str3ds.extend([map_data.id for map_data in maps_in_str3ds]) - if open_func_params["str3ds_list"] is None: - open_func_params["str3ds_list"] = [str3ds_name] - else: - open_func_params["str3ds_list"].append(str3ds_name) - # rasters not in strds - open_func_params["raster_list"] = [ - name for name in grass_objects["raster"] if name not in rasters_in_strds - ] - # rasters 3d not in str3ds - open_func_params["raster_3d_list"] = [ - name - for name in grass_objects["raster_3d"] - if name not in raster3ds_in_str3ds - ] + self._list_all_mapset(open_func_params) else: # Format str inputs into list for object_type, elem in open_func_params.items(): @@ -110,108 +92,58 @@ def open_dataset( name for name in grass_obj_name_list if name not in drop_variables ] - return open_grass_maps(filename_or_obj, **open_func_params) + return self._open_grass_maps(filename_or_obj, **open_func_params) def guess_can_open(self, filename_or_obj) -> bool: """infer if the path is a GRASS mapset. TODO: add support for whole project.""" return dir_is_grass_mapset(filename_or_obj) + def _list_all_mapset(self, open_func_params): + """List map objects in the whole mapset. + If a map is part of a STDS, do not list it as a single map. + """ + grass_objects = self.grass_interface.list_grass_objects() + # strds + rasters_in_strds = [] + for strds_name in grass_objects["strds"]: + maps_in_strds = self.grass_interface.list_maps_in_strds(strds_name) + rasters_in_strds.extend([map_data.id for map_data in maps_in_strds]) + if open_func_params["strds_list"] is None: + open_func_params["strds_list"] = [strds_name] + else: + open_func_params["strds_list"].append(strds_name) + raster3ds_in_str3ds = [] + # str3ds + for str3ds_name in grass_objects["str3ds"]: + maps_in_str3ds = self.grass_interface.list_maps_in_str3ds(str3ds_name) + raster3ds_in_str3ds.extend([map_data.id for map_data in maps_in_str3ds]) + if open_func_params["str3ds_list"] is None: + open_func_params["str3ds_list"] = [str3ds_name] + else: + open_func_params["str3ds_list"].append(str3ds_name) + # rasters not in strds + open_func_params["raster_list"] = [ + name for name in grass_objects["raster"] if name not in rasters_in_strds + ] + # rasters 3d not in str3ds + open_func_params["raster_3d_list"] = [ + name + for name in grass_objects["raster_3d"] + if name not in raster3ds_in_str3ds + ] -def dir_is_grass_mapset(filename_or_obj: str | Path) -> bool: - """ - Check if the given path is a GRASS mapset. - """ - try: - dirpath = Path(filename_or_obj) - except TypeError: - return False - if dirpath.is_dir(): - wind_file = dirpath / Path("WIND") - # A newly created mapset might only have WIND, VAR appears later. - if wind_file.exists(): - return True - return False - - -def dir_is_grass_project(filename_or_obj: str | Path) -> bool: - """Return True if a subdir named PERMANENT is present.""" - try: + def check_accessible_mapset(self, filename_or_obj): dirpath = Path(filename_or_obj) - except TypeError: - return False - if dirpath.is_dir(): - return (dirpath / Path("PERMANENT")).is_dir() - else: - return False - - -def get_coordinates(grass_i: GrassInterface, raster_3d: bool) -> dict: - """return xarray coordinates from GRASS region.""" - current_region = grass_i.get_region() - lim_e = current_region.e - lim_w = current_region.w - lim_n = current_region.n - lim_s = current_region.s - lim_t = current_region.t - lim_b = current_region.b - dz = current_region.tbres - if raster_3d: - dx = current_region.ewres3 - dy = current_region.nsres3 - else: - dx = current_region.ewres - dy = current_region.nsres - # GRASS limits are at the edge of the region. - # In the exported DataArray, coordinates are at the center of the cell - # Stop not changed to include it in the range - start_w = lim_w + dx / 2 - stop_e = lim_e - start_s = lim_s + dy / 2 - stop_n = lim_n - start_b = lim_b + dz / 2 - stop_t = lim_t - x_coords = np.arange(start=start_w, stop=stop_e, step=dx, dtype=np.float32) - y_coords = np.arange(start=start_s, stop=stop_n, step=dy, dtype=np.float32) - z_coords = np.arange(start=start_b, stop=stop_t, step=dz, dtype=np.float32) - return {"x": x_coords, "y": y_coords, "z": z_coords} - - -def open_grass_maps( - filename_or_obj: str | Path, - raster_list: Iterable[str] = None, - raster_3d_list: Iterable[str] = None, - strds_list: Iterable[str] = None, - str3ds_list: Iterable[str] = None, - raise_on_not_found: bool = True, -) -> xr.Dataset: - """ - Open a GRASS mapset and return an xarray dataset. - """ - dirpath = Path(filename_or_obj) - if not dir_is_grass_mapset(dirpath): - raise ValueError(f"{filename_or_obj} is not a GRASS mapset") - mapset = dirpath.stem - project = dirpath.parent.stem - gisdb = dirpath.parent.parent - - # Check if we're already in a GRASS session - session = None - if "GISRC" not in os.environ: - # No existing session, create a new one - session = grass_session.Session( - gisdb=str(gisdb), location=str(project), mapset=str(mapset) - ) - session.__enter__() - gi = GrassInterface() - - else: - # We're in an existing session, check if it matches the requested path - gi = GrassInterface() - gisenv = gi.get_gisenv() + mapset = dirpath.stem + project_path = dirpath.parent + gisdb_path = project_path.parent + project = project_path.stem + gisdb = gisdb_path + gisenv = self.grass_interface.get_gisenv() current_gisdb = gisenv["GISDBASE"] current_location = gisenv["LOCATION_NAME"] - accessible_mapsets = gi.get_accessible_mapsets() + accessible_mapsets = self.grass_interface.get_accessible_mapsets() requested_path = Path(gisdb) / Path(project) current_path = Path(current_gisdb) / Path(current_location) @@ -223,31 +155,42 @@ def open_grass_maps( f"{current_gisdb}/{current_location}. " f"Accessible mapsets: {accessible_mapsets}." ) - try: + + def _open_grass_maps( + self, + filename_or_obj: str | Path, + raster_list: Iterable[str] = None, + raster_3d_list: Iterable[str] = None, + strds_list: Iterable[str] = None, + str3ds_list: Iterable[str] = None, + ) -> xr.Dataset: + """ + Open a GRASS mapset and return an xarray dataset. + """ # Configuration for processing different GRASS map types map_processing_configs = [ { "input_list": raster_list, - "existence_check_method": gi.name_is_raster, - "open_function": open_grass_raster, + "existence_check_method": self.grass_interface.name_is_raster, + "open_function": self._open_grass_raster, "not_found_key": "raster", }, { "input_list": raster_3d_list, - "existence_check_method": gi.name_is_raster_3d, - "open_function": open_grass_raster_3d, + "existence_check_method": self.grass_interface.name_is_raster_3d, + "open_function": self._open_grass_raster_3d, "not_found_key": "raster_3d", }, { "input_list": strds_list, - "existence_check_method": gi.name_is_strds, - "open_function": open_grass_strds, + "existence_check_method": self.grass_interface.name_is_strds, + "open_function": self._open_grass_strds, "not_found_key": "strds", }, { "input_list": str3ds_list, - "existence_check_method": gi.name_is_str3ds, - "open_function": open_grass_str3ds, + "existence_check_method": self.grass_interface.name_is_str3ds, + "open_function": self._open_grass_str3ds, "not_found_key": "str3ds", }, ] @@ -260,271 +203,309 @@ def open_grass_maps( 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) + data_array = config["open_function"](map_name) raw_coords_list.append(data_array.coords) data_array_list.append(data_array) - if raise_on_not_found and any(not_found.values()): + if any(not_found.values()): raise ValueError(f"Objects not found: {not_found}") - crs_wkt = gi.get_crs_wkt_str() - finally: - if session is not None: - session.__exit__(None, None, None) - - coords_dict = {} - for coords in 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": crs_wkt, - "Conventions": "CF-1.13-draft", - # "title": "", - "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": "", - } - dataset = xr.Dataset(data_vars=data_array_dict, coords=coords_dict, attrs=attrs) - return dataset - - -def open_grass_raster(raster_name: str, grass_i: GrassInterface) -> xr.DataArray: - """Open a single raster map.""" - x_coords, y_coords, _ = get_coordinates(grass_i, raster_3d=False).values() - dims = ["y", "x"] - coordinates = dict.fromkeys(dims) - coordinates["x"] = x_coords - coordinates["y"] = y_coords - raster_array = grass_i.read_raster_map(raster_name) - data_array = xr.DataArray( - raster_array, - coords=coordinates, - dims=dims, - name=grass_i.get_name_from_id(raster_name), - ) - # Add CF attributes - r_infos = grass_i.get_raster_info(raster_name) - da_with_attrs = set_cf_coordinates(data_array, gi=grass_i, is_3d=False) - da_with_attrs.attrs["long_name"] = r_infos.get("title", "") - da_with_attrs.attrs["source"] = ",".join([r_infos["source1"], r_infos["source2"]]) - da_with_attrs.attrs["units"] = r_infos.get("units", "") - da_with_attrs.attrs["comment"] = r_infos.get("comments", "") - # CF attributes "institution" and "references" - # Do not correspond to a direct GRASS value. - return da_with_attrs - - -def open_grass_raster_3d(raster_3d_name: str, grass_i: GrassInterface) -> xr.DataArray: - """Open a single 3D raster map.""" - x_coords, y_coords, z_coords = get_coordinates(grass_i, raster_3d=True).values() - dims = ["z", "y_3d", "x_3d"] - coordinates = dict.fromkeys(dims) - coordinates["x_3d"] = x_coords - coordinates["y_3d"] = y_coords - coordinates["z"] = z_coords - raster_array = grass_i.read_raster3d_map(raster_3d_name) - - data_array = xr.DataArray( - raster_array, - coords=coordinates, - dims=dims, - name=grass_i.get_name_from_id(raster_3d_name), - ) - # Add CF attributes - r3_infos = grass_i.get_raster3d_info(raster_3d_name) - da_with_attrs = set_cf_coordinates( - data_array, gi=grass_i, is_3d=True, z_unit=r3_infos["vertical_units"] - ) - da_with_attrs.attrs["long_name"] = r3_infos.get("title", "") - da_with_attrs.attrs["source"] = ",".join([r3_infos["source1"], r3_infos["source2"]]) - da_with_attrs.attrs["units"] = r3_infos.get("units", "") - da_with_attrs.attrs["comment"] = r3_infos.get("comments", "") - # CF attributes "institution" and "references" - # Do not correspond to a direct GRASS value. - return da_with_attrs - - -def open_grass_strds(strds_name: str, grass_i: GrassInterface) -> xr.DataArray: - """Open a STRDS with lazy loading - data is only loaded when accessed""" - 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_id, stds_type="strds") - if strds_infos.temporal_type == "absolute": - time_unit = None - else: - time_unit = strds_infos.time_unit - start_time_dim = f"start_time_{strds_name}" - end_time_dim = f"end_time_{strds_name}" - - map_list = grass_i.list_maps_in_strds(strds_id) - region = grass_i.get_region() - - # Create a single backend array for the entire STRDS - backend_array = GrassSTDSBackendArray( - shape=(len(map_list), region.rows, region.cols), - dtype=map_list[0].dtype, - map_list=map_list, - map_type="raster", - grass_interface=grass_i, - ) - lazy_array = xr.core.indexing.LazilyIndexedArray(backend_array) - - # Create Variable with lazy array - var = xr.Variable(dims=[start_time_dim, "y", "x"], data=lazy_array) - - # Extract time coordinates - start_times = [map_data.start_time for map_data in map_list] - end_times = [map_data.end_time for map_data in map_list] - - # Create coordinates - coordinates = { - "x": x_coords, - "y": y_coords, - start_time_dim: start_times, - end_time_dim: (start_time_dim, end_times), - } - - # Convert to DataArray - data_array = xr.DataArray( - var, - coords=coordinates, - name=strds_name, - ) - - # Add CF attributes - r_infos = grass_i.get_raster_info(map_list[0].id) - da_with_attrs = set_cf_coordinates( - data_array, - 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"]]) - da_with_attrs.attrs["units"] = r_infos.get("units", "") - da_with_attrs.attrs["comment"] = r_infos.get("comments", "") - # CF attributes "institution" and "references" - # Do not correspond to a direct GRASS value. - return da_with_attrs - - -def open_grass_str3ds(str3ds_name: str, grass_i: GrassInterface) -> xr.DataArray: - """Open a STR3DS with lazy loading - data is only loaded when accessed - 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_id, stds_type="str3ds") - if strds_infos.temporal_type == "absolute": - time_unit = None - else: - time_unit = strds_infos.time_unit - start_time_dim = f"start_time_{str3ds_name}" - end_time_dim = f"end_time_{str3ds_name}" - - map_list = grass_i.list_maps_in_str3ds(str3ds_id) - region = grass_i.get_region() - - # Create a single backend array for the entire STR3DS - backend_array = GrassSTDSBackendArray( - shape=(len(map_list), region.depths, region.rows3, region.cols3), - dtype=map_list[0].dtype, - map_list=map_list, - map_type="raster3d", - grass_interface=grass_i, - ) - lazy_array = xr.core.indexing.LazilyIndexedArray(backend_array) - - # Create Variable with lazy array - var = xr.Variable(dims=[start_time_dim, "z", "y_3d", "x_3d"], data=lazy_array) - - # Extract time coordinates - start_times = [map_data.start_time for map_data in map_list] - end_times = [map_data.end_time for map_data in map_list] - - # Create coordinates - coordinates = { - "x_3d": x_coords, - "y_3d": y_coords, - "z": z_coords, - start_time_dim: start_times, - end_time_dim: (start_time_dim, end_times), - } - - # Convert to DataArray - data_array = xr.DataArray( - var, - coords=coordinates, - name=str3ds_name, - ) - - # Add CF attributes - r3_infos = grass_i.get_raster3d_info(map_list[0].id) - da_with_attrs = set_cf_coordinates( - data_array, - gi=grass_i, - is_3d=True, - z_unit=r3_infos["vertical_units"], - 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([r3_infos["source1"], r3_infos["source2"]]) - da_with_attrs.attrs["units"] = r3_infos.get("units", "") - da_with_attrs.attrs["comment"] = r3_infos.get("comments", "") - # CF attributes "institution" and "references" - # Do not correspond to a direct GRASS value. - return da_with_attrs - - -def set_cf_coordinates( - da: xr.DataArray, - gi: GrassInterface, - is_3d: bool, - z_unit: str = "", - time_dims: Optional[list[str, str]] = None, - time_unit: str = "", -): - """Set coordinate attributes according to CF conventions""" - spatial_unit = gi.get_spatial_units() - if is_3d: - # da["z"].attrs["positive"] = "up" # Not defined by grass - da["z"].attrs["axis"] = "Z" - da["z"].attrs["units"] = z_unit - y_coord = "y_3d" - x_coord = "x_3d" - else: - y_coord = "y" - x_coord = "x" - 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(): - da[x_coord].attrs["long_name"] = "longitude" - da[x_coord].attrs["units"] = "degrees_east" - da[x_coord].attrs["standard_name"] = "longitude" - da[y_coord].attrs["long_name"] = "latitude" - da[y_coord].attrs["units"] = "degrees_north" - da[y_coord].attrs["standard_name"] = "latitude" - else: - da[x_coord].attrs["long_name"] = "x-coordinate in Cartesian system" - da[y_coord].attrs["long_name"] = "y-coordinate in Cartesian system" - if gi.is_xy(): - da[x_coord].attrs["standard_name"] = "x_coordinate" - da[y_coord].attrs["standard_name"] = "y_coordinate" + crs_wkt = self.grass_interface.get_crs_wkt_str() + + coords_dict = {} + for coords in 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} + + # Set attributes + dirpath = Path(filename_or_obj) + mapset = dirpath.stem + project = dirpath.parent.stem + attrs = { + "crs_wkt": crs_wkt, + "Conventions": "CF-1.13-draft", + # "title": "", + "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": "", + } + dataset = xr.Dataset(data_vars=data_array_dict, coords=coords_dict, attrs=attrs) + return dataset + + def _set_cf_coordinates_attributes( + self, + da: xr.DataArray, + is_3d: bool, + z_unit: str = "", + time_dims: Optional[list[str, str]] = None, + time_unit: str = "", + ): + """Set coordinate attributes according to CF conventions""" + spatial_unit = self.grass_interface.get_spatial_units() + if is_3d: + # da["z"].attrs["positive"] = "up" # Not defined by grass + da["z"].attrs["axis"] = "Z" + da["z"].attrs["units"] = z_unit + y_coord = "y_3d" + x_coord = "x_3d" + else: + y_coord = "y" + x_coord = "x" + 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 self.grass_interface.is_latlon(): + da[x_coord].attrs["long_name"] = "longitude" + da[x_coord].attrs["units"] = "degrees_east" + da[x_coord].attrs["standard_name"] = "longitude" + da[y_coord].attrs["long_name"] = "latitude" + da[y_coord].attrs["units"] = "degrees_north" + da[y_coord].attrs["standard_name"] = "latitude" + else: + da[x_coord].attrs["long_name"] = "x-coordinate in Cartesian system" + da[y_coord].attrs["long_name"] = "y-coordinate in Cartesian system" + if self.grass_interface.is_xy(): + da[x_coord].attrs["standard_name"] = "x_coordinate" + da[y_coord].attrs["standard_name"] = "y_coordinate" + else: + da[x_coord].attrs["standard_name"] = "projection_x_coordinate" + da[y_coord].attrs["standard_name"] = "projection_y_coordinate" + da[x_coord].attrs["units"] = spatial_unit + da[y_coord].attrs["units"] = spatial_unit + return da + + def _open_grass_raster(self, raster_name: str) -> xr.DataArray: + """Open a single raster map.""" + x_coords, y_coords, _ = self.grass_interface.get_coordinates( + raster_3d=False + ).values() + dims = ["y", "x"] + coordinates = dict.fromkeys(dims) + coordinates["x"] = x_coords + coordinates["y"] = y_coords + raster_array = self.grass_interface.read_raster_map(raster_name) + data_array = xr.DataArray( + raster_array, + coords=coordinates, + dims=dims, + name=self.grass_interface.get_name_from_id(raster_name), + ) + # Add CF attributes + r_infos = self.grass_interface.get_raster_info(raster_name) + da_with_attrs = self._set_cf_coordinates_attributes(data_array, is_3d=False) + da_with_attrs.attrs["long_name"] = r_infos.get("title", "") + da_with_attrs.attrs["source"] = ",".join( + [r_infos["source1"], r_infos["source2"]] + ) + da_with_attrs.attrs["units"] = r_infos.get("units", "") + da_with_attrs.attrs["comment"] = r_infos.get("comments", "") + # CF attributes "institution" and "references" + # Do not correspond to a direct GRASS value. + return da_with_attrs + + def _open_grass_raster_3d(self, raster_3d_name: str) -> xr.DataArray: + """Open a single 3D raster map.""" + x_coords, y_coords, z_coords = self.grass_interface.get_coordinates( + raster_3d=True + ).values() + dims = ["z", "y_3d", "x_3d"] + coordinates = dict.fromkeys(dims) + coordinates["x_3d"] = x_coords + coordinates["y_3d"] = y_coords + coordinates["z"] = z_coords + raster_array = self.grass_interface.read_raster3d_map(raster_3d_name) + + data_array = xr.DataArray( + raster_array, + coords=coordinates, + dims=dims, + name=self.grass_interface.get_name_from_id(raster_3d_name), + ) + # Add CF attributes + r3_infos = self.grass_interface.get_raster3d_info(raster_3d_name) + da_with_attrs = self._set_cf_coordinates_attributes( + data_array, is_3d=True, z_unit=r3_infos["vertical_units"] + ) + da_with_attrs.attrs["long_name"] = r3_infos.get("title", "") + da_with_attrs.attrs["source"] = ",".join( + [r3_infos["source1"], r3_infos["source2"]] + ) + da_with_attrs.attrs["units"] = r3_infos.get("units", "") + da_with_attrs.attrs["comment"] = r3_infos.get("comments", "") + # CF attributes "institution" and "references" + # Do not correspond to a direct GRASS value. + return da_with_attrs + + def _open_grass_strds(self, strds_name: str) -> xr.DataArray: + """Open a STRDS with lazy loading - data is only loaded when accessed""" + strds_id = self.grass_interface.get_id_from_name(strds_name) + strds_name = self.grass_interface.get_name_from_id(strds_id) + x_coords, y_coords, _ = self.grass_interface.get_coordinates( + raster_3d=False + ).values() + strds_infos = self.grass_interface.get_stds_infos(strds_id, stds_type="strds") + if strds_infos.temporal_type == "absolute": + time_unit = None + else: + time_unit = strds_infos.time_unit + start_time_dim = f"start_time_{strds_name}" + end_time_dim = f"end_time_{strds_name}" + + map_list = self.grass_interface.list_maps_in_strds(strds_id) + region = self.grass_interface.get_region() + + # Create a single backend array for the entire STRDS + backend_array = GrassSTDSBackendArray( + shape=(len(map_list), region.rows, region.cols), + dtype=map_list[0].dtype, + map_list=map_list, + map_type="raster", + grass_interface=self.grass_interface, + ) + lazy_array = xr.core.indexing.LazilyIndexedArray(backend_array) + + # Create Variable with lazy array + var = xr.Variable(dims=[start_time_dim, "y", "x"], data=lazy_array) + + # Extract time coordinates + start_times = [map_data.start_time for map_data in map_list] + end_times = [map_data.end_time for map_data in map_list] + + # Create coordinates + coordinates = { + "x": x_coords, + "y": y_coords, + start_time_dim: start_times, + end_time_dim: (start_time_dim, end_times), + } + + # Convert to DataArray + data_array = xr.DataArray( + var, + coords=coordinates, + name=strds_name, + ) + + # Add CF attributes + r_infos = self.grass_interface.get_raster_info(map_list[0].id) + da_with_attrs = self._set_cf_coordinates_attributes( + data_array, + 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"]] + ) + da_with_attrs.attrs["units"] = r_infos.get("units", "") + da_with_attrs.attrs["comment"] = r_infos.get("comments", "") + # CF attributes "institution" and "references" + # Do not correspond to a direct GRASS value. + return da_with_attrs + + def _open_grass_str3ds(self, str3ds_name: str) -> xr.DataArray: + """Open a STR3DS with lazy loading - data is only loaded when accessed + TODO: Figure out what to do when the z value of the maps is time.""" + str3ds_id = self.grass_interface.get_id_from_name(str3ds_name) + str3ds_name = self.grass_interface.get_name_from_id(str3ds_id) + x_coords, y_coords, z_coords = self.grass_interface.get_coordinates( + raster_3d=True + ).values() + strds_infos = self.grass_interface.get_stds_infos(str3ds_id, stds_type="str3ds") + if strds_infos.temporal_type == "absolute": + time_unit = None else: - da[x_coord].attrs["standard_name"] = "projection_x_coordinate" - da[y_coord].attrs["standard_name"] = "projection_y_coordinate" - da[x_coord].attrs["units"] = spatial_unit - da[y_coord].attrs["units"] = spatial_unit - return da + time_unit = strds_infos.time_unit + start_time_dim = f"start_time_{str3ds_name}" + end_time_dim = f"end_time_{str3ds_name}" + + map_list = self.grass_interface.list_maps_in_str3ds(str3ds_id) + region = self.grass_interface.get_region() + + # Create a single backend array for the entire STR3DS + backend_array = GrassSTDSBackendArray( + shape=(len(map_list), region.depths, region.rows3, region.cols3), + dtype=map_list[0].dtype, + map_list=map_list, + map_type="raster3d", + grass_interface=self.grass_interface, + ) + lazy_array = xr.core.indexing.LazilyIndexedArray(backend_array) + + # Create Variable with lazy array + var = xr.Variable(dims=[start_time_dim, "z", "y_3d", "x_3d"], data=lazy_array) + + # Extract time coordinates + start_times = [map_data.start_time for map_data in map_list] + end_times = [map_data.end_time for map_data in map_list] + + # Create coordinates + coordinates = { + "x_3d": x_coords, + "y_3d": y_coords, + "z": z_coords, + start_time_dim: start_times, + end_time_dim: (start_time_dim, end_times), + } + + # Convert to DataArray + data_array = xr.DataArray( + var, + coords=coordinates, + name=str3ds_name, + ) + + # Add CF attributes + r3_infos = self.grass_interface.get_raster3d_info(map_list[0].id) + da_with_attrs = self._set_cf_coordinates_attributes( + data_array, + is_3d=True, + z_unit=r3_infos["vertical_units"], + 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( + [r3_infos["source1"], r3_infos["source2"]] + ) + da_with_attrs.attrs["units"] = r3_infos.get("units", "") + da_with_attrs.attrs["comment"] = r3_infos.get("comments", "") + # CF attributes "institution" and "references" + # Do not correspond to a direct GRASS value. + return da_with_attrs + + +def dir_is_grass_mapset(filename_or_obj: str | Path) -> bool: + """ + Check if the given path is a GRASS mapset. + """ + try: + dirpath = Path(filename_or_obj) + except TypeError: + return False + if dirpath.is_dir(): + wind_file = dirpath / Path("WIND") + # A newly created mapset might only have WIND, VAR appears later. + if wind_file.exists(): + return True + return False + + +def dir_is_grass_project(filename_or_obj: str | Path) -> bool: + """Return True if a subdir named PERMANENT is present.""" + try: + dirpath = Path(filename_or_obj) + except TypeError: + return False + if dirpath.is_dir(): + return (dirpath / Path("PERMANENT")).is_dir() + else: + return False diff --git a/tests/test_tograss.py b/tests/test_tograss.py index 47e82ad..55ea024 100644 --- a/tests/test_tograss.py +++ b/tests/test_tograss.py @@ -68,10 +68,7 @@ def test_dataarray_2d_conversion( ) target_mapset_name = temp_gisdb.mapset # Use PERMANENT mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = mapset_path_obj if mapset_is_path_obj else str(mapset_path_obj) + grass_raster_name_full = ( f"{sample_da.name}@{target_mapset_name}" # Moved and defined ) @@ -79,8 +76,6 @@ def test_dataarray_2d_conversion( try: to_grass( dataset=sample_da, - mapset=mapset_arg, - create=False, ) available_rasters = grass_i.list_raster(mapset=target_mapset_name) assert grass_raster_name_full in available_rasters, ( @@ -175,17 +170,11 @@ def test_dataarray_3d_conversion( ) target_mapset_name = temp_gisdb.mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = mapset_path_obj if mapset_is_path_obj else str(mapset_path_obj) # Try statement for file cleanup try: to_grass( dataset=sample_da, - mapset=mapset_arg, - create=False, ) grass_raster_name_full = f"{sample_da.name}@{target_mapset_name}" available_rasters_3d = grass_i.list_raster3d(mapset=target_mapset_name) @@ -281,15 +270,9 @@ def test_dataarray_to_strds_conversion( sample_da["time"].attrs["units"] = "days" target_mapset_name = temp_gisdb.mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = mapset_path_obj if mapset_is_path_obj else str(mapset_path_obj) to_grass( dataset=sample_da, - mapset=mapset_arg, - create=False, dims={"test_strds": {"start_time": "time"}}, ) strds_id = f"{sample_da.name}@{target_mapset_name}" @@ -426,16 +409,8 @@ def test_dataarray_to_str3ds_conversion( if time_dim_type == "relative": sample_da["time"].attrs["units"] = "days" - target_mapset_name = temp_gisdb.mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = mapset_path_obj if mapset_is_path_obj else str(mapset_path_obj) - to_grass( dataset=sample_da, - mapset=mapset_arg, - create=False, dims={"test_str3ds_vol": {"start_time": "time"}}, ) try: @@ -526,11 +501,6 @@ def test_dataset_conversion_mixed_types( ): """Test conversion of an xr.Dataset with mixed DataArray types.""" session_crs_wkt = grass_i.get_crs_wkt_str() - target_mapset_name = temp_gisdb.mapset # Use PERMANENT mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = mapset_path_obj if mapset_is_path_obj else str(mapset_path_obj) # 2D Raster Spec da_2d_spec = { @@ -592,8 +562,6 @@ def test_dataset_conversion_mixed_types( to_grass( dataset=sample_ds, - mapset=mapset_arg, - create=False, dims={ "strds_var": {"start_time": "time"}, "str3ds_var": {"start_time": "time"}, @@ -713,7 +681,7 @@ def test_mapset_creation_true(self, temp_gisdb, grass_i: GrassInterface): name="data_for_mapset_creation", ) - to_grass(dataset=sample_da, mapset=str(mapset_path), create=True) + to_grass(dataset=sample_da) assert mapset_path.exists() and mapset_path.is_dir(), ( f"Mapset directory {mapset_path} was not created." @@ -788,7 +756,7 @@ def test_mapset_creation_false_existing_mapset( name="data_for_existing_mapset", ) - to_grass(dataset=sample_da, mapset=str(mapset_path), create=False) + to_grass(dataset=sample_da) available_rasters = grass_i.list_raster(mapset=existing_mapset_name) assert sample_da.name in available_rasters, ( @@ -812,7 +780,6 @@ def test_dims_mapping( """Test 'dims' mapping functionality.""" session_crs_wkt = grass_i.get_crs_wkt_str() target_mapset_name = temp_gisdb.mapset - mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name da_name = "dims_test_raster" img_height, img_width = 3, 2 @@ -832,8 +799,6 @@ def test_dims_mapping( sample_da = sample_da.rename(rename_map[da_name]) to_grass( dataset=sample_da, - mapset=str(mapset_path), - create=False, dims=rename_map, ) available_rasters = grass_i.list_raster(mapset=target_mapset_name) @@ -856,11 +821,6 @@ def test_dimension_transposition( and verifies they are correctly transposed when written to GRASS. """ session_crs_wkt = grass_i.get_crs_wkt_str() - target_mapset_name = temp_gisdb.mapset - mapset_path_obj = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / target_mapset_name - ) - mapset_arg = str(mapset_path_obj) # Define expected dimensions for verification height_2d, width_2d = 5, 7 @@ -971,22 +931,20 @@ def test_dimension_transposition( try: # Write 2D raster - to_grass(dataset=da_2d, mapset=mapset_arg) + to_grass(dataset=da_2d) # Write 3D raster - to_grass(dataset=da_3d, mapset=mapset_arg) + to_grass(dataset=da_3d) # Write STRDS to_grass( dataset=da_strds, - mapset=mapset_arg, dims={strds_name: {"start_time": "time"}}, ) # Write STR3DS to_grass( dataset=da_str3ds, - mapset=mapset_arg, dims={str3ds_name: {"start_time": "time"}}, ) diff --git a/tests/test_tograss_error_handling.py b/tests/test_tograss_error_handling.py index 9538eed..622c019 100644 --- a/tests/test_tograss_error_handling.py +++ b/tests/test_tograss_error_handling.py @@ -14,7 +14,6 @@ """ from pathlib import Path -import tempfile import xarray as xr import numpy as np @@ -28,8 +27,6 @@ class TestToGrassErrorHandling: def test_missing_crs_wkt_attribute(self, temp_gisdb, grass_i: GrassInterface): """Test error handling when input xarray object is missing 'crs_wkt' attribute.""" - mapset_name = temp_gisdb.mapset - mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / mapset_name # Create a DataArray without crs_wkt # The helper function always adds it, so we create one manually here. sample_da_no_crs = xr.DataArray( @@ -44,12 +41,10 @@ def test_missing_crs_wkt_attribute(self, temp_gisdb, grass_i: GrassInterface): (KeyError, AttributeError, ValueError), match=r"(crs_wkt|CRS mismatch|has no attribute 'attrs')", ): - to_grass(dataset=sample_da_no_crs, mapset=str(mapset_path), create=False) + to_grass(dataset=sample_da_no_crs) def test_incompatible_crs_wkt(self, temp_gisdb, grass_i: GrassInterface): """Test error handling with an incompatible 'crs_wkt' attribute.""" - mapset_name = temp_gisdb.mapset - mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / mapset_name session_crs_wkt = grass_i.get_crs_wkt_str() # Create an incompatible CRS WKT string incompatible_crs = CRS.from_epsg(4326) # WGS 84 @@ -69,113 +64,7 @@ def test_incompatible_crs_wkt(self, temp_gisdb, grass_i: GrassInterface): ValueError, match=r"CRS mismatch", ): - to_grass(dataset=sample_da, mapset=str(mapset_path), create=False) - - def test_invalid_mapset_path_non_existent_parent( - self, temp_gisdb, grass_i: GrassInterface - ): - """Test error with mapset path having a non-existent parent directory.""" - pytest.skip("Skipping mapset creation test due to GRASS <8.5 tgis.init() bug.") - session_crs_wkt = grass_i.get_crs_wkt_str() - # Path to a non-existent directory, then the mapset - mapset_path = ( - Path(temp_gisdb.gisdb) - / temp_gisdb.project - / "non_existent_parent_dir" - / "my_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_invalid_parent", - ) - - with pytest.raises( - ValueError, - match=r"Mapset.*not found and its parent directory.*is not a valid GRASS project", - ): - to_grass(dataset=sample_da, mapset=str(mapset_path), create=True) - - def test_invalid_mapset_path_is_file(self, temp_gisdb, grass_i: GrassInterface): - """Test error with mapset path being an existing file.""" - pytest.skip("Skipping mapset creation test due to GRASS <8.5 tgis.init() bug.") - session_crs_wkt = grass_i.get_crs_wkt_str() - - # Create an empty file where the mapset directory would be - file_as_mapset_path = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / "file_instead_of_mapset" - ) - with open(file_as_mapset_path, "w") as f: - f.write("This is a file.") - - assert file_as_mapset_path.is_file() - - 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_mapset_is_file", - ) - - with pytest.raises(ValueError, match=r"not a directory"): - to_grass(dataset=sample_da, mapset=str(file_as_mapset_path), create=True) - - file_as_mapset_path.unlink() # Clean up the created file - - def test_parent_dir_not_grass_location(self, grass_i: GrassInterface): - """Test error when parent of mapset is not a GRASS Location (create=True).""" - pytest.skip("Skipping mapset creation test due to GRASS <8.5 tgis.init() bug.") - session_crs_wkt = grass_i.get_crs_wkt_str() - - with tempfile.TemporaryDirectory( - prefix="not_a_grass_loc_" - ) as tmp_non_grass_dir: - mapset_path_in_non_grass_loc = Path(tmp_non_grass_dir) / "my_mapset_here" - - 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_non_grass_parent", - ) - - # This relies on to_grass checking if the parent is a GRASS location. - # The exact error message might vary based on implementation. - with pytest.raises( - ValueError, - match=r"(not a valid GRASS project|Parent directory.*not a GRASS location|Invalid GIS database)", - ): - to_grass( - dataset=sample_da, - mapset=str(mapset_path_in_non_grass_loc), - create=True, - ) - - def test_create_false_mapset_not_exists(self, temp_gisdb, grass_i: GrassInterface): - """Test error when create=False and mapset does not exist.""" - session_crs_wkt = grass_i.get_crs_wkt_str() - non_existent_mapset_path = ( - Path(temp_gisdb.gisdb) / temp_gisdb.project / "mapset_does_not_exist_at_all" - ) - - assert not non_existent_mapset_path.exists() # Ensure it really doesn't exist - - 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_create_false_no_mapset", - ) - - with pytest.raises( - ValueError, - match=r"is not a valid directory", - ): - to_grass( - dataset=sample_da, mapset=str(non_existent_mapset_path), create=False - ) + to_grass(dataset=sample_da) def test_mapset_not_accessible_simplified(self, grass_i: GrassInterface): """Test simplified 'mapset not accessible' by providing a syntactically valid but unrelated path.""" @@ -210,10 +99,8 @@ def test_mapset_not_accessible_simplified(self, grass_i: GrassInterface): ValueError, match=r"not found and .* is not a valid GRASS project", ): - to_grass(dataset=sample_da, mapset=unrelated_path, create=True) - to_grass( - dataset=sample_da, mapset=unrelated_path, create=False - ) # Also test with create=False + to_grass(dataset=sample_da) + to_grass(dataset=sample_da) @pytest.mark.usefixtures("grass_session_fixture") @@ -221,22 +108,17 @@ class TestToGrassInputValidation: def test_invalid_dataset_type(self, temp_gisdb, grass_i: GrassInterface): """Test error handling for invalid 'dataset' parameter type. That a first try. Let's see how it goes considering that the tested code uses duck typing.""" - mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / temp_gisdb.mapset - invalid_datasets = [123, "a string", [1, 2, 3], {"data": np.array([1])}, None] for invalid_ds in invalid_datasets: with pytest.raises( TypeError, match=r"'dataset must be either an Xarray DataArray or Dataset", ): - to_grass(dataset=invalid_ds, mapset=str(mapset_path), create=False) + to_grass(dataset=invalid_ds) def test_invalid_dims_parameter_type(self, temp_gisdb, grass_i: GrassInterface): """Test error handling for invalid 'dims' parameter type or content.""" session_crs_wkt = grass_i.get_crs_wkt_str() - # Use PERMANENT mapset to avoid issues with tgis.init() for newly created mapsets - mapset_path = Path(temp_gisdb.gisdb) / temp_gisdb.project / temp_gisdb.mapset - # No need to create PERMANENT mapset as it's guaranteed by temp_gisdb fixture sample_da = create_sample_dataarray( dims_spec={"y": np.arange(2.0), "x": np.arange(2.0)}, @@ -257,37 +139,11 @@ def test_invalid_dims_parameter_type(self, temp_gisdb, grass_i: GrassInterface): ): to_grass( dataset=sample_da, - mapset=str(mapset_path), dims=invalid_dims, - create=False, ) - def test_invalid_mapset_parameter_type(self, temp_gisdb, grass_i: GrassInterface): - """Test error handling for invalid 'mapset' parameter type.""" - session_crs_wkt = grass_i.get_crs_wkt_str() - 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_mapset_type_validation", - ) - - invalid_mapset_params = [ - 123, - None, - ["path", "to", "mapset"], - {"path": "mapset_dir"}, - ] - for invalid_mapset in invalid_mapset_params: - with pytest.raises( - TypeError, - 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), @@ -300,7 +156,5 @@ def test_invalid_dims_invalid_var(self, temp_gisdb, grass_i: GrassInterface): ): to_grass( dataset=sample_da, - mapset=str(mapset_path), dims={"invalid_var_name": {"x": "my_x"}}, - create=False, ) diff --git a/uv.lock b/uv.lock index 2d527a1..8007a42 100644 --- a/uv.lock +++ b/uv.lock @@ -493,7 +493,6 @@ wheels = [ name = "xarray-grass" source = { editable = "." } dependencies = [ - { name = "grass-session" }, { name = "numpy" }, { name = "pandas" }, { name = "pyproj" }, @@ -502,6 +501,7 @@ dependencies = [ [package.dev-dependencies] dev = [ + { name = "grass-session" }, { name = "pre-commit" }, { name = "pytest" }, { name = "pytest-random-order" }, @@ -511,7 +511,6 @@ dev = [ [package.metadata] requires-dist = [ - { name = "grass-session", specifier = ">=0.5" }, { name = "numpy", specifier = ">=2.2.5" }, { name = "pandas", specifier = ">=2.2.3" }, { name = "pyproj", specifier = ">=3.7.1" }, @@ -520,6 +519,7 @@ requires-dist = [ [package.metadata.requires-dev] dev = [ + { name = "grass-session", specifier = ">=0.5" }, { name = "pre-commit", specifier = ">=4.2.0" }, { name = "pytest", specifier = ">=8.3.5" }, { name = "pytest-random-order", specifier = ">=1.1.1" },