Data I/O and GISData container

This module provides consolidated I/O for preprocessed MOBIDIC data and a container class for managing the complete preprocessed dataset.

Overview

After preprocessing GIS data, processing the river network, and optionally processing reservoirs, MOBIDICpy consolidates all data into a single GISData object that can be saved to and loaded from disk. This approach:

Simplifies data management: Single object contains all preprocessed data
Ensures consistency: Grids, network, and metadata stay synchronized
Enables caching: Save expensive preprocessing results for reuse
Facilitates sharing: Package preprocessed data for model runs

Classes

GISData container

Container for preprocessed GIS data.

This class holds all preprocessed spatial data including grids, river network, reservoirs, and hillslope-reach mapping. It provides methods to save/load consolidated preprocessed data.

Attributes:

Name	Type	Description
`grids`		Dictionary of 2D numpy arrays containing raster data
`metadata`		Dictionary containing grid metadata (transform, CRS, resolution, etc.)
`network`		GeoDataFrame with processed river network
`reservoirs`		Reservoirs object with reservoir data (optional)
`hillslope_reach_map`		2D array mapping each cell to its downstream reach
`config`		MOBIDIC configuration used for preprocessing

Source code in mobidic/preprocessing/preprocessor.py

class GISData:
    """Container for preprocessed GIS data.

    This class holds all preprocessed spatial data including grids, river network,
    reservoirs, and hillslope-reach mapping. It provides methods to save/load
    consolidated preprocessed data.

    Attributes:
        grids: Dictionary of 2D numpy arrays containing raster data
        metadata: Dictionary containing grid metadata (transform, CRS, resolution, etc.)
        network: GeoDataFrame with processed river network
        reservoirs: Reservoirs object with reservoir data (optional)
        hillslope_reach_map: 2D array mapping each cell to its downstream reach
        config: MOBIDIC configuration used for preprocessing
    """

    def __init__(
        self,
        grids: dict[str, np.ndarray],
        metadata: dict[str, Any],
        network: gpd.GeoDataFrame,
        hillslope_reach_map: np.ndarray,
        config: MOBIDICConfig,
        reservoirs: Optional[Reservoirs] = None,
    ):
        """Initialize GISData container.

        Args:
            grids: Dictionary of 2D numpy arrays with grid data
            metadata: Dictionary with grid metadata
            network: Processed river network GeoDataFrame
            hillslope_reach_map: 2D array with reach assignments
            config: MOBIDIC configuration
            reservoirs: Optional Reservoirs object
        """
        self.grids = grids
        self.metadata = metadata
        self.network = network
        self.reservoirs = reservoirs
        self.hillslope_reach_map = hillslope_reach_map
        self.config = config

    def save(
        self, gisdata_path: str | Path, network_path: str | Path, reservoirs_path: Optional[str | Path] = None
    ) -> None:
        """Save preprocessed data to files.

        Args:
            gisdata_path: Path to save gridded data (NetCDF format)
            network_path: Path to save river network (GeoParquet format)
            reservoirs_path: Optional path to save reservoirs (GeoParquet format)
        """
        from mobidic.preprocessing.io import save_gisdata, save_network, save_reservoirs

        save_gisdata(self, gisdata_path)
        save_network(self.network, network_path, format="parquet")

        if self.reservoirs is not None and reservoirs_path is not None:
            save_reservoirs(self.reservoirs, reservoirs_path, format="parquet")

    @classmethod
    def load(
        cls, gisdata_path: str | Path, network_path: str | Path, reservoirs_path: Optional[str | Path] = None
    ) -> "GISData":
        """Load preprocessed data from files.

        Args:
            gisdata_path: Path to gridded data file (NetCDF)
            network_path: Path to river network file (GeoParquet)
            reservoirs_path: Optional path to reservoirs file (GeoParquet)

        Returns:
            GISData object with loaded data
        """
        from mobidic.preprocessing.io import load_gisdata, load_reservoirs

        gisdata = load_gisdata(gisdata_path, network_path)

        # Load reservoirs if path provided
        if reservoirs_path is not None:
            reservoirs_path = Path(reservoirs_path)
            if reservoirs_path.exists():
                gisdata.reservoirs = load_reservoirs(reservoirs_path)
            else:
                logger.warning(f"Reservoirs file not found: {reservoirs_path}")
                gisdata.reservoirs = None

        return gisdata

`init(grids, metadata, network, hillslope_reach_map, config, reservoirs=None)`

Initialize GISData container.

Parameters:

Name	Type	Description	Default
`grids`	`dict[str, ndarray]`	Dictionary of 2D numpy arrays with grid data	required
`metadata`	`dict[str, Any]`	Dictionary with grid metadata	required
`network`	`GeoDataFrame`	Processed river network GeoDataFrame	required
`hillslope_reach_map`	`ndarray`	2D array with reach assignments	required
`config`	`MOBIDICConfig`	MOBIDIC configuration	required
`reservoirs`	`Optional[Reservoirs]`	Optional Reservoirs object	`None`

Source code in mobidic/preprocessing/preprocessor.py

def __init__(
    self,
    grids: dict[str, np.ndarray],
    metadata: dict[str, Any],
    network: gpd.GeoDataFrame,
    hillslope_reach_map: np.ndarray,
    config: MOBIDICConfig,
    reservoirs: Optional[Reservoirs] = None,
):
    """Initialize GISData container.

    Args:
        grids: Dictionary of 2D numpy arrays with grid data
        metadata: Dictionary with grid metadata
        network: Processed river network GeoDataFrame
        hillslope_reach_map: 2D array with reach assignments
        config: MOBIDIC configuration
        reservoirs: Optional Reservoirs object
    """
    self.grids = grids
    self.metadata = metadata
    self.network = network
    self.reservoirs = reservoirs
    self.hillslope_reach_map = hillslope_reach_map
    self.config = config

`load(gisdata_path, network_path, reservoirs_path=None)` `classmethod`

Load preprocessed data from files.

Parameters:

Name	Type	Description	Default
`gisdata_path`	`str \| Path`	Path to gridded data file (NetCDF)	required
`network_path`	`str \| Path`	Path to river network file (GeoParquet)	required
`reservoirs_path`	`Optional[str \| Path]`	Optional path to reservoirs file (GeoParquet)	`None`

Returns:

Type	Description
`GISData`	GISData object with loaded data

Source code in mobidic/preprocessing/preprocessor.py

@classmethod
def load(
    cls, gisdata_path: str | Path, network_path: str | Path, reservoirs_path: Optional[str | Path] = None
) -> "GISData":
    """Load preprocessed data from files.

    Args:
        gisdata_path: Path to gridded data file (NetCDF)
        network_path: Path to river network file (GeoParquet)
        reservoirs_path: Optional path to reservoirs file (GeoParquet)

    Returns:
        GISData object with loaded data
    """
    from mobidic.preprocessing.io import load_gisdata, load_reservoirs

    gisdata = load_gisdata(gisdata_path, network_path)

    # Load reservoirs if path provided
    if reservoirs_path is not None:
        reservoirs_path = Path(reservoirs_path)
        if reservoirs_path.exists():
            gisdata.reservoirs = load_reservoirs(reservoirs_path)
        else:
            logger.warning(f"Reservoirs file not found: {reservoirs_path}")
            gisdata.reservoirs = None

    return gisdata

`save(gisdata_path, network_path, reservoirs_path=None)`

Save preprocessed data to files.

Parameters:

Name	Type	Description	Default
`gisdata_path`	`str \| Path`	Path to save gridded data (NetCDF format)	required
`network_path`	`str \| Path`	Path to save river network (GeoParquet format)	required
`reservoirs_path`	`Optional[str \| Path]`	Optional path to save reservoirs (GeoParquet format)	`None`

Source code in mobidic/preprocessing/preprocessor.py

def save(
    self, gisdata_path: str | Path, network_path: str | Path, reservoirs_path: Optional[str | Path] = None
) -> None:
    """Save preprocessed data to files.

    Args:
        gisdata_path: Path to save gridded data (NetCDF format)
        network_path: Path to save river network (GeoParquet format)
        reservoirs_path: Optional path to save reservoirs (GeoParquet format)
    """
    from mobidic.preprocessing.io import save_gisdata, save_network, save_reservoirs

    save_gisdata(self, gisdata_path)
    save_network(self.network, network_path, format="parquet")

    if self.reservoirs is not None and reservoirs_path is not None:
        save_reservoirs(self.reservoirs, reservoirs_path, format="parquet")

Functions

Save and load functions

Save preprocessed gridded data to NetCDF file.

This function saves all grids (DTM, flow direction, soil parameters, etc.), metadata, and hillslope-reach mapping to a single NetCDF file using xarray.

Parameters:

Name	Type	Description	Default
`gisdata`	`GISData`	GISData object containing preprocessed data	required
`output_path`	`str \| Path`	Path to output NetCDF file	required

Examples:

>>> from mobidic import run_preprocessing, load_config
>>> config = load_config("Arno.yaml")
>>> gisdata = run_preprocessing(config)
>>> from mobidic.preprocessing.io import save_gisdata
>>> save_gisdata(gisdata, "Arno_gisdata.nc")

Source code in mobidic/preprocessing/io.py

def save_gisdata(gisdata: "GISData", output_path: str | Path) -> None:
    """Save preprocessed gridded data to NetCDF file.

    This function saves all grids (DTM, flow direction, soil parameters, etc.),
    metadata, and hillslope-reach mapping to a single NetCDF file using xarray.

    Args:
        gisdata: GISData object containing preprocessed data
        output_path: Path to output NetCDF file

    Examples:
        >>> from mobidic import run_preprocessing, load_config
        >>> config = load_config("Arno.yaml")
        >>> gisdata = run_preprocessing(config)
        >>> from mobidic.preprocessing.io import save_gisdata
        >>> save_gisdata(gisdata, "Arno_gisdata.nc")
    """
    output_path = Path(output_path)
    output_path.parent.mkdir(parents=True, exist_ok=True)

    logger.info(f"Saving gridded data to NetCDF: {output_path}")

    # Get dimensions from metadata
    nrows, ncols = gisdata.metadata["shape"]
    resolution = gisdata.metadata["resolution"]

    # Create coordinate arrays
    # Use cell centers following MOBIDIC convention (xllcorner and yllcorner are already at cell centers)
    x = gisdata.metadata["xllcorner"] + np.arange(ncols) * resolution[0]
    y = gisdata.metadata["yllcorner"] + np.arange(nrows) * resolution[1]

    # Create xarray Dataset
    data_vars = {}

    # Add all grids as data variables
    for name, grid in gisdata.grids.items():
        data_vars[name] = (["y", "x"], grid)

    # Add hillslope-reach mapping
    data_vars["hillslope_reach_map"] = (["y", "x"], gisdata.hillslope_reach_map)

    # Add grid mapping variable for CRS (CF-1.12 compliance)
    data_vars["crs"] = ([], 0)  # Scalar variable to hold CRS information

    # Create dataset
    ds = xr.Dataset(
        data_vars=data_vars,
        coords={
            "x": (["x"], x),
            "y": (["y"], y),
        },
    )

    # Add grid mapping variable attributes (CF-1.12 compliance)
    ds["crs"].attrs = crs_to_cf_attrs(gisdata.metadata["crs"])

    # Add coordinate variable attributes (CF-1.12 compliance)
    ds["x"].attrs = {
        "standard_name": "projection_x_coordinate",
        "long_name": "x coordinate of projection",
        "units": "m",
        "axis": "X",
    }
    ds["y"].attrs = {
        "standard_name": "projection_y_coordinate",
        "long_name": "y coordinate of projection",
        "units": "m",
        "axis": "Y",
    }

    # Add metadata as attributes
    ds.attrs["basin_id"] = gisdata.config.basin.id
    ds.attrs["paramset_id"] = gisdata.config.basin.paramset_id
    ds.attrs["basin_baricenter_lon"] = gisdata.config.basin.baricenter.lon
    ds.attrs["basin_baricenter_lat"] = gisdata.config.basin.baricenter.lat
    ds.attrs["resolution_x"] = resolution[0]
    ds.attrs["resolution_y"] = resolution[1]
    ds.attrs["decimation_factor"] = gisdata.config.simulation.decimation
    ds.attrs["flow_dir_notation"] = "MOBIDIC"
    ds.attrs["Conventions"] = "CF-1.12"
    ds.attrs["title"] = "MOBIDIC preprocessed GIS data"
    ds.attrs["source"] = "MOBIDICpy preprocessing"
    ds.attrs["history"] = f"Created by MOBIDICpy version {__version__}"

    # Add variable metadata
    ds["dtm"].attrs = {
        "long_name": "Digital Terrain Model",
        "units": "m",
        "description": "Elevation above sea level",
        "grid_mapping": "crs",
    }
    ds["flow_dir"].attrs = {
        "long_name": "Flow Direction",
        "units": "1",
        "description": "Flow direction in MOBIDIC notation (1-8)",
        "notation": "MOBIDIC: 1=NE, 2=N, 3=NW, 4=W, 5=SW, 6=S, 7=SE, 8=E",
        "grid_mapping": "crs",
    }
    ds["flow_acc"].attrs = {
        "long_name": "Flow Accumulation",
        "units": "cells",
        "description": "Number of upstream cells draining to each cell",
        "grid_mapping": "crs",
    }
    ds["Wc0"].attrs = {
        "long_name": "Capillary Water Holding Capacity",
        "units": "m",
        "description": "Maximum water holding capacity in soil small pores",
        "grid_mapping": "crs",
    }
    ds["Wg0"].attrs = {
        "long_name": "Gravitational Water Holding Capacity",
        "units": "m",
        "description": "Maximum water holding capacity in soil large pores",
        "grid_mapping": "crs",
    }
    ds["ks"].attrs = {
        "long_name": "Hydraulic Conductivity",
        "units": "m/s",
        "description": "Soil hydraulic conductivity",
        "grid_mapping": "crs",
    }
    ds["hillslope_reach_map"].attrs = {
        "long_name": "Hillslope to Reach Mapping",
        "units": "1",
        "description": "MOBIDIC ID of downstream reach for each cell",
        "grid_mapping": "crs",
        "_FillValue": np.nan,
    }

    ds["alpsur"].attrs = {
        "long_name": "Surface Routing Coefficient",
        "units": "1",
        "description": "Coefficient for surface routing",
        "grid_mapping": "crs",
        "_FillValue": np.nan,
    }

    # Add optional parameter metadata
    if "kf" in gisdata.grids:
        ds["kf"].attrs = {
            "long_name": "Aquifer Conductivity",
            "units": "m/s",
            "description": "Aquifer hydraulic conductivity",
            "grid_mapping": "crs",
        }
    if "CH" in gisdata.grids:
        ds["CH"].attrs = {
            "long_name": "Heat Exchange Coefficient",
            "units": "1",
            "description": "Turbulent exchange coefficient for heat",
            "grid_mapping": "crs",
        }
    if "Alb" in gisdata.grids:
        ds["Alb"].attrs = {
            "long_name": "Albedo",
            "units": "1",
            "description": "Surface albedo",
            "grid_mapping": "crs",
        }
    if "gamma" in gisdata.grids:
        ds["gamma"].attrs = {
            "long_name": "Percolation Coefficient",
            "units": "1/s",
            "description": "Percolation coefficient from capillary to gravitational reservoir",
            "grid_mapping": "crs",
        }
    if "kappa" in gisdata.grids:
        ds["kappa"].attrs = {
            "long_name": "Adsorption Coefficient",
            "units": "1/s",
            "description": "Adsorption coefficient from gravitational to capillary reservoir",
            "grid_mapping": "crs",
        }
    if "beta" in gisdata.grids:
        ds["beta"].attrs = {
            "long_name": "Hypodermic Flow Coefficient",
            "units": "1/s",
            "description": "Hypodermic (subsurface) flow coefficient",
            "grid_mapping": "crs",
        }
    if "alpha" in gisdata.grids:
        ds["alpha"].attrs = {
            "long_name": "Hillslope Flow Coefficient",
            "units": "1/s",
            "description": "Surface hillslope flow coefficient",
            "grid_mapping": "crs",
        }
    if "Ma" in gisdata.grids:
        ds["Ma"].attrs = {
            "long_name": "Artesian Aquifer Mask",
            "units": "1",
            "description": "Binary mask defining artesian aquifer extension (0=no, 1=yes)",
            "grid_mapping": "crs",
        }
    if "Mf" in gisdata.grids:
        ds["Mf"].attrs = {
            "long_name": "Freatic Aquifer Mask",
            "units": "1",
            "description": "Binary mask defining freatic aquifer extension (0=no, 1=yes)",
            "grid_mapping": "crs",
        }

    # Save to NetCDF
    encoding = {var: {"zlib": True, "complevel": 4} for var in ds.data_vars}
    ds.to_netcdf(output_path, encoding=encoding, engine="netcdf4")

    logger.success(f"Gridded data saved to {output_path}")
    logger.debug(f"File size: {output_path.stat().st_size / 1024 / 1024:.2f} MB")

Load preprocessed GIS data from NetCDF and GeoParquet files.

Parameters:

Name	Type	Description	Default
`gisdata_path`	`str \| Path`	Path to gridded data NetCDF file	required
`network_path`	`str \| Path`	Path to river network GeoParquet file	required

Returns:

Type	Description
`GISData`	GISData object containing loaded data

Raises:

Type	Description
`FileNotFoundError`	If either file does not exist

Examples:

>>> from mobidic.preprocessing.io import load_gisdata
>>> gisdata = load_gisdata("Arno_gisdata.nc", "Arno_network.parquet")

Source code in mobidic/preprocessing/io.py

def load_gisdata(gisdata_path: str | Path, network_path: str | Path) -> "GISData":
    """Load preprocessed GIS data from NetCDF and GeoParquet files.

    Args:
        gisdata_path: Path to gridded data NetCDF file
        network_path: Path to river network GeoParquet file

    Returns:
        GISData object containing loaded data

    Raises:
        FileNotFoundError: If either file does not exist

    Examples:
        >>> from mobidic.preprocessing.io import load_gisdata
        >>> gisdata = load_gisdata("Arno_gisdata.nc", "Arno_network.parquet")
    """
    from mobidic.preprocessing.preprocessor import GISData

    gisdata_path = Path(gisdata_path)
    network_path = Path(network_path)

    # Check files exist
    if not gisdata_path.exists():
        raise FileNotFoundError(f"GIS data file not found: {gisdata_path}")
    if not network_path.exists():
        raise FileNotFoundError(f"Network file not found: {network_path}")

    logger.info(f"Loading gridded data from NetCDF: {gisdata_path}")

    # Load NetCDF dataset
    ds = xr.open_dataset(gisdata_path)

    # Extract grids
    grids = {}
    grid_vars = [
        "dtm",
        "flow_dir",
        "flow_acc",
        "Wc0",
        "Wg0",
        "ks",
        "kf",
        "CH",
        "Alb",
        "Ma",
        "Mf",
        "gamma",
        "kappa",
        "beta",
        "alpha",
        "alpsur",
    ]

    for var in grid_vars:
        if var in ds:
            grids[var] = ds[var].values

    # Extract hillslope-reach mapping
    hillslope_reach_map = ds["hillslope_reach_map"].values

    # Extract metadata
    # Get CRS from grid mapping variable (CF-1.12 compliant)
    crs_value = ds["crs"].attrs.get("crs_wkt") if "crs" in ds else ds.attrs.get("crs")

    # Reconstruct xllcorner and yllcorner from coordinates
    resolution_x = ds.attrs.get("resolution_x")
    resolution_y = ds.attrs.get("resolution_y")
    xllcorner = float(ds.x.values[0])
    yllcorner = float(ds.y.values[0])

    metadata = {
        "shape": (len(ds.y), len(ds.x)),
        "resolution": (resolution_x, resolution_y),
        "crs": crs_value,
        "nodata": np.nan,  # CF-1.12: use variable-level _FillValue instead
        "xllcorner": xllcorner,
        "yllcorner": yllcorner,
        "bounds": None,  # Not stored, can reconstruct if needed
        "transform": None,  # Not stored, can reconstruct if needed
    }

    # Close dataset
    ds.close()

    logger.success(f"Gridded data loaded: {len(grids)} grids")

    # Load river network
    network = load_network(network_path)

    # Create GISData object
    # Note: config will be None since we only have partial config from attributes
    gisdata = GISData(
        grids=grids,
        metadata=metadata,
        network=network,
        hillslope_reach_map=hillslope_reach_map,
        config=None,  # type: ignore
    )

    return gisdata

Save processed river network to file.

Saves the river network to either GeoParquet (default, recommended) or Shapefile format. GeoParquet is more efficient, has no field name limitations, and preserves data types better.

Parameters:

Name	Type	Description	Default
`network`	`GeoDataFrame`	Processed river network GeoDataFrame	required
`output_path`	`str \| Path`	Path to output file	required
`format`	`str`	Output format, either ‘parquet’ (default) or ‘shapefile’	`'parquet'`

Raises:

Type	Description
`ValueError`	If format is not supported

Examples:

>>> from mobidic import process_river_network
>>> network = process_river_network("river_network.shp")
>>> from mobidic.preprocessing.io import save_network
>>> save_network(network, "river_network.parquet")  # Default: parquet
>>> save_network(network, "river_network.shp", format="shapefile")  # Shapefile

Source code in mobidic/preprocessing/io.py

def save_network(network: gpd.GeoDataFrame, output_path: str | Path, format: str = "parquet") -> None:
    """Save processed river network to file.

    Saves the river network to either GeoParquet (default, recommended) or Shapefile format.
    GeoParquet is more efficient, has no field name limitations, and preserves data types better.

    Args:
        network: Processed river network GeoDataFrame
        output_path: Path to output file
        format: Output format, either 'parquet' (default) or 'shapefile'

    Raises:
        ValueError: If format is not supported

    Examples:
        >>> from mobidic import process_river_network
        >>> network = process_river_network("river_network.shp")
        >>> from mobidic.preprocessing.io import save_network
        >>> save_network(network, "river_network.parquet")  # Default: parquet
        >>> save_network(network, "river_network.shp", format="shapefile")  # Shapefile
    """
    output_path = Path(output_path)
    output_path.parent.mkdir(parents=True, exist_ok=True)

    if format == "parquet":
        logger.info(f"Saving river network to GeoParquet: {output_path}")
        network.to_parquet(output_path)
        logger.success(f"Network saved to GeoParquet: {output_path}")
        logger.debug(f"File size: {output_path.stat().st_size / 1024:.2f} KB")
    elif format == "shapefile":
        import warnings

        logger.info(f"Saving river network to Shapefile: {output_path}")
        # Suppress shapefile field name truncation warnings (10 char limit is expected)
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", message=".*Column names longer than 10 characters.*")
            warnings.filterwarnings("ignore", message=".*Normalized/laundered field name.*")
            network.to_file(output_path)
        logger.success(f"Network saved to Shapefile: {output_path}")
        logger.debug(f"File size: {output_path.stat().st_size / 1024:.2f} KB")
    else:
        raise ValueError(f"Unsupported format: {format}. Use 'parquet' or 'shapefile'")

Load processed river network from GeoParquet file.

Loads river network data from GeoParquet format only (the recommended format). For loading shapefiles, use geopandas.read_file() directly.

Parameters:

Name	Type	Description	Default
`network_path`	`str \| Path`	Path to network GeoParquet file (.parquet)	required

Returns:

Type	Description
`GeoDataFrame`	GeoDataFrame with river network

Raises:

Type	Description
`FileNotFoundError`	If network file does not exist
`ValueError`	If file is not a .parquet file

Examples:

>>> from mobidic.preprocessing.io import load_network
>>> network = load_network("river_network.parquet")

Source code in mobidic/preprocessing/io.py

def load_network(network_path: str | Path) -> gpd.GeoDataFrame:
    """Load processed river network from GeoParquet file.

    Loads river network data from GeoParquet format only (the recommended format).
    For loading shapefiles, use geopandas.read_file() directly.

    Args:
        network_path: Path to network GeoParquet file (.parquet)

    Returns:
        GeoDataFrame with river network

    Raises:
        FileNotFoundError: If network file does not exist
        ValueError: If file is not a .parquet file

    Examples:
        >>> from mobidic.preprocessing.io import load_network
        >>> network = load_network("river_network.parquet")
    """
    network_path = Path(network_path)

    if not network_path.exists():
        raise FileNotFoundError(f"Network file not found: {network_path}")

    if network_path.suffix != ".parquet":
        raise ValueError(
            f"Only .parquet files supported. Got: {network_path.suffix}. "
            "Use geopandas.read_file() to load other formats."
        )

    logger.info(f"Loading river network from GeoParquet: {network_path}")
    network = gpd.read_parquet(network_path)
    logger.success(f"River network loaded: {len(network)} reaches")

    return network

Save reservoirs data to file.

Saves reservoir data to GeoParquet format, including all reservoir metadata, stage-storage curves, regulation curves, and spatial information.

Parameters:

Name	Type	Description	Default
`reservoirs`	`Reservoirs`	Reservoirs object containing reservoir data	required
`output_path`	`str \| Path`	Path to output file	required
`format`	`str`	Output format (only ‘parquet’ is supported)	`'parquet'`

Raises:

Type	Description
`ValueError`	If format is not supported

Examples:

>>> from mobidic import process_reservoirs
>>> reservoirs = process_reservoirs(...)
>>> from mobidic.preprocessing.io import save_reservoirs
>>> save_reservoirs(reservoirs, "reservoirs.parquet")

Source code in mobidic/preprocessing/io.py

def save_reservoirs(reservoirs: "Reservoirs", output_path: str | Path, format: str = "parquet") -> None:
    """Save reservoirs data to file.

    Saves reservoir data to GeoParquet format, including all reservoir metadata,
    stage-storage curves, regulation curves, and spatial information.

    Args:
        reservoirs: Reservoirs object containing reservoir data
        output_path: Path to output file
        format: Output format (only 'parquet' is supported)

    Raises:
        ValueError: If format is not supported

    Examples:
        >>> from mobidic import process_reservoirs
        >>> reservoirs = process_reservoirs(...)
        >>> from mobidic.preprocessing.io import save_reservoirs
        >>> save_reservoirs(reservoirs, "reservoirs.parquet")
    """
    from mobidic.preprocessing.reservoirs import Reservoirs

    if not isinstance(reservoirs, Reservoirs):
        raise TypeError(f"Expected Reservoirs object, got {type(reservoirs)}")

    output_path = Path(output_path)
    output_path.parent.mkdir(parents=True, exist_ok=True)

    if format.lower() != "parquet":
        raise ValueError(f"Unsupported format: {format}. Only 'parquet' is supported.")

    logger.info(f"Saving reservoirs to GeoParquet: {output_path}")
    reservoirs.save(output_path, format=format)
    logger.debug(f"File size: {output_path.stat().st_size / 1024:.2f} KB")

Load reservoirs data from file.

Loads reservoir data from GeoParquet format, including all reservoir metadata, stage-storage curves, regulation curves, and spatial information.

Parameters:

Name	Type	Description	Default
`input_path`	`str \| Path`	Path to reservoirs Parquet file	required

Returns:

Type	Description
`Reservoirs`	Reservoirs object with loaded data

Raises:

Type	Description
`FileNotFoundError`	If reservoir file does not exist
`ValueError`	If file is not a .parquet file

Examples:

>>> from mobidic.preprocessing.io import load_reservoirs
>>> reservoirs = load_reservoirs("reservoirs.parquet")

Source code in mobidic/preprocessing/io.py

def load_reservoirs(input_path: str | Path) -> "Reservoirs":
    """Load reservoirs data from file.

    Loads reservoir data from GeoParquet format, including all reservoir metadata,
    stage-storage curves, regulation curves, and spatial information.

    Args:
        input_path: Path to reservoirs Parquet file

    Returns:
        Reservoirs object with loaded data

    Raises:
        FileNotFoundError: If reservoir file does not exist
        ValueError: If file is not a .parquet file

    Examples:
        >>> from mobidic.preprocessing.io import load_reservoirs
        >>> reservoirs = load_reservoirs("reservoirs.parquet")
    """
    from mobidic.preprocessing.reservoirs import Reservoirs

    input_path = Path(input_path)

    if not input_path.exists():
        raise FileNotFoundError(f"Reservoirs file not found: {input_path}")

    if input_path.suffix != ".parquet":
        raise ValueError(
            f"Only .parquet files supported. Got: {input_path.suffix}. "
            "Use geopandas.read_parquet() to load other formats."
        )

    logger.info(f"Loading reservoirs from GeoParquet: {input_path}")
    reservoirs = Reservoirs.load(input_path)

    return reservoirs

Usage examples

Example 1: complete preprocessing workflow

from mobidic import load_config, run_preprocessing, GISData

# Load configuration
config = load_config("config.yaml")

# Run preprocessing
gisdata = run_preprocessing(config)

# Save preprocessed data (including reservoirs if configured)
gisdata.save(
    gisdata_path="output/gisdata.nc",
    network_path="output/network.parquet",
    reservoirs_path="output/reservoirs.parquet"  # Optional
)

# Later, load preprocessed data
loaded_gisdata = GISData.load(
    gisdata_path="output/gisdata.nc",
    network_path="output/network.parquet",
    reservoirs_path="output/reservoirs.parquet"  # Optional
)

# Access components
print(f"Grid variables: {list(loaded_gisdata.grids.keys())}")
print(f"Network reaches: {len(loaded_gisdata.network)}")
print(f"Grid shape: {loaded_gisdata.metadata['shape']}")
if loaded_gisdata.reservoirs:
    print(f"Reservoirs: {len(loaded_gisdata.reservoirs)}")

Example 2: working with GISData

from mobidic import GISData
import numpy as np

# Create GISData object
gisdata = GISData()

# Add grid data
gisdata.grids['dtm'] = np.random.rand(100, 100)
gisdata.grids['flow_dir'] = np.random.randint(1, 9, (100, 100))
gisdata.grids['ks'] = np.random.rand(100, 100) * 10

# Add network
gisdata.network = network_gdf  # GeoDataFrame from process_river_network()

# Add metadata
gisdata.metadata = {
    'shape': (100, 100),
    'resolution': 100.0,
    'crs': 'EPSG:32632',
    'transform': affine_transform,
}

# Save
gisdata.save("output/gisdata.nc", "output/network.parquet")

Example 3: saving network only

from mobidic import process_river_network, save_network, load_network

# Process network
network = process_river_network(
    shapefile_path="data/river_network.shp",
    join_single_tributaries=True,
)

# Save as GeoParquet (recommended)
save_network(network, "output/network.parquet", format="parquet")

# Or save as Shapefile
save_network(network, "output/network.shp", format="shapefile")

# Load network later
loaded_network = load_network("output/network.parquet")

Example 4: working with reservoirs

from mobidic import process_reservoirs, save_reservoirs, load_reservoirs, GISData

# Process reservoir data from shapefiles and CSVs
reservoirs = process_reservoirs(
    res_shape_path="reservoirs/reservoirs.shp",
    stage_storage_path="reservoirs/stage_storage.csv",
    regulation_curves_path="reservoirs/regulation_curves.csv",
    regulation_schedule_path="reservoirs/regulation_schedule.csv",
    initial_volumes_path="reservoirs/initial_volumes.csv",  # Optional
    grid_transform=gisdata.metadata['transform'],
    grid_shape=gisdata.metadata['shape'],
    network=gisdata.network,
)

# Save reservoirs to GeoParquet
save_reservoirs(reservoirs, "output/reservoirs.parquet")

# Load reservoirs later
loaded_reservoirs = load_reservoirs("output/reservoirs.parquet")

# Access reservoir data
for reservoir in loaded_reservoirs:
    print(f"Reservoir {reservoir.id}: {reservoir.name}")
    print(f"  Basin pixels: {len(reservoir.basin_pixels)}")
    print(f"  Inlet reaches: {reservoir.inlet_reaches}")
    print(f"  Outlet reach: {reservoir.outlet_reach}")
    print(f"  Initial volume: {reservoir.initial_volume} m³")
    print(f"  Stage-storage curve: {len(reservoir.stage_storage_curve)} points")

# Include in GISData
gisdata.reservoirs = reservoirs
gisdata.save(
    gisdata_path="output/gisdata.nc",
    network_path="output/network.parquet",
    reservoirs_path="output/reservoirs.parquet"
)

File formats

NetCDF for grid data

Grid data is saved in NetCDF4 format with:

Structure: - Each grid variable is a 2D data variable (y, x dimensions) - Coordinate variables for x and y - Comprehensive metadata (CRS, transform, resolution)

Compression: - zlib compression (level 4 by default) - Chunking optimized for spatial access patterns

Advantages: - Self-describing format with embedded metadata - Efficient compression for large grids - CF-compliant for interoperability - Supports NaN for nodata values

GeoParquet for network data

River network data is saved in GeoParquet format with:

Advantages: - Very fast read/write performance - Excellent compression ratios - Preserves all attribute types (int, float, list, etc.) - Native support for complex geometries - Column-oriented storage for efficient queries

Requirements: - Requires pyarrow package: pip install pyarrow

Fallback: - If pyarrow not available, can use Shapefile format - Shapefile has limitations (attribute names, data types)

GeoParquet for reservoir data

Reservoir data is saved in GeoParquet format with:

Structure: - One row per reservoir with all associated data - Polygon geometry for reservoir basin - Stage-storage curve as nested DataFrame - Regulation curves and schedules as dictionaries - Basin pixels as list of linear indices - Inlet/outlet reach references

Advantages: - Efficient storage of nested data structures - Preserves Python data types (lists, dicts, DataFrames) - Fast read/write performance - Column-oriented storage for selective loading

Special handling: - Dictionary keys are zero-padded strings (“000”, “001”, etc.) for Parquet compatibility - DataFrames are serialized to Parquet bytes and stored as binary columns

Data consistency

The GISData class ensures consistency between components:

Grid validation

All grids must have the same shape:

gisdata = GISData()
gisdata.grids['dtm'] = np.zeros((100, 100))
gisdata.grids['ks'] = np.zeros((100, 150))  # Different shape - will raise error on save

Metadata requirements

Required metadata fields: - shape: Tuple (nrows, ncols) - resolution: Float or tuple (x_res, y_res) - crs: String (e.g., “EPSG:32632”) - transform: Affine transform object

Network validation

The network must be a GeoDataFrame with specific required columns: - mobidic_id: Integer reach identifiers - geometry: LineString geometries - upstream_1, upstream_2, downstream: Topology references - strahler_order, calc_order: Ordering information - length_m, width_m: Geometric parameters

Spatial reference handling

CRS representation

The CRS is stored in multiple places:

NetCDF grids: As global attribute crs (WKT or EPSG code)
GeoDataFrame: Native CRS property
Metadata dict: As string for convenience

These should all be consistent and are validated on save/load.

Affine transform

The affine transform maps pixel coordinates to geographic coordinates:

from affine import Affine

# Example: 100m resolution, origin at (600000, 4800000)
transform = Affine(100.0, 0.0, 600000.0,
                   0.0, -100.0, 4900000.0)

# Transform pixel (row, col) to (x, y)
x, y = transform * (col, row)

The transform is stored in the metadata dict and embedded in NetCDF grid variables.

Integration with preprocessing

The preprocessing workflow automatically creates and populates GISData:

from mobidic import load_config, run_preprocessing

config = load_config("config.yaml")

# This function internally:
# 1. Creates GISData object
# 2. Reads all raster files specified in config
# 3. Processes river network
# 4. Computes hillslope-reach mapping
# 5. Processes reservoirs (if configured)
# 6. Populates grids, network, reservoirs, and metadata
gisdata = run_preprocessing(config)

# Save for later use
gisdata.save(config.paths.gisdata, config.paths.network)

Data I/O and GISData container

Overview

Classes

GISData container

__init__(grids, metadata, network, hillslope_reach_map, config, reservoirs=None)

load(gisdata_path, network_path, reservoirs_path=None) classmethod

save(gisdata_path, network_path, reservoirs_path=None)

Functions

Save and load functions

Usage examples

Example 1: complete preprocessing workflow

Example 2: working with GISData

Example 3: saving network only

Example 4: working with reservoirs

File formats

NetCDF for grid data

GeoParquet for network data

GeoParquet for reservoir data

Data consistency

Grid validation

Metadata requirements

Network validation

Spatial reference handling

CRS representation

Affine transform

Integration with preprocessing

`init(grids, metadata, network, hillslope_reach_map, config, reservoirs=None)`

`load(gisdata_path, network_path, reservoirs_path=None)` `classmethod`

`save(gisdata_path, network_path, reservoirs_path=None)`