IO

DHyDAMO IO API

common

ExtendedDataFrame

Bases: DataFrame

Source code in hydrolib/dhydamo/io/common.py

class ExtendedDataFrame(pd.DataFrame):
    _metadata = ["required_columns"] + pd.DataFrame._metadata

    def __init__(self, required_columns=None, *args, **kwargs):
        super().__init__(*args, **kwargs)

        if required_columns is None:
            required_columns = []

        self.required_columns = (
            required_columns[:]
            if isinstance(required_columns, list)
            else [required_columns]
        )

    def delete_all(self):
        """
        Empty the dataframe
        """
        if not self.empty:
            self.iloc[:, 0] = np.nan
            # inplace: returning a copy of an ExtendedDataFrame gives a plain DataFrame
            self.dropna(inplace=True)  # noqa: PD002

    def set_data(self, df, index_col):
        if not self.empty:
            self.delete_all()

        # Copy content
        for col, values in df.items():
            self[col] = values.to_numpy()

        if index_col is None:
            self.index = df.index
            self.index.name = df.index.name

        else:
            self.index = df[index_col]
            self.index.name = index_col

        # Check columns and types
        self._check_columns()

    def add_data(self, df):
        if not np.isin(df.columns, self.columns).all():
            raise KeyError(
                "The new df contains columns that are not present in the current df."
            )

        # Concatenate data
        current = pd.DataFrame(self.values, index=self.index, columns=self.columns)
        newdf = pd.concat([current, df], ignore_index=False, sort=False)

        # Empty current df
        self.delete_all()

        # Add values again
        self.set_data(newdf, index_col=self.index.name)

    def _check_columns(self):
        """
        Check presence of columns in geodataframe
        """
        present_columns = self.columns.tolist()
        for i, column in enumerate(self.required_columns):
            if column not in present_columns:
                raise KeyError(
                    'Column "{}" not found. Got {}, Expected at least {}'.format(
                        column,
                        ", ".join(present_columns),
                        ", ".join(self.required_columns),
                    )
                )

    def read_gpkg_layer(
        self,
        gpkg_path: str | Path,
        layer_name: str = None,
        column_mapping: dict = None,
        index_col: str = None,
    ):
        """Deprecated. Use ``HyDAMO.load_from_gpkg()`` for DAMO data loading."""
        warnings.warn(
            "read_gpkg_layer() is deprecated and will be removed in a future release. "
            "Use HyDAMO.load_from_gpkg() to load DAMO-formatted GeoPackages.",
            DeprecationWarning,
            stacklevel=2,
        )
        self._read_gpkg_layer(
            gpkg_path=gpkg_path,
            layer_name=layer_name,
            column_mapping=column_mapping,
            index_col=index_col,
        )

    def _read_gpkg_layer(
        self,
        gpkg_path: str | Path,
        layer_name: str = None,
        column_mapping: dict = None,
        index_col: str = None,
        **kwargs,
    ):
        """
        Read GML file to GeoDataFrame.

        This function has the option to group Points into LineStrings. To do so,
        specify the groupby column (which the set has in common) and the order_column,
        the column which indicates the order of the grouping.

        A mask file can be specified to clip the selection.

        Parameters
        ----------
        gml_path : str
            Path to the GML file
        layer_name: str
            Layer name of the desired layer in the package
        index_col : str
            Optional, column to be set as index
        """

        if not Path(gpkg_path).exists():
            raise OSError(f'File not found: "{gpkg_path}"')

        if isinstance(gpkg_path, Path):
            gpkg_path = str(gpkg_path)

        if layer_name.lower() not in map(str.lower, gpd.list_layers(gpkg_path).name.tolist()):
            raise ValueError(f'Layer "{layer_name}" does not exist in: "{gpkg_path}"')

        layer = gpd.read_file(gpkg_path, layer=layer_name, engine='pyogrio')
        columns = [col.lower() for col in layer.columns]
        layer.columns = columns
        if 'geometry' in layer.columns:
            layer = layer.drop("geometry", axis=1)

        df = layer
        if column_mapping is not None:
            df = df.rename(columns=column_mapping)

        # Add data to class GeoDataFrame
        self.set_data(df, index_col=index_col)

delete_all()

Empty the dataframe

Source code in hydrolib/dhydamo/io/common.py

def delete_all(self):
    """
    Empty the dataframe
    """
    if not self.empty:
        self.iloc[:, 0] = np.nan
        # inplace: returning a copy of an ExtendedDataFrame gives a plain DataFrame
        self.dropna(inplace=True)  # noqa: PD002

read_gpkg_layer(gpkg_path: str | Path, layer_name: str = None, column_mapping: dict = None, index_col: str = None)

Deprecated. Use HyDAMO.load_from_gpkg() for DAMO data loading.

Source code in hydrolib/dhydamo/io/common.py

def read_gpkg_layer(
    self,
    gpkg_path: str | Path,
    layer_name: str = None,
    column_mapping: dict = None,
    index_col: str = None,
):
    """Deprecated. Use ``HyDAMO.load_from_gpkg()`` for DAMO data loading."""
    warnings.warn(
        "read_gpkg_layer() is deprecated and will be removed in a future release. "
        "Use HyDAMO.load_from_gpkg() to load DAMO-formatted GeoPackages.",
        DeprecationWarning,
        stacklevel=2,
    )
    self._read_gpkg_layer(
        gpkg_path=gpkg_path,
        layer_name=layer_name,
        column_mapping=column_mapping,
        index_col=index_col,
    )

ExtendedGeoDataFrame

Bases: GeoDataFrame

Source code in hydrolib/dhydamo/io/common.py

class ExtendedGeoDataFrame(gpd.GeoDataFrame):
    # normal properties
    _metadata = ["required_columns", "geotype"] + gpd.GeoDataFrame._metadata

    def __init__(self, geotype, required_columns=None, logger=logging, *args, **kwargs):
        # Check type
        if required_columns is None:
            required_columns = []
        elif not isinstance(required_columns, list):
            required_columns = [required_columns]

        # Add required columns to column list
        if "columns" in kwargs.keys():
            kwargs["columns"] += required_columns
        else:
            kwargs["columns"] = required_columns

        super().__init__(*args, **kwargs)

        self.required_columns = required_columns[:]
        self.geotype = geotype

    def copy(self, deep=True):
        """
        Create a copy
        """

        index = self.index.tolist() if not deep else deepcopy(self.index.tolist())
        columns = self.columns.tolist() if not deep else deepcopy(self.columns.tolist())

        edf = ExtendedGeoDataFrame(
            geotype=self.geotype,
            required_columns=[],
            index=index,
            columns=columns,
        )
        edf.required_columns.extend(self.required_columns[:])
        edf.loc[:, :] = self.values if not deep else deepcopy(self.values)

        return edf

    def delete_all(self):
        """
        Empty the dataframe
        """
        if not self.empty:
            # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
            self.drop(self.index, inplace=True)  # noqa: PD002

    def read_shp(
        self,
        path: str | Path,
        index_col: str = None,
        column_mapping: dict = None,
        check_columns: bool = True,
        proj_crs=None,
        clip: Polygon | MultiPolygon = None,
        check_geotype: bool = True,
        id_col: str = "code",
        filter_cols: bool = False,
        filter_rows=None,
    ):
        """
        Import function, extended with type checks. Does not destroy reference to object.
        """
        # Read GeoDataFrame
        gdf = gpd.read_file(path)

        # Only keep required columns
        if filter_cols:
            logger.info("Filtering required column keys")
            gdf = gdf.drop(
                columns=gdf.columns[~gdf.columns.isin(self.required_columns)]
            )

        # filter out rows on key/value pairs if required
        if filter_rows is not None:
            logger.info("Filter rows using key value pairs")
            filtered = (gdf[list(filter_rows)] == pd.Series(filter_rows)).all(axis=1)
            gdf = gdf[filtered]

        # Drop features without geometry
        total_features = len(gdf)
        gdf = gdf.drop(gdf.index[gdf.geometry.isna()])  # temporary fix
        missing_features = len(gdf.index[gdf.geometry.isna()])
        logger.debug(
            f"{missing_features} out of {total_features} do not have a geometry"
        )

        # Rename columns:
        if column_mapping is not None:
            gdf = gdf.rename(columns=column_mapping)

        # Check number of entries
        if gdf.empty:
            raise OSError("Imported shapefile contains no rows.")

        # Add data to class GeoDataFrame
        self.set_data(
            gdf,
            index_col=index_col,
            check_columns=check_columns,
            check_geotype=check_geotype,
        )

        # Clip if extent is provided
        if clip is not None:
            self.clip(clip)

        # To re-project CRS system to projected CRS, first all empty geometries should be dropped
        if proj_crs is not None:
            self.check_projection(proj_crs)
        else:
            logger.debug("No projected CRS is given in ini-file")

    def set_data(self, gdf, index_col=None, check_columns=True, check_geotype=True):
        if not self.empty:
            self.delete_all()

        # Check columns
        if check_columns:
            self._check_columns(gdf)

        # Copy content
        for col, values in gdf.items():
            if str(values.dtype) == "geometry":
                # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
                self.set_geometry(values.to_numpy(), inplace=True)  # noqa: PD002
            else:
                self[col] = values.to_numpy()

        if index_col is None:
            self.index = gdf.index
            self.index.name = gdf.index.name

        else:
            self.index = gdf[index_col]
            self.index.name = index_col

        # Check geometry types
        if check_geotype:
            self._check_geotype()

    def _check_columns(self, gdf):
        """
        Check presence of columns in geodataframe
        """
        present_columns = gdf.columns.tolist()
        for column in self.required_columns:
            if column not in present_columns:
                # geopackages don't support very long names. Only give an error if the first 25 characters don't occur, otherwise lengthen the column name
                if column not in present_columns:
                    raise KeyError(
                        'Column "{}" not found. Got {}, Expected at least {}'.format(
                            column,
                            ", ".join(present_columns),
                            ", ".join(self.required_columns),
                        )
                    )

    def _check_geotype(self):
        """
        Check geometry type
        """
        if not all(isinstance(geo, self.geotype) for geo in self.geometry):
            raise TypeError(
                'Geometrytype "{}" required. The input shapefile has geometry type(s) {}.'.format(
                    re.findall("([A-Z].*)'", repr(self.geotype))[0],
                    self.geometry.type.unique().tolist(),
                )
            )

    def show_gpkg(self, gpkg_path: str | Path):
        if not Path(gpkg_path).exists():
            raise OSError(f'File not found: "{gpkg_path}"')

        if isinstance(gpkg_path, Path):
            gpkg_path = str(gpkg_path)

        layerlist = gpd.list_layers(gpkg_path).name.tolist()
        logger.info(
            "Content of gpkg-file %s, containing %d layers:",
            gpkg_path,
            len(layerlist),
        )
        logger.info(
            "\tINDEX\t|\tNAME                        \t|\tGEOM_TYPE      \t|\t NFEATURES\t|\t   NFIELDS"
        )
        for laynum, layer_name in enumerate(layerlist):
            layer = gpd.read_file(gpkg_path, layer=layer_name)
            if layer.empty:
                logger.warning('Layer "%s" is empty.', layer_name)
                continue

            nfields = len(layer.columns)
            nfeatures = layer.shape[0]

            if isinstance(layer, gpd.GeoDataFrame):
                geom_type = layer.geom_type.iloc[0]
                if geom_type is None:
                    geom_type = "None"                
            else:
                geom_type = "None"

            logger.info(
                "\t%5d\t|\t%30s\t|\t%s\t|\t%10d\t|\t%10d",
                laynum,
                layer_name,
                geom_type,
                nfeatures,
                nfields,
            )

    def read_gpkg_layer(
        self,
        gpkg_path: str | Path,
        layer_name: str,
        index_col: str = None,
        groupby_column: str = None,
        order_column: str = None,
        id_col: str = "code",
        column_mapping: dict = None,
        check_columns: bool = True,
        check_geotype: bool = True,
        clip: Polygon | MultiPolygon = None,
        cliptype: str = "select",
        check_3d: bool = True,
    ):
        """Deprecated. Use ``HyDAMO.load_from_gpkg()`` for DAMO data loading."""
        warnings.warn(
            "read_gpkg_layer() is deprecated and will be removed in a future release. "
            "Use HyDAMO.load_from_gpkg() to load DAMO-formatted GeoPackages.",
            DeprecationWarning,
            stacklevel=2,
        )
        self._read_gpkg_layer(
            gpkg_path=gpkg_path,
            layer_name=layer_name,
            index_col=index_col,
            groupby_column=groupby_column,
            order_column=order_column,
            id_col=id_col,
            column_mapping=column_mapping,
            check_columns=check_columns,
            check_geotype=check_geotype,
            clip=clip,
            cliptype=cliptype,
            check_3d=check_3d,
        )

    def _read_gpkg_layer(
        self,
        gpkg_path: str | Path,
        layer_name: str,
        index_col: str = None,
        groupby_column: str = None,
        order_column: str = None,
        id_col: str = "code",
        column_mapping: dict = None,
        check_columns: bool = True,
        check_geotype: bool = True,
        clip: Polygon | MultiPolygon = None,
        cliptype: str = "select",
        check_3d: bool = True,
    ):
        if not Path(gpkg_path).exists():
            raise OSError(f'File not found: "{gpkg_path}"')

        if isinstance(gpkg_path, Path):
            gpkg_path = str(gpkg_path)

        if layer_name.lower() not in map(str.lower, gpd.list_layers(gpkg_path).name.tolist()):
            raise ValueError(f'Layer "{layer_name}" does not exist in: "{gpkg_path}"')

        layer = gpd.read_file(gpkg_path, layer=layer_name, engine='pyogrio')
        columns = [col.lower() for col in layer.columns]
        layer.columns = columns

        if groupby_column is not None:
            fields, geometries = _points_to_lines(layer, groupby_column, order_column, check_3d)
        else:
            fields = layer
            geometries = layer.geometry

        # Create geodataframe
        gdf = gpd.GeoDataFrame(fields, columns=columns, geometry=geometries)

        if column_mapping is not None:
            gdf = gdf.rename(columns=column_mapping)

        if index_col is not None:
            _make_index_unique(gdf, index_col)

        # Add data to class GeoDataFrame
        self.set_data(
            gdf,
            index_col=index_col,
            check_columns=check_columns,
            check_geotype=check_geotype,
        )

        if clip is not None:
            self.clip(geometry=clip, cliptype=cliptype)

    def clip(self, geometry: Polygon | MultiPolygon, cliptype: str = "select", clip_and_drop: bool = True):
        """
        Clip geometry
        """
        if not isinstance(geometry, (Polygon, MultiPolygon)):
            raise TypeError("Expected geometry of type Polygon or MultiPolygon")

        # Clip if needed          
        if cliptype == "clip":
            pre_geomtypes = self.geom_type.unique().tolist()
            if "Polygon" in pre_geomtypes and "MultiPolygon" not in pre_geomtypes:
                pre_geomtypes.append("MultiPolygon")
            if "MultiPolygon" in pre_geomtypes and "Polygon" not in pre_geomtypes:
                pre_geomtypes.append("Polygon")
            gdf = gpd.clip(self, gpd.GeoDataFrame(geometry=[geometry], crs=self.crs))
            if clip_and_drop:
                # Reduce to allowed 
                gdf = gpd.GeoDataFrame(gdf, crs=gdf.crs)
                rem = gdf[~gdf.geom_type.isin(pre_geomtypes)]
                # Explode to deal with geometry collections
                rem = rem.explode(ignore_index=False, index_parts=False)
                rem = rem[rem.geom_type.isin(pre_geomtypes)]
                gdf = pd.concat([gdf[gdf.geom_type.isin(pre_geomtypes)], rem], ignore_index=False)
                # Merge duplicate indices to single multigeometry
                keep = []
                for idx, group in gdf.groupby(level=0, sort=False):
                    if len(group) == 1:
                        keep.append(group.iloc[0].copy())
                    elif len(group) > 1:
                        row = group.iloc[0].copy()
                        geoms = group.geometry.tolist()
                        geom = shapely.union_all(geoms)
                        if geom.geom_type not in pre_geomtypes:
                            logger.warning(
                                "Cannot deduplicate index '%s': clipping resulted in split geometries, using the first geometry.",
                                idx,
                            )
                            geom = geoms[0]
                        row.at["geometry"] = geom
                        keep.append(row)
                gdf = gpd.GeoDataFrame(keep, crs=gdf.crs)
        elif cliptype == "select":
            gdf = self.loc[self.intersects(geometry).to_numpy()]
        else:
            raise ValueError(f"Cliptype {cliptype} not recognized. Use 'clip' or 'select'.")
        if gdf.empty:
            raise ValueError("Found no features within extent geometry.")

        self.set_data(gdf)

    def check_projection(self, crs_out):
        """
        Check if reprojection is required
        """
        if crs_out != self.crs:
            # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
            self.to_crs(crs_out, inplace=True)  # noqa: PD002
        else:
            logger.info("OSM data has same projection as projected crs in ini-file")

    def branch_to_prof(
        self, offset=0.0, vertex_end=False, rename_col=None, prefix="", suffix=""
    ):
        """Create profiles on branches from branch data"""

        gdf_out = self.copy()

        # interpolate over feature geometries
        if vertex_end:
            chainage = self.length - offset
            p = self.interpolate(chainage)
        else:
            chainage = offset
            p = self.interpolate(chainage)
        gdf_out.geometry = p
        gdf_out["offset"] = chainage

        if rename_col is not None:
            try:
                gdf_out["branch_id"] = gdf_out[rename_col]
                gdf_out[rename_col] = [
                    f"{prefix}{g[1][rename_col]}{suffix}" for g in self.iterrows()
                ]
            except Exception:
                raise ValueError(f"Column rename with '{rename_col}' did not succeed.")

        return gdf_out

    def merge_columns(self, col1, col2, rename_col):
        """merge columns"""

        if col1 or col2 in self.columns.to_numpy():
            try:
                self[rename_col] = self[col1] + self[col2]
            except Exception:
                raise ValueError(
                    f"Merge of two profile columns'{col1}' and '{col2}' did not succeed."
                )

    def snap_to_branch(self, branches, snap_method, maxdist=5):
        """Snap the geometries to the branch"""
        spatial.find_nearest_branch(
            branches=branches, geometries=self, method=snap_method, maxdist=maxdist
        )

branch_to_prof(offset=0.0, vertex_end=False, rename_col=None, prefix='', suffix='')

Create profiles on branches from branch data

Source code in hydrolib/dhydamo/io/common.py

def branch_to_prof(
    self, offset=0.0, vertex_end=False, rename_col=None, prefix="", suffix=""
):
    """Create profiles on branches from branch data"""

    gdf_out = self.copy()

    # interpolate over feature geometries
    if vertex_end:
        chainage = self.length - offset
        p = self.interpolate(chainage)
    else:
        chainage = offset
        p = self.interpolate(chainage)
    gdf_out.geometry = p
    gdf_out["offset"] = chainage

    if rename_col is not None:
        try:
            gdf_out["branch_id"] = gdf_out[rename_col]
            gdf_out[rename_col] = [
                f"{prefix}{g[1][rename_col]}{suffix}" for g in self.iterrows()
            ]
        except Exception:
            raise ValueError(f"Column rename with '{rename_col}' did not succeed.")

    return gdf_out

check_projection(crs_out)

Check if reprojection is required

Source code in hydrolib/dhydamo/io/common.py

def check_projection(self, crs_out):
    """
    Check if reprojection is required
    """
    if crs_out != self.crs:
        # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
        self.to_crs(crs_out, inplace=True)  # noqa: PD002
    else:
        logger.info("OSM data has same projection as projected crs in ini-file")

clip(geometry: Polygon | MultiPolygon, cliptype: str = 'select', clip_and_drop: bool = True)

Clip geometry

Source code in hydrolib/dhydamo/io/common.py

def clip(self, geometry: Polygon | MultiPolygon, cliptype: str = "select", clip_and_drop: bool = True):
    """
    Clip geometry
    """
    if not isinstance(geometry, (Polygon, MultiPolygon)):
        raise TypeError("Expected geometry of type Polygon or MultiPolygon")

    # Clip if needed          
    if cliptype == "clip":
        pre_geomtypes = self.geom_type.unique().tolist()
        if "Polygon" in pre_geomtypes and "MultiPolygon" not in pre_geomtypes:
            pre_geomtypes.append("MultiPolygon")
        if "MultiPolygon" in pre_geomtypes and "Polygon" not in pre_geomtypes:
            pre_geomtypes.append("Polygon")
        gdf = gpd.clip(self, gpd.GeoDataFrame(geometry=[geometry], crs=self.crs))
        if clip_and_drop:
            # Reduce to allowed 
            gdf = gpd.GeoDataFrame(gdf, crs=gdf.crs)
            rem = gdf[~gdf.geom_type.isin(pre_geomtypes)]
            # Explode to deal with geometry collections
            rem = rem.explode(ignore_index=False, index_parts=False)
            rem = rem[rem.geom_type.isin(pre_geomtypes)]
            gdf = pd.concat([gdf[gdf.geom_type.isin(pre_geomtypes)], rem], ignore_index=False)
            # Merge duplicate indices to single multigeometry
            keep = []
            for idx, group in gdf.groupby(level=0, sort=False):
                if len(group) == 1:
                    keep.append(group.iloc[0].copy())
                elif len(group) > 1:
                    row = group.iloc[0].copy()
                    geoms = group.geometry.tolist()
                    geom = shapely.union_all(geoms)
                    if geom.geom_type not in pre_geomtypes:
                        logger.warning(
                            "Cannot deduplicate index '%s': clipping resulted in split geometries, using the first geometry.",
                            idx,
                        )
                        geom = geoms[0]
                    row.at["geometry"] = geom
                    keep.append(row)
            gdf = gpd.GeoDataFrame(keep, crs=gdf.crs)
    elif cliptype == "select":
        gdf = self.loc[self.intersects(geometry).to_numpy()]
    else:
        raise ValueError(f"Cliptype {cliptype} not recognized. Use 'clip' or 'select'.")
    if gdf.empty:
        raise ValueError("Found no features within extent geometry.")

    self.set_data(gdf)

copy(deep=True)

Create a copy

Source code in hydrolib/dhydamo/io/common.py

def copy(self, deep=True):
    """
    Create a copy
    """

    index = self.index.tolist() if not deep else deepcopy(self.index.tolist())
    columns = self.columns.tolist() if not deep else deepcopy(self.columns.tolist())

    edf = ExtendedGeoDataFrame(
        geotype=self.geotype,
        required_columns=[],
        index=index,
        columns=columns,
    )
    edf.required_columns.extend(self.required_columns[:])
    edf.loc[:, :] = self.values if not deep else deepcopy(self.values)

    return edf

delete_all()

Empty the dataframe

Source code in hydrolib/dhydamo/io/common.py

def delete_all(self):
    """
    Empty the dataframe
    """
    if not self.empty:
        # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
        self.drop(self.index, inplace=True)  # noqa: PD002

merge_columns(col1, col2, rename_col)

merge columns

Source code in hydrolib/dhydamo/io/common.py

def merge_columns(self, col1, col2, rename_col):
    """merge columns"""

    if col1 or col2 in self.columns.to_numpy():
        try:
            self[rename_col] = self[col1] + self[col2]
        except Exception:
            raise ValueError(
                f"Merge of two profile columns'{col1}' and '{col2}' did not succeed."
            )

read_gpkg_layer(gpkg_path: str | Path, layer_name: str, index_col: str = None, groupby_column: str = None, order_column: str = None, id_col: str = 'code', column_mapping: dict = None, check_columns: bool = True, check_geotype: bool = True, clip: Polygon | MultiPolygon = None, cliptype: str = 'select', check_3d: bool = True)

Deprecated. Use HyDAMO.load_from_gpkg() for DAMO data loading.

Source code in hydrolib/dhydamo/io/common.py

def read_gpkg_layer(
    self,
    gpkg_path: str | Path,
    layer_name: str,
    index_col: str = None,
    groupby_column: str = None,
    order_column: str = None,
    id_col: str = "code",
    column_mapping: dict = None,
    check_columns: bool = True,
    check_geotype: bool = True,
    clip: Polygon | MultiPolygon = None,
    cliptype: str = "select",
    check_3d: bool = True,
):
    """Deprecated. Use ``HyDAMO.load_from_gpkg()`` for DAMO data loading."""
    warnings.warn(
        "read_gpkg_layer() is deprecated and will be removed in a future release. "
        "Use HyDAMO.load_from_gpkg() to load DAMO-formatted GeoPackages.",
        DeprecationWarning,
        stacklevel=2,
    )
    self._read_gpkg_layer(
        gpkg_path=gpkg_path,
        layer_name=layer_name,
        index_col=index_col,
        groupby_column=groupby_column,
        order_column=order_column,
        id_col=id_col,
        column_mapping=column_mapping,
        check_columns=check_columns,
        check_geotype=check_geotype,
        clip=clip,
        cliptype=cliptype,
        check_3d=check_3d,
    )

read_shp(path: str | Path, index_col: str = None, column_mapping: dict = None, check_columns: bool = True, proj_crs=None, clip: Polygon | MultiPolygon = None, check_geotype: bool = True, id_col: str = 'code', filter_cols: bool = False, filter_rows=None)

Import function, extended with type checks. Does not destroy reference to object.

Source code in hydrolib/dhydamo/io/common.py

def read_shp(
    self,
    path: str | Path,
    index_col: str = None,
    column_mapping: dict = None,
    check_columns: bool = True,
    proj_crs=None,
    clip: Polygon | MultiPolygon = None,
    check_geotype: bool = True,
    id_col: str = "code",
    filter_cols: bool = False,
    filter_rows=None,
):
    """
    Import function, extended with type checks. Does not destroy reference to object.
    """
    # Read GeoDataFrame
    gdf = gpd.read_file(path)

    # Only keep required columns
    if filter_cols:
        logger.info("Filtering required column keys")
        gdf = gdf.drop(
            columns=gdf.columns[~gdf.columns.isin(self.required_columns)]
        )

    # filter out rows on key/value pairs if required
    if filter_rows is not None:
        logger.info("Filter rows using key value pairs")
        filtered = (gdf[list(filter_rows)] == pd.Series(filter_rows)).all(axis=1)
        gdf = gdf[filtered]

    # Drop features without geometry
    total_features = len(gdf)
    gdf = gdf.drop(gdf.index[gdf.geometry.isna()])  # temporary fix
    missing_features = len(gdf.index[gdf.geometry.isna()])
    logger.debug(
        f"{missing_features} out of {total_features} do not have a geometry"
    )

    # Rename columns:
    if column_mapping is not None:
        gdf = gdf.rename(columns=column_mapping)

    # Check number of entries
    if gdf.empty:
        raise OSError("Imported shapefile contains no rows.")

    # Add data to class GeoDataFrame
    self.set_data(
        gdf,
        index_col=index_col,
        check_columns=check_columns,
        check_geotype=check_geotype,
    )

    # Clip if extent is provided
    if clip is not None:
        self.clip(clip)

    # To re-project CRS system to projected CRS, first all empty geometries should be dropped
    if proj_crs is not None:
        self.check_projection(proj_crs)
    else:
        logger.debug("No projected CRS is given in ini-file")

snap_to_branch(branches, snap_method, maxdist=5)

Snap the geometries to the branch

Source code in hydrolib/dhydamo/io/common.py

def snap_to_branch(self, branches, snap_method, maxdist=5):
    """Snap the geometries to the branch"""
    spatial.find_nearest_branch(
        branches=branches, geometries=self, method=snap_method, maxdist=maxdist
    )

damo_converters

base

BaseDamoConverter

Bases: ABC

Abstract base class for versioned HyDAMO GeoPackage converters.

Subclasses declare hydamo_version and layer_specs, and optionally override postprocess for any column-level canonicalization required after all layers are loaded.

Source code in hydrolib/dhydamo/io/damo_converters/base.py

class BaseDamoConverter(ABC):
    """Abstract base class for versioned HyDAMO GeoPackage converters.

    Subclasses declare ``hydamo_version`` and ``layer_specs``, and optionally
    override ``postprocess`` for any column-level canonicalization required
    after all layers are loaded.
    """

    hydamo_version: str
    layer_specs: tuple[LayerSpec, ...]

    def load_into(
        self,
        hydamo,
        gpkg_path: str | Path,
        clip_layers: dict | None = None,
        check_3d: bool = True,
    ) -> None:
        """Read all declared layer specs from *gpkg_path* into *hydamo*.

        Parameters
        ----------
        hydamo : HyDAMO
            The ``HyDAMO`` instance whose attributes are populated.
        gpkg_path : str or Path
            Path to the source GeoPackage.
        clip_layers : dict, optional
            Per-layer spatial clip. Maps ``target_attr`` to a
            ``(geometry, cliptype)`` tuple. Layers absent from the dict are
            loaded without clipping.
        check_3d : bool
            Whether to require Z coordinates on profile points when grouping
            them into lines. Other layers are not checked for Z coordinates.
            Pass ``False`` to allow 2D profile points. Default is ``True``.
        """
        gpkg_path = Path(gpkg_path)
        available_layers = {
            layer_name.lower(): layer_name
            for layer_name in gpd.list_layers(gpkg_path).name.tolist()
        }

        for spec in self.layer_specs:
            actual_layer = available_layers.get(spec.source_layer.lower())
            if actual_layer is None:
                if spec.optional:
                    logger.info(
                        'Skipping optional HyDAMO layer "%s" for DAMO %s.',
                        spec.source_layer,
                        self.hydamo_version,
                    )
                    continue
                raise ValueError(
                    f'Missing required HyDAMO layer "{spec.source_layer}" for DAMO {self.hydamo_version}.'
                )

            target = getattr(hydamo, spec.target_attr)
            kwargs = {
                "gpkg_path": gpkg_path,
                "layer_name": actual_layer,
                "index_col": spec.index_col,
                "column_mapping": spec.column_mapping if spec.column_mapping else None,
                "check_geotype": spec.check_geotype,
                "check_3d": check_3d,
            }
            if spec.groupby_column is not None:
                kwargs["groupby_column"] = spec.groupby_column
            if spec.order_column is not None:
                kwargs["order_column"] = spec.order_column
            if clip_layers and spec.target_attr in clip_layers:
                kwargs["clip"], kwargs["cliptype"] = clip_layers[spec.target_attr]

            target._read_gpkg_layer(**kwargs)

        self.postprocess(hydamo)

    def postprocess(self, hydamo) -> None:
        """Hook for version-specific canonicalization after all layers are loaded.

        Override in subclasses to rename, copy, or derive columns on any
        ``HyDAMO`` attribute that requires post-load normalisation.
        """

load_into(hydamo, gpkg_path: str | Path, clip_layers: dict | None = None, check_3d: bool = True) -> None

Read all declared layer specs from gpkg_path into hydamo.

Parameters

hydamo : HyDAMO The HyDAMO instance whose attributes are populated. gpkg_path : str or Path Path to the source GeoPackage. clip_layers : dict, optional Per-layer spatial clip. Maps target_attr to a (geometry, cliptype) tuple. Layers absent from the dict are loaded without clipping. check_3d : bool Whether to require Z coordinates on profile points when grouping them into lines. Other layers are not checked for Z coordinates. Pass False to allow 2D profile points. Default is True.

Source code in hydrolib/dhydamo/io/damo_converters/base.py

def load_into(
    self,
    hydamo,
    gpkg_path: str | Path,
    clip_layers: dict | None = None,
    check_3d: bool = True,
) -> None:
    """Read all declared layer specs from *gpkg_path* into *hydamo*.

    Parameters
    ----------
    hydamo : HyDAMO
        The ``HyDAMO`` instance whose attributes are populated.
    gpkg_path : str or Path
        Path to the source GeoPackage.
    clip_layers : dict, optional
        Per-layer spatial clip. Maps ``target_attr`` to a
        ``(geometry, cliptype)`` tuple. Layers absent from the dict are
        loaded without clipping.
    check_3d : bool
        Whether to require Z coordinates on profile points when grouping
        them into lines. Other layers are not checked for Z coordinates.
        Pass ``False`` to allow 2D profile points. Default is ``True``.
    """
    gpkg_path = Path(gpkg_path)
    available_layers = {
        layer_name.lower(): layer_name
        for layer_name in gpd.list_layers(gpkg_path).name.tolist()
    }

    for spec in self.layer_specs:
        actual_layer = available_layers.get(spec.source_layer.lower())
        if actual_layer is None:
            if spec.optional:
                logger.info(
                    'Skipping optional HyDAMO layer "%s" for DAMO %s.',
                    spec.source_layer,
                    self.hydamo_version,
                )
                continue
            raise ValueError(
                f'Missing required HyDAMO layer "{spec.source_layer}" for DAMO {self.hydamo_version}.'
            )

        target = getattr(hydamo, spec.target_attr)
        kwargs = {
            "gpkg_path": gpkg_path,
            "layer_name": actual_layer,
            "index_col": spec.index_col,
            "column_mapping": spec.column_mapping if spec.column_mapping else None,
            "check_geotype": spec.check_geotype,
            "check_3d": check_3d,
        }
        if spec.groupby_column is not None:
            kwargs["groupby_column"] = spec.groupby_column
        if spec.order_column is not None:
            kwargs["order_column"] = spec.order_column
        if clip_layers and spec.target_attr in clip_layers:
            kwargs["clip"], kwargs["cliptype"] = clip_layers[spec.target_attr]

        target._read_gpkg_layer(**kwargs)

    self.postprocess(hydamo)

postprocess(hydamo) -> None

Hook for version-specific canonicalization after all layers are loaded.

Override in subclasses to rename, copy, or derive columns on any HyDAMO attribute that requires post-load normalisation.

Source code in hydrolib/dhydamo/io/damo_converters/base.py

def postprocess(self, hydamo) -> None:
    """Hook for version-specific canonicalization after all layers are loaded.

    Override in subclasses to rename, copy, or derive columns on any
    ``HyDAMO`` attribute that requires post-load normalisation.
    """

LayerSpec dataclass

Schema descriptor for a single HyDAMO GeoPackage layer.

Parameters

target_attr : str Name of the HyDAMO instance attribute to populate. source_layer : str GeoPackage layer name to read (case-insensitive match). index_col : str, optional Column to set as the DataFrame index after loading. groupby_column : str, optional Column used to group point rows into line geometries. order_column : str, optional Column that defines point order within each group. column_mapping : dict[str, str], optional {source_col: target_col} renames applied before writing to the attribute. check_geotype : bool Whether to validate that geometries match the expected type declared on the target ExtendedGeoDataFrame. Default is True. optional : bool When True, a missing GeoPackage layer is logged and skipped rather than raising an error. Default is True.

Source code in hydrolib/dhydamo/io/damo_converters/base.py

@dataclass(frozen=True)
class LayerSpec:
    """Schema descriptor for a single HyDAMO GeoPackage layer.

    Parameters
    ----------
    target_attr : str
        Name of the ``HyDAMO`` instance attribute to populate.
    source_layer : str
        GeoPackage layer name to read (case-insensitive match).
    index_col : str, optional
        Column to set as the DataFrame index after loading.
    groupby_column : str, optional
        Column used to group point rows into line geometries.
    order_column : str, optional
        Column that defines point order within each group.
    column_mapping : dict[str, str], optional
        ``{source_col: target_col}`` renames applied before writing to the
        attribute.
    check_geotype : bool
        Whether to validate that geometries match the expected type declared
        on the target ``ExtendedGeoDataFrame``. Default is ``True``.
    optional : bool
        When ``True``, a missing GeoPackage layer is logged and skipped
        rather than raising an error. Default is ``True``.
    """

    target_attr: str
    source_layer: str
    index_col: str | None = None
    groupby_column: str | None = None
    order_column: str | None = None
    column_mapping: dict[str, str] = field(default_factory=dict)
    check_geotype: bool = True
    optional: bool = True

replace_layer_specs(base_specs: tuple[LayerSpec, ...], *overrides: LayerSpec) -> tuple[LayerSpec, ...]

Merge overrides into base_specs, replacing entries by target_attr.

Specs present in base_specs are replaced in-place by the matching override if one exists. Overrides with no matching base spec are appended at the end.

Parameters

base_specs : tuple[LayerSpec, ...] The canonical spec tuple from the parent converter. *overrides : LayerSpec Replacement LayerSpec objects keyed by target_attr.

Returns

tuple[LayerSpec, ...] New tuple with overrides applied.

Source code in hydrolib/dhydamo/io/damo_converters/base.py

def replace_layer_specs(
    base_specs: tuple[LayerSpec, ...], *overrides: LayerSpec
) -> tuple[LayerSpec, ...]:
    """Merge *overrides* into *base_specs*, replacing entries by ``target_attr``.

    Specs present in *base_specs* are replaced in-place by the matching
    override if one exists. Overrides with no matching base spec are
    appended at the end.

    Parameters
    ----------
    base_specs : tuple[LayerSpec, ...]
        The canonical spec tuple from the parent converter.
    *overrides : LayerSpec
        Replacement ``LayerSpec`` objects keyed by ``target_attr``.

    Returns
    -------
    tuple[LayerSpec, ...]
        New tuple with overrides applied.
    """
    override_map = {spec.target_attr: spec for spec in overrides}
    merged_specs = [override_map.pop(spec.target_attr, spec) for spec in base_specs]
    merged_specs.extend(override_map.values())
    return tuple(merged_specs)

v22

Damo22Converter

Bases: BaseDamoConverter

DAMO v2.2 converter.

Defines the canonical layer-to-attribute mapping used as the baseline for all other version converters. Post-processing copies soortregelmiddel to soortafsluitmiddel on management_device to normalise the column name.

Source code in hydrolib/dhydamo/io/damo_converters/v22.py

class Damo22Converter(BaseDamoConverter):
    """DAMO v2.2 converter.

    Defines the canonical layer-to-attribute mapping used as the baseline for
    all other version converters. Post-processing copies soortregelmiddel to
    soortafsluitmiddel on management_device to normalise the column name.
    """

    hydamo_version = "2.2"
    layer_specs = (
        LayerSpec("branches", "HydroObject", index_col="code"),
        LayerSpec(
            "profile",
            "ProfielPunt",
            groupby_column="profiellijnid",
            order_column="codevolgnummer",
        ),
        LayerSpec("profile_roughness", "RuwheidProfiel"),
        LayerSpec("profile_line", "profiellijn"),
        LayerSpec("profile_group", "profielgroep"),
        LayerSpec("param_profile", "NormgeparamProfiel"),
        LayerSpec("param_profile_values", "NormgeparamProfielWaarde"),
        LayerSpec("weirs", "Stuw"),
        LayerSpec("opening", "Kunstwerkopening"),
        LayerSpec("closing_device", "Afsluitmiddel"),
        LayerSpec("management_device", "Regelmiddel"),
        LayerSpec("bridges", "Brug", index_col="code"),
        LayerSpec("culverts", "DuikerSifonHevel", index_col="code"),
        LayerSpec("pumpstations", "Gemaal", index_col="code"),
        LayerSpec("pumps", "Pomp", index_col="code"),
        LayerSpec("management", "Sturing", index_col="code"),
        LayerSpec(
            "boundary_conditions",
            "hydrologischerandvoorwaarde",
            index_col="code",
        ),
        LayerSpec(
            "catchments",
            "afvoergebiedaanvoergebied",
            index_col="code",
            check_geotype=False,
        ),
        LayerSpec("laterals", "lateraleknoop"),
    )

    def postprocess(self, hydamo) -> None:
        device = hydamo.management_device
        if (
            not device.empty
            and "soortregelmiddel" in device.columns
            and "soortafsluitmiddel" not in device.columns
        ):
            device["soortafsluitmiddel"] = device["soortregelmiddel"]
        if not hydamo.profile.empty and "profiellijnid" in hydamo.profile.columns:
            # After groupby, each row represents a profile line. Replace the
            # per-point code with the profile line's own ID.
            hydamo.profile["code"] = hydamo.profile["profiellijnid"]
        if not hydamo.catchments.empty and not hydamo.laterals.empty:
            # boundary_node is not a DAMO schema column, but it is required by
            # drrreader.py and drrwriter.py for every RR workflow. It resolves
            # the catchment→lateral FK (lateraleknoopid → laterals.globalid)
            # to the lateral's code, which is used as the RR boundary node ID.
            catchment_lateral = RelationsRR.CATCHMENT_LATERAL
            mapping = hydamo.laterals.set_index(catchment_lateral.parent_key)["code"]
            hydamo.catchments["boundary_node"] = hydamo.catchments[
                catchment_lateral.foreign_key
            ].map(mapping)

v23

Damo23Converter

Bases: Damo22Converter

DAMO v2.3 converter.

Schema is identical to v2.2. This class exists as a version sentinel so that files tagged 2.3 load correctly without raising an unsupported-version error. If v2.3 introduces schema differences in the future, override layer_specs here.

Source code in hydrolib/dhydamo/io/damo_converters/v23.py

class Damo23Converter(Damo22Converter):
    """DAMO v2.3 converter.

    Schema is identical to v2.2. This class exists as a version sentinel so that
    files tagged 2.3 load correctly without raising an unsupported-version error.
    If v2.3 introduces schema differences in the future, override layer_specs here.
    """

    hydamo_version = "2.3"

v24

Damo24Converter

Bases: Damo22Converter

DAMO v2.4 converter.

Schema is identical to v2.2. This class exists as a version sentinel so that files tagged 2.4 load correctly without raising an unsupported-version error. If v2.4 introduces schema differences in the future, override layer_specs here.

Source code in hydrolib/dhydamo/io/damo_converters/v24.py

class Damo24Converter(Damo22Converter):
    """DAMO v2.4 converter.

    Schema is identical to v2.2. This class exists as a version sentinel so that
    files tagged 2.4 load correctly without raising an unsupported-version error.
    If v2.4 introduces schema differences in the future, override layer_specs here.
    """

    hydamo_version = "2.4"

v25

Damo25Converter

Bases: Damo22Converter

DAMO v2.5 converter.

Schema changes relative to v2.2: - HydroObject: categorieoppervlaktewater → categorieoppwaterlichaam - Kunstwerkopening: afsluitmiddelid → regelmiddelid - Afsluitmiddel replaces Regelmiddel as the management device layer; typeafsluitmiddel → soortafsluitmiddel, typeregelbaarheid → soortregelbaarheid - Sturing: afsluitmiddelid → regelmiddelid, typestuurvariabele → stuurvariabele

Source code in hydrolib/dhydamo/io/damo_converters/v25.py

class Damo25Converter(Damo22Converter):
    """DAMO v2.5 converter.

    Schema changes relative to v2.2:
    - HydroObject: categorieoppervlaktewater → categorieoppwaterlichaam
    - Kunstwerkopening: afsluitmiddelid → regelmiddelid
    - Afsluitmiddel replaces Regelmiddel as the management device layer;
      typeafsluitmiddel → soortafsluitmiddel, typeregelbaarheid → soortregelbaarheid
    - Sturing: afsluitmiddelid → regelmiddelid, typestuurvariabele → stuurvariabele
    """

    hydamo_version = "2.5"
    layer_specs = replace_layer_specs(
        Damo22Converter.layer_specs,
        LayerSpec(
            "branches",
            "HydroObject",
            index_col="code",
            column_mapping={"categorieoppervlaktewater": "categorieoppwaterlichaam"},
        ),
        LayerSpec(
            "opening",
            "Kunstwerkopening",
            column_mapping={"afsluitmiddelid": "regelmiddelid"},
        ),
        # Override the inherited Afsluitmiddel -> closing_device mapping:
        # in 2.5, Afsluitmiddel populates the internal management_device below.
        # The internal closing_device attribute is retained, but Regelmiddel
        # is absent from 2.5 input, so this optional spec is skipped.
        # Removing this override alone would load Afsluitmiddel into both tables.
        LayerSpec("closing_device", "Regelmiddel"),
        LayerSpec(
            "management_device",
            "Afsluitmiddel",
            column_mapping={
                "typeafsluitmiddel": "soortafsluitmiddel",
                "typeregelbaarheid": "soortregelbaarheid",
            },
        ),
        LayerSpec(
            "management",
            "Sturing",
            index_col="code",
            column_mapping={
                "afsluitmiddelid": "regelmiddelid",
                "typestuurvariabele": "stuurvariabele",
            },
        ),
    )

dimrwriter

DIMRWriter

Temporary files for DIMR-configuration - to be superseded by Hydrolib-core functionality.

Source code in hydrolib/dhydamo/io/dimrwriter.py

class DIMRWriter:
    """Temporary files for DIMR-configuration - to be superseded by Hydrolib-core functionality."""

    @validate_arguments
    def __init__(
        self, dimr_path: str = None, output_path: str | Path = None
    ) -> None:
        """Class initialization for the DIMR writer. This is a temporary fix to include writers for RR and RTC, that are not yet fully included in Hydrolib-core, or have bugs. Eventually, this file should be redundant as this functionality is fully included in Hydrolib-core

        Args:
            dimr_path (Union[str, Path], optional): _description_. Defaults to None.
            output_path (Union[str, Path], optional): _description_. Defaults to None.
        """
        if dimr_path is None:
            self.run_dimr = r"C:\Program Files\Deltares\D-HYDRO Suite 2022.04 1D2D\plugins\DeltaShell.Dimr\kernels\x64\dimr\scripts\run_dimr.bat"
        else:
            self.run_dimr = dimr_path

        if output_path is None:
            self.output_path = Path(os.path.abspath("."))
        else:
            self.output_path = output_path

        self.template_dir = Path(os.path.abspath("."))

    @validate_arguments
    def write_runbat(self, debuglevel=6, runlog=None) -> None:
        """Generates a run.bat to run the model in DIMR. The path to the executable is provided by the user in the class initialization, or set to the default D-Hydro installation folder."""

        with open(os.path.join(self.output_path, "run.bat"), "w") as f:
            f.write("@ echo off\n")
            f.write("set OMP_NUM_THREADS=2\n")
            if runlog is None:
                f.write('call "' + str(self.run_dimr) + '" -d ' + str(debuglevel) + '\n')
            else:
                f.write('call "' + str(self.run_dimr) + '" -d ' + str(debuglevel) + ' > ' + str(runlog) + '\n')

    @validate_arguments
    def add_crs(self, netcdf_path: str | Path = None) -> None:
        """Reads the Netcdf file and addes the required attributes for a valid CRS."""        
        if netcdf_path is None:
            netfile = list((self.output_path / 'dflowfm').glob('*.nc'))[0]
            if not netfile.exists():
                raise FileNotFoundError(f"Netcdf file not found in {self.output_path / 'dflowfm'}. Provide the correct path via the netcdf_path argument.")
        else:
            if netcdf_path.is_dir():
                netfile = list(netcdf_path.glob('*.nc'))[0]
            elif netcdf_path.is_file():
                netfile = netcdf_path  
            else:        
                raise FileNotFoundError(f"Netcdf file not found at {netcdf_path}.")   
        netf = nc.Dataset(netfile, 'r+')
        netf.Conventions =  'CF-1.8 UGRID-1.0 Deltares-0.10'
        proj = netf.createVariable('projected_coordinate_system','i4')
        proj.epsg = 28992
        proj.grid_mapping_name = "Amersfoort / RD New" #"Rijksdriehoeksstelsel"
        proj.longitude_of_prime_meridian = 0.0#; // double
        proj.semi_major_axis = 6377397.155#; // double
        proj.semi_minor_axis = 6356078.962818189#; // double
        proj.inverse_flattening = 299.1528128#; // double
        proj.proj4_params = "+prdoj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
        proj.EPSG_code = "EPSG:28992"
        proj.value = "value is equal to EPSG code"
        proj.wkt = "PROJCS[\"Amersfoort / RD New\",\n    GEOGCS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,\n"             
        netf.close()   

    @staticmethod
    def _add_unique_coupler_item(
        coupler: ET.Element,
        source_name: str,
        target_name: str,
        seen_pairs: set[tuple[str, str]],
        gn_brackets: str,
    ) -> bool:
        if source_name is None or target_name is None:
            return False

        pair = (source_name, target_name)
        if pair in seen_pairs:
            return False
        seen_pairs.add(pair)

        item = ET.SubElement(coupler, gn_brackets + "item")
        item.text = ""
        item.tail = "\n"

        source = ET.SubElement(item, gn_brackets + "sourceName")
        source.text = source_name
        source.tail = "\n"

        target = ET.SubElement(item, gn_brackets + "targetName")
        target.text = target_name
        target.tail = "\n"
        return True

    @staticmethod
    def _get_complex_coupler_blocks(
        rtc_model: DRTCModel, name_prefix: str
    ) -> list[ET.Element]:
        """Return the <coupler> blocks from the complex controllers' dimr_config whose name starts with name_prefix."""
        if rtc_model.complex_controllers is None:
            return []

        complex_config = rtc_model.complex_controllers["dimr_config"][0]
        blocks = []
        for block in complex_config:
            block_tag = block.tag.split("}", 1)[-1] if block.tag.startswith("{") else block.tag
            if block_tag == "coupler" and block.attrib["name"].lower().startswith(name_prefix):
                blocks.append(block)
        return blocks

    #@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def write_dimrconfig(
        self, fm: FMModel, rr_model: DRRModel = None, rtc_model: DRTCModel = None
    ) -> None:
        """Construct the contents of a DIMR-config file, using elements from FM, RTC and RR models.

        Args:
            fm (instance of class FM-model): Hydrolib-core model object. When only fm is set, only a 1d2d model is generated.
            rr_model (instance of a DRRModel object, optional): Class containig all information from the RR-model. Defaults to None.
            rtc_model (instance of a DRTCmodel, optional): Class containig all information from the RTC-model. Defaults to None.
        """
        # initiate flags for RR and RTC blocks
        RR = False
        RTC = False
        if rr_model is not None:
            RR = True
        if rtc_model is not None:
            RTC = True

        generalname = "http://schemas.deltares.nl/dimr"  # Config
        xsi_name = "http://www.w3.org/2001/XMLSchema-instance"
        gn_brackets = "{" + generalname + "}"

        # registering namespaces
        ET.register_namespace("", generalname)
        ET.register_namespace("xsi", xsi_name)

        # Parsing xml file
        configfile = ET.parse(self.template_dir / "dimr_config.xml")

        myroot = configfile.getroot()
        myroot[1][1].text = fm.filepath.parts[-2]
        myroot[1][2].text = fm.filepath.name

        # control blocks
        control = ET.Element(gn_brackets + "control")
        control.text = ""
        control.tail = "\n"
        myroot.insert(1, control)

        if RR or RTC:
            parallel = ET.SubElement(control, gn_brackets + "parallel")
            parallel.tail = "\n"
            parallel.text = ""

        if RTC:
            # a flow_to_rtc coupler is only needed for PID/interval controllers,
            # or complex controllers whose dimr_config defines one themselves
            complex_flow_to_rtc_blocks = self._get_complex_coupler_blocks(
                rtc_model, "flow_to_rtc"
            )
            needs_flow_to_rtc_coupler = (
                bool(rtc_model.pid_controllers)
                or bool(rtc_model.interval_controllers)
                or bool(complex_flow_to_rtc_blocks)
            )

            startgrouprtc = ET.SubElement(parallel, gn_brackets + "startGroup")
            startgrouprtc.text = ""
            startgrouprtc.tail = "\n"

            timertc = ET.SubElement(startgrouprtc, gn_brackets + "time")
            timertc.text = f"{0} {rtc_model.time_settings['step']} {fm.time.tstop}"
            timertc.tail = "\n"
            if needs_flow_to_rtc_coupler:
                couplerrtcrr = ET.SubElement(startgrouprtc, gn_brackets + "coupler")
                couplerrtcrr.attrib = {"name": "flow_to_rtc"}
                couplerrtcrr.tail = "\n"

            startrtc = ET.SubElement(startgrouprtc, gn_brackets + "start")
            startrtc.attrib = {"name": "Real_Time_Control"}
            startrtc.tail = "\n"
            couplerrtc = ET.SubElement(startgrouprtc, gn_brackets + "coupler")
            couplerrtc.attrib = {"name": "rtc_to_flow"}
            couplerrtc.tail = "\n"

        if RR:
            startgrouprr = ET.SubElement(parallel, gn_brackets + "startGroup")
            startgrouprr.text = ""
            startgrouprr.tail = "\n"

            timerr = ET.SubElement(startgrouprr, gn_brackets + "time")
            timerr.text = (
                f"{0} {rr_model.d3b_parameters['Timestepsize']} {fm.time.tstop}"
            )
            timerr.tail = "\n"
            couplerrrf = ET.SubElement(startgrouprr, gn_brackets + "coupler")
            couplerrrf.attrib = {"name": "flow_to_rr"}
            couplerrrf.tail = "\n"
            startrr = ET.SubElement(startgrouprr, gn_brackets + "start")
            startrr.attrib = {"name": "RR"}
            startrr.tail = "\n"
            couplerfrr = ET.SubElement(startgrouprr, gn_brackets + "coupler")
            couplerfrr.attrib = {"name": "rr_to_flow"}
            couplerfrr.tail = "\n"
        if RR or RTC:
            startfm = ET.SubElement(parallel, gn_brackets + "start")
        else:
            startfm = ET.SubElement(control, gn_brackets + "start")
        startfm.attrib = {"name": "DFM"}
        startfm.tail = "\n"

        # component blocks
        if RTC:
            comprtc = ET.Element(gn_brackets + "component")
            comprtc.attrib = {"name": "Real_Time_Control"}
            comprtc.tail = "\n"
            librtc = ET.SubElement(comprtc, gn_brackets + "library")
            librtc.text = "FBCTools_BMI"
            librtc.tail = "\n"
            wdrtc = ET.SubElement(comprtc, gn_brackets + "workingDir")
            wdrtc.text = "rtc"
            wdrtc.tail = "\n"
            frtc = ET.SubElement(comprtc, gn_brackets + "inputFile")
            frtc.text = "."
            frtc.tail = "\n"
            myroot.insert(2, comprtc)

        if RR:
            comprr = ET.Element(gn_brackets + "component")
            comprr.attrib = {"name": "RR"}
            comprr.tail = "\n"
            librr = ET.SubElement(comprr, gn_brackets + "library")
            librr.text = "rr_dll"
            librr.tail = "\n"
            wdrr = ET.SubElement(comprr, gn_brackets + "workingDir")
            wdrr.text = "rr"
            wdrr.tail = "\n"
            frr = ET.SubElement(comprr, gn_brackets + "inputFile")
            frr.text = "Sobek_3b.fnm"
            frr.tail = "\n"
            myroot.insert(4, comprr)

            # FM to RR coupler
            couplerfmrr = ET.Element(gn_brackets + "coupler")
            couplerfmrr.attrib = {"name": "flow_to_rr"}
            couplerfmrr.tail = "\n"

            sourcefmrr = ET.SubElement(couplerfmrr, gn_brackets + "sourceComponent")
            sourcefmrr.text = "DFM"
            sourcefmrr.tail = "\n"

            targetfmrr = ET.SubElement(couplerfmrr, gn_brackets + "targetComponent")
            targetfmrr.text = "RR"
            targetfmrr.tail = "\n"

            for i in rr_model.external_forcings.boundary_nodes.keys():
                item = ET.SubElement(couplerfmrr, gn_brackets + "item")
                item.text = ""
                item.tail = "\n"

                source = ET.SubElement(item, gn_brackets + "sourceName")
                source.text = f"laterals/{i}/water_level"
                source.tail = "\n"

                target = ET.SubElement(item, gn_brackets + "targetName")
                target.text = f"catchments/{i}/water_level"
                target.tail = "\n"

            logfmrr = ET.SubElement(couplerfmrr, gn_brackets + "logger")
            logfmrr.text = ""
            logfmrr.tail = "\n"

            wdlogfmrr = ET.SubElement(logfmrr, gn_brackets + "workingDir")
            wdlogfmrr.text = "."
            wdlogfmrr.tail = "\n"

            flogfmrr = ET.SubElement(logfmrr, gn_brackets + "outputFile")
            flogfmrr.text = "dflowfm_to_rr.nc"
            flogfmrr.tail = "\n"
            myroot.append(couplerfmrr)

        if RTC:
            # RTC to FM coupler
            couplerrtcfm = ET.Element(gn_brackets + "coupler")
            couplerrtcfm.attrib = {"name": "rtc_to_flow"}
            couplerrtcfm.tail = "\n"

            sourcertcfm = ET.SubElement(couplerrtcfm, gn_brackets + "sourceComponent")
            sourcertcfm.text = "Real_Time_Control"
            sourcertcfm.tail = "\n"

            targetrtcfm = ET.SubElement(couplerrtcfm, gn_brackets + "targetComponent")
            targetrtcfm.text = "DFM"
            targetrtcfm.tail = "\n"

            rtcdict = {
                      "Crest level (s)": ["weirs","crestLevel"],
                      "Gate lower edge level (s)": ["orifices","gateLowerEdgeLevel"],
                      "Capacity (p)": ["pumps","capacity"],
                      "Valve opening (s)": ["culverts", "valveOpeningHeight"]
                      }

            rtc_to_flow_seen_pairs = set()
            rtc_to_flow_simple_count = 0
            rtc_to_flow_complex_count = 0

            for i in rtc_model.all_controllers.keys():
                svar=rtc_model.all_controllers[i]['steering_variable']
                if self._add_unique_coupler_item(
                    couplerrtcfm,
                    f"[Output]{i}/{svar}",
                    f"{rtcdict[svar][0]}/{i}/{rtcdict[svar][1]}",
                    rtc_to_flow_seen_pairs,
                    gn_brackets,
                ):
                    rtc_to_flow_simple_count += 1

            # check if there are user-specified controller that should be included
            if rtc_model.complex_controllers is not None:
                complex_config = rtc_model.complex_controllers["dimr_config"][0]
                for block in complex_config:
                    block_tag = block.tag.split("}", 1)[-1] if block.tag.startswith("{") else block.tag
                    if block_tag == "coupler":
                        if block.attrib["name"].lower().startswith("rtc_to_flow"):
                            for iblock in block:
                                iblock_tag = iblock.tag.split("}", 1)[-1] if iblock.tag.startswith("{") else iblock.tag
                                if iblock_tag == "item":
                                    source = iblock.find(".//{*}sourceName")
                                    target = iblock.find(".//{*}targetName")
                                    if self._add_unique_coupler_item(
                                        couplerrtcfm,
                                        source.text if source is not None else None,
                                        target.text if target is not None else None,
                                        rtc_to_flow_seen_pairs,
                                        gn_brackets,
                                    ):
                                        rtc_to_flow_complex_count += 1

            logger.info(
                "dimr_config.xml: Wrote rtc_to_flow coupler with %d items "
                "(%d simple, %d complex).",
                len(rtc_to_flow_seen_pairs),
                rtc_to_flow_simple_count,
                rtc_to_flow_complex_count,
            )

            logrtcfm = ET.SubElement(couplerrtcfm, gn_brackets + "logger")
            logrtcfm.text = ""
            logrtcfm.tail = "\n"

            wdlogrtcfm = ET.SubElement(logrtcfm, gn_brackets + "workingDir")
            wdlogrtcfm.text = "."
            wdlogrtcfm.tail = "\n"

            flogrtcfm = ET.SubElement(logrtcfm, gn_brackets + "outputFile")
            flogrtcfm.text = "rtc_to_dflowfm.nc"
            flogrtcfm.tail = "\n"

            myroot.append(couplerrtcfm)

            # the Fm to RTC coupler is not always needed. It is for PID controllers.
            coupler_exists = False
            flow_to_rtc_seen_pairs = set()
            flow_to_rtc_simple_count = 0
            flow_to_rtc_complex_count = 0
            if bool(rtc_model.pid_controllers) or bool(rtc_model.interval_controllers):
                couplerfmrtc = ET.Element(gn_brackets + "coupler")
                couplerfmrtc.attrib = {"name": "flow_to_rtc"}
                couplerfmrtc.tail = "\n"

                sourcefmrtc = ET.SubElement(
                    couplerfmrtc, gn_brackets + "sourceComponent"
                )
                sourcefmrtc.text = "DFM"
                sourcefmrtc.tail = "\n"

                targetfmrtc = ET.SubElement(
                    couplerfmrtc, gn_brackets + "targetComponent"
                )
                targetfmrtc.text = "Real_Time_Control"
                targetfmrtc.tail = "\n"

                for i in rtc_model.pid_controllers.keys():
                    if rtc_model.pid_controllers[i]['target_variable'] == 'Discharge (op)':
                        source_name = f"observations/{rtc_model.pid_controllers[i]['observation_point']}/discharge"
                    elif rtc_model.pid_controllers[i]['target_variable'] == 'Water level (op)':
                        source_name = f"observations/{rtc_model.pid_controllers[i]['observation_point']}/water_level"
                    else:
                        raise ValueError('Invalid target variable in controller: should bo discharge or water level.')
                    target_name = f"[Input]{rtc_model.pid_controllers[i]['observation_point']}/{rtc_model.pid_controllers[i]['target_variable']}"
                    if self._add_unique_coupler_item(
                        couplerfmrtc,
                        source_name,
                        target_name,
                        flow_to_rtc_seen_pairs,
                        gn_brackets,
                    ):
                        flow_to_rtc_simple_count += 1

                # Loop through all interval controllers
                for i in rtc_model.interval_controllers.keys():
                    if rtc_model.interval_controllers[i]["target_variable"] == "Discharge (op)":
                        source_name = f"observations/{rtc_model.interval_controllers[i]['observation_point']}/discharge"
                    elif rtc_model.interval_controllers[i]["target_variable"] == "Water level (op)":
                        source_name = f"observations/{rtc_model.interval_controllers[i]['observation_point']}/water_level"
                    else:
                        raise ValueError("Invalid target variable in controller: should be discharge or water level.")
                    target_name = f"[Input]{rtc_model.interval_controllers[i]['observation_point']}/{rtc_model.interval_controllers[i]['target_variable']}"
                    if self._add_unique_coupler_item(
                        couplerfmrtc,
                        source_name,
                        target_name,
                        flow_to_rtc_seen_pairs,
                        gn_brackets,
                    ):
                        flow_to_rtc_simple_count += 1

                coupler_exists = True
            # it could be that there are no PID/interval controllers, but there are complex controllers
            for block in complex_flow_to_rtc_blocks:
                if not coupler_exists:
                    couplerfmrtc = ET.Element(gn_brackets + "coupler")
                    couplerfmrtc.attrib = {"name": "flow_to_rtc"}
                    couplerfmrtc.tail = "\n"

                    sourcefmrtc = ET.SubElement(
                        couplerfmrtc, gn_brackets + "sourceComponent"
                    )
                    sourcefmrtc.text = "DFM"
                    sourcefmrtc.tail = "\n"

                    targetfmrtc = ET.SubElement(
                        couplerfmrtc, gn_brackets + "targetComponent"
                    )
                    targetfmrtc.text = "Real_Time_Control"
                    targetfmrtc.tail = "\n"

                    coupler_exists = True
                for iblock in block:
                    iblock_tag = iblock.tag.split("}", 1)[-1] if iblock.tag.startswith("{") else iblock.tag
                    if iblock_tag == "item":
                        source = iblock.find(".//{*}sourceName")
                        target = iblock.find(".//{*}targetName")
                        if self._add_unique_coupler_item(
                            couplerfmrtc,
                            source.text if source is not None else None,
                            target.text if target is not None else None,
                            flow_to_rtc_seen_pairs,
                            gn_brackets,
                        ):
                            flow_to_rtc_complex_count += 1
            # in any of those two cases, add the controllers and the logger
            if coupler_exists:
                logger.info(
                    "dimr_config.xml: Wrote flow_to_rtc coupler with %d items "
                    "(%d simple, %d complex).",
                    len(flow_to_rtc_seen_pairs),
                    flow_to_rtc_simple_count,
                    flow_to_rtc_complex_count,
                )

                logrtc = ET.SubElement(couplerfmrtc, gn_brackets + "logger")
                logrtc.text = ""
                logrtc.tail = "\n"

                wdlogrtc = ET.SubElement(logrtc, gn_brackets + "workingDir")
                wdlogrtc.text = "."
                wdlogrtc.tail = "\n"

                flogrtc = ET.SubElement(logrtc, gn_brackets + "outputFile")
                flogrtc.text = "dflowfm_to_rtc.nc"
                flogrtc.tail = "\n"

                myroot.append(couplerfmrtc)

        if RR:
            # finally the RR to FM coupler
            couplerrrfm = ET.Element(gn_brackets + "coupler")
            couplerrrfm.attrib = {"name": "rr_to_flow"}
            couplerrrfm.tail = "\n"

            sourcerrfm = ET.SubElement(couplerrrfm, gn_brackets + "sourceComponent")
            sourcerrfm.text = "RR"
            sourcerrfm.tail = "\n"

            targetrrfm = ET.SubElement(couplerrrfm, gn_brackets + "targetComponent")
            targetrrfm.text = "DFM"
            targetrrfm.tail = "\n"

            for i in rr_model.external_forcings.boundary_nodes.keys():
                item = ET.SubElement(couplerrrfm, gn_brackets + "item")
                item.text = ""
                item.tail = "\n"

                source = ET.SubElement(item, gn_brackets + "sourceName")
                source.text = f"catchments/{i}/water_discharge"
                source.tail = "\n"

                target = ET.SubElement(item, gn_brackets + "targetName")
                target.text = f"laterals/{i}/water_discharge"
                target.tail = "\n"

            logrrfm = ET.SubElement(couplerrrfm, gn_brackets + "logger")
            logrrfm.text = ""
            logrrfm.tail = "\n"

            wdlogrrfm = ET.SubElement(logrrfm, gn_brackets + "workingDir")
            wdlogrrfm.text = "."
            wdlogrrfm.tail = "\n"

            flogrrfm = ET.SubElement(logrrfm, gn_brackets + "outputFile")
            flogrrfm.text = "rr_to_dflowfm.nc"
            flogrrfm.tail = "\n"

            myroot.append(couplerrrfm)
        # use the function from DRTC to finish the file with the correct namespace
        DRTCModel.finish_file(myroot, configfile, self.output_path / "dimr_config.xml")

__init__(dimr_path: str = None, output_path: str | Path = None) -> None

Class initialization for the DIMR writer. This is a temporary fix to include writers for RR and RTC, that are not yet fully included in Hydrolib-core, or have bugs. Eventually, this file should be redundant as this functionality is fully included in Hydrolib-core

Parameters:

Name	Type	Description	Default
dimr_path	Union[str, Path]	description. Defaults to None.	None
output_path	Union[str, Path]	description. Defaults to None.	None

Source code in hydrolib/dhydamo/io/dimrwriter.py

@validate_arguments
def __init__(
    self, dimr_path: str = None, output_path: str | Path = None
) -> None:
    """Class initialization for the DIMR writer. This is a temporary fix to include writers for RR and RTC, that are not yet fully included in Hydrolib-core, or have bugs. Eventually, this file should be redundant as this functionality is fully included in Hydrolib-core

    Args:
        dimr_path (Union[str, Path], optional): _description_. Defaults to None.
        output_path (Union[str, Path], optional): _description_. Defaults to None.
    """
    if dimr_path is None:
        self.run_dimr = r"C:\Program Files\Deltares\D-HYDRO Suite 2022.04 1D2D\plugins\DeltaShell.Dimr\kernels\x64\dimr\scripts\run_dimr.bat"
    else:
        self.run_dimr = dimr_path

    if output_path is None:
        self.output_path = Path(os.path.abspath("."))
    else:
        self.output_path = output_path

    self.template_dir = Path(os.path.abspath("."))

add_crs(netcdf_path: str | Path = None) -> None

Reads the Netcdf file and addes the required attributes for a valid CRS.

Source code in hydrolib/dhydamo/io/dimrwriter.py

@validate_arguments
def add_crs(self, netcdf_path: str | Path = None) -> None:
    """Reads the Netcdf file and addes the required attributes for a valid CRS."""        
    if netcdf_path is None:
        netfile = list((self.output_path / 'dflowfm').glob('*.nc'))[0]
        if not netfile.exists():
            raise FileNotFoundError(f"Netcdf file not found in {self.output_path / 'dflowfm'}. Provide the correct path via the netcdf_path argument.")
    else:
        if netcdf_path.is_dir():
            netfile = list(netcdf_path.glob('*.nc'))[0]
        elif netcdf_path.is_file():
            netfile = netcdf_path  
        else:        
            raise FileNotFoundError(f"Netcdf file not found at {netcdf_path}.")   
    netf = nc.Dataset(netfile, 'r+')
    netf.Conventions =  'CF-1.8 UGRID-1.0 Deltares-0.10'
    proj = netf.createVariable('projected_coordinate_system','i4')
    proj.epsg = 28992
    proj.grid_mapping_name = "Amersfoort / RD New" #"Rijksdriehoeksstelsel"
    proj.longitude_of_prime_meridian = 0.0#; // double
    proj.semi_major_axis = 6377397.155#; // double
    proj.semi_minor_axis = 6356078.962818189#; // double
    proj.inverse_flattening = 299.1528128#; // double
    proj.proj4_params = "+prdoj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
    proj.EPSG_code = "EPSG:28992"
    proj.value = "value is equal to EPSG code"
    proj.wkt = "PROJCS[\"Amersfoort / RD New\",\n    GEOGCS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,\n"             
    netf.close()   

write_dimrconfig(fm: FMModel, rr_model: DRRModel = None, rtc_model: DRTCModel = None) -> None

Construct the contents of a DIMR-config file, using elements from FM, RTC and RR models.

Parameters:

Name	Type	Description	Default
fm	instance of class FM-model	Hydrolib-core model object. When only fm is set, only a 1d2d model is generated.	required
rr_model	instance of a DRRModel object	Class containig all information from the RR-model. Defaults to None.	None
rtc_model	instance of a DRTCmodel	Class containig all information from the RTC-model. Defaults to None.	None

Source code in hydrolib/dhydamo/io/dimrwriter.py

def write_dimrconfig(
    self, fm: FMModel, rr_model: DRRModel = None, rtc_model: DRTCModel = None
) -> None:
    """Construct the contents of a DIMR-config file, using elements from FM, RTC and RR models.

    Args:
        fm (instance of class FM-model): Hydrolib-core model object. When only fm is set, only a 1d2d model is generated.
        rr_model (instance of a DRRModel object, optional): Class containig all information from the RR-model. Defaults to None.
        rtc_model (instance of a DRTCmodel, optional): Class containig all information from the RTC-model. Defaults to None.
    """
    # initiate flags for RR and RTC blocks
    RR = False
    RTC = False
    if rr_model is not None:
        RR = True
    if rtc_model is not None:
        RTC = True

    generalname = "http://schemas.deltares.nl/dimr"  # Config
    xsi_name = "http://www.w3.org/2001/XMLSchema-instance"
    gn_brackets = "{" + generalname + "}"

    # registering namespaces
    ET.register_namespace("", generalname)
    ET.register_namespace("xsi", xsi_name)

    # Parsing xml file
    configfile = ET.parse(self.template_dir / "dimr_config.xml")

    myroot = configfile.getroot()
    myroot[1][1].text = fm.filepath.parts[-2]
    myroot[1][2].text = fm.filepath.name

    # control blocks
    control = ET.Element(gn_brackets + "control")
    control.text = ""
    control.tail = "\n"
    myroot.insert(1, control)

    if RR or RTC:
        parallel = ET.SubElement(control, gn_brackets + "parallel")
        parallel.tail = "\n"
        parallel.text = ""

    if RTC:
        # a flow_to_rtc coupler is only needed for PID/interval controllers,
        # or complex controllers whose dimr_config defines one themselves
        complex_flow_to_rtc_blocks = self._get_complex_coupler_blocks(
            rtc_model, "flow_to_rtc"
        )
        needs_flow_to_rtc_coupler = (
            bool(rtc_model.pid_controllers)
            or bool(rtc_model.interval_controllers)
            or bool(complex_flow_to_rtc_blocks)
        )

        startgrouprtc = ET.SubElement(parallel, gn_brackets + "startGroup")
        startgrouprtc.text = ""
        startgrouprtc.tail = "\n"

        timertc = ET.SubElement(startgrouprtc, gn_brackets + "time")
        timertc.text = f"{0} {rtc_model.time_settings['step']} {fm.time.tstop}"
        timertc.tail = "\n"
        if needs_flow_to_rtc_coupler:
            couplerrtcrr = ET.SubElement(startgrouprtc, gn_brackets + "coupler")
            couplerrtcrr.attrib = {"name": "flow_to_rtc"}
            couplerrtcrr.tail = "\n"

        startrtc = ET.SubElement(startgrouprtc, gn_brackets + "start")
        startrtc.attrib = {"name": "Real_Time_Control"}
        startrtc.tail = "\n"
        couplerrtc = ET.SubElement(startgrouprtc, gn_brackets + "coupler")
        couplerrtc.attrib = {"name": "rtc_to_flow"}
        couplerrtc.tail = "\n"

    if RR:
        startgrouprr = ET.SubElement(parallel, gn_brackets + "startGroup")
        startgrouprr.text = ""
        startgrouprr.tail = "\n"

        timerr = ET.SubElement(startgrouprr, gn_brackets + "time")
        timerr.text = (
            f"{0} {rr_model.d3b_parameters['Timestepsize']} {fm.time.tstop}"
        )
        timerr.tail = "\n"
        couplerrrf = ET.SubElement(startgrouprr, gn_brackets + "coupler")
        couplerrrf.attrib = {"name": "flow_to_rr"}
        couplerrrf.tail = "\n"
        startrr = ET.SubElement(startgrouprr, gn_brackets + "start")
        startrr.attrib = {"name": "RR"}
        startrr.tail = "\n"
        couplerfrr = ET.SubElement(startgrouprr, gn_brackets + "coupler")
        couplerfrr.attrib = {"name": "rr_to_flow"}
        couplerfrr.tail = "\n"
    if RR or RTC:
        startfm = ET.SubElement(parallel, gn_brackets + "start")
    else:
        startfm = ET.SubElement(control, gn_brackets + "start")
    startfm.attrib = {"name": "DFM"}
    startfm.tail = "\n"

    # component blocks
    if RTC:
        comprtc = ET.Element(gn_brackets + "component")
        comprtc.attrib = {"name": "Real_Time_Control"}
        comprtc.tail = "\n"
        librtc = ET.SubElement(comprtc, gn_brackets + "library")
        librtc.text = "FBCTools_BMI"
        librtc.tail = "\n"
        wdrtc = ET.SubElement(comprtc, gn_brackets + "workingDir")
        wdrtc.text = "rtc"
        wdrtc.tail = "\n"
        frtc = ET.SubElement(comprtc, gn_brackets + "inputFile")
        frtc.text = "."
        frtc.tail = "\n"
        myroot.insert(2, comprtc)

    if RR:
        comprr = ET.Element(gn_brackets + "component")
        comprr.attrib = {"name": "RR"}
        comprr.tail = "\n"
        librr = ET.SubElement(comprr, gn_brackets + "library")
        librr.text = "rr_dll"
        librr.tail = "\n"
        wdrr = ET.SubElement(comprr, gn_brackets + "workingDir")
        wdrr.text = "rr"
        wdrr.tail = "\n"
        frr = ET.SubElement(comprr, gn_brackets + "inputFile")
        frr.text = "Sobek_3b.fnm"
        frr.tail = "\n"
        myroot.insert(4, comprr)

        # FM to RR coupler
        couplerfmrr = ET.Element(gn_brackets + "coupler")
        couplerfmrr.attrib = {"name": "flow_to_rr"}
        couplerfmrr.tail = "\n"

        sourcefmrr = ET.SubElement(couplerfmrr, gn_brackets + "sourceComponent")
        sourcefmrr.text = "DFM"
        sourcefmrr.tail = "\n"

        targetfmrr = ET.SubElement(couplerfmrr, gn_brackets + "targetComponent")
        targetfmrr.text = "RR"
        targetfmrr.tail = "\n"

        for i in rr_model.external_forcings.boundary_nodes.keys():
            item = ET.SubElement(couplerfmrr, gn_brackets + "item")
            item.text = ""
            item.tail = "\n"

            source = ET.SubElement(item, gn_brackets + "sourceName")
            source.text = f"laterals/{i}/water_level"
            source.tail = "\n"

            target = ET.SubElement(item, gn_brackets + "targetName")
            target.text = f"catchments/{i}/water_level"
            target.tail = "\n"

        logfmrr = ET.SubElement(couplerfmrr, gn_brackets + "logger")
        logfmrr.text = ""
        logfmrr.tail = "\n"

        wdlogfmrr = ET.SubElement(logfmrr, gn_brackets + "workingDir")
        wdlogfmrr.text = "."
        wdlogfmrr.tail = "\n"

        flogfmrr = ET.SubElement(logfmrr, gn_brackets + "outputFile")
        flogfmrr.text = "dflowfm_to_rr.nc"
        flogfmrr.tail = "\n"
        myroot.append(couplerfmrr)

    if RTC:
        # RTC to FM coupler
        couplerrtcfm = ET.Element(gn_brackets + "coupler")
        couplerrtcfm.attrib = {"name": "rtc_to_flow"}
        couplerrtcfm.tail = "\n"

        sourcertcfm = ET.SubElement(couplerrtcfm, gn_brackets + "sourceComponent")
        sourcertcfm.text = "Real_Time_Control"
        sourcertcfm.tail = "\n"

        targetrtcfm = ET.SubElement(couplerrtcfm, gn_brackets + "targetComponent")
        targetrtcfm.text = "DFM"
        targetrtcfm.tail = "\n"

        rtcdict = {
                  "Crest level (s)": ["weirs","crestLevel"],
                  "Gate lower edge level (s)": ["orifices","gateLowerEdgeLevel"],
                  "Capacity (p)": ["pumps","capacity"],
                  "Valve opening (s)": ["culverts", "valveOpeningHeight"]
                  }

        rtc_to_flow_seen_pairs = set()
        rtc_to_flow_simple_count = 0
        rtc_to_flow_complex_count = 0

        for i in rtc_model.all_controllers.keys():
            svar=rtc_model.all_controllers[i]['steering_variable']
            if self._add_unique_coupler_item(
                couplerrtcfm,
                f"[Output]{i}/{svar}",
                f"{rtcdict[svar][0]}/{i}/{rtcdict[svar][1]}",
                rtc_to_flow_seen_pairs,
                gn_brackets,
            ):
                rtc_to_flow_simple_count += 1

        # check if there are user-specified controller that should be included
        if rtc_model.complex_controllers is not None:
            complex_config = rtc_model.complex_controllers["dimr_config"][0]
            for block in complex_config:
                block_tag = block.tag.split("}", 1)[-1] if block.tag.startswith("{") else block.tag
                if block_tag == "coupler":
                    if block.attrib["name"].lower().startswith("rtc_to_flow"):
                        for iblock in block:
                            iblock_tag = iblock.tag.split("}", 1)[-1] if iblock.tag.startswith("{") else iblock.tag
                            if iblock_tag == "item":
                                source = iblock.find(".//{*}sourceName")
                                target = iblock.find(".//{*}targetName")
                                if self._add_unique_coupler_item(
                                    couplerrtcfm,
                                    source.text if source is not None else None,
                                    target.text if target is not None else None,
                                    rtc_to_flow_seen_pairs,
                                    gn_brackets,
                                ):
                                    rtc_to_flow_complex_count += 1

        logger.info(
            "dimr_config.xml: Wrote rtc_to_flow coupler with %d items "
            "(%d simple, %d complex).",
            len(rtc_to_flow_seen_pairs),
            rtc_to_flow_simple_count,
            rtc_to_flow_complex_count,
        )

        logrtcfm = ET.SubElement(couplerrtcfm, gn_brackets + "logger")
        logrtcfm.text = ""
        logrtcfm.tail = "\n"

        wdlogrtcfm = ET.SubElement(logrtcfm, gn_brackets + "workingDir")
        wdlogrtcfm.text = "."
        wdlogrtcfm.tail = "\n"

        flogrtcfm = ET.SubElement(logrtcfm, gn_brackets + "outputFile")
        flogrtcfm.text = "rtc_to_dflowfm.nc"
        flogrtcfm.tail = "\n"

        myroot.append(couplerrtcfm)

        # the Fm to RTC coupler is not always needed. It is for PID controllers.
        coupler_exists = False
        flow_to_rtc_seen_pairs = set()
        flow_to_rtc_simple_count = 0
        flow_to_rtc_complex_count = 0
        if bool(rtc_model.pid_controllers) or bool(rtc_model.interval_controllers):
            couplerfmrtc = ET.Element(gn_brackets + "coupler")
            couplerfmrtc.attrib = {"name": "flow_to_rtc"}
            couplerfmrtc.tail = "\n"

            sourcefmrtc = ET.SubElement(
                couplerfmrtc, gn_brackets + "sourceComponent"
            )
            sourcefmrtc.text = "DFM"
            sourcefmrtc.tail = "\n"

            targetfmrtc = ET.SubElement(
                couplerfmrtc, gn_brackets + "targetComponent"
            )
            targetfmrtc.text = "Real_Time_Control"
            targetfmrtc.tail = "\n"

            for i in rtc_model.pid_controllers.keys():
                if rtc_model.pid_controllers[i]['target_variable'] == 'Discharge (op)':
                    source_name = f"observations/{rtc_model.pid_controllers[i]['observation_point']}/discharge"
                elif rtc_model.pid_controllers[i]['target_variable'] == 'Water level (op)':
                    source_name = f"observations/{rtc_model.pid_controllers[i]['observation_point']}/water_level"
                else:
                    raise ValueError('Invalid target variable in controller: should bo discharge or water level.')
                target_name = f"[Input]{rtc_model.pid_controllers[i]['observation_point']}/{rtc_model.pid_controllers[i]['target_variable']}"
                if self._add_unique_coupler_item(
                    couplerfmrtc,
                    source_name,
                    target_name,
                    flow_to_rtc_seen_pairs,
                    gn_brackets,
                ):
                    flow_to_rtc_simple_count += 1

            # Loop through all interval controllers
            for i in rtc_model.interval_controllers.keys():
                if rtc_model.interval_controllers[i]["target_variable"] == "Discharge (op)":
                    source_name = f"observations/{rtc_model.interval_controllers[i]['observation_point']}/discharge"
                elif rtc_model.interval_controllers[i]["target_variable"] == "Water level (op)":
                    source_name = f"observations/{rtc_model.interval_controllers[i]['observation_point']}/water_level"
                else:
                    raise ValueError("Invalid target variable in controller: should be discharge or water level.")
                target_name = f"[Input]{rtc_model.interval_controllers[i]['observation_point']}/{rtc_model.interval_controllers[i]['target_variable']}"
                if self._add_unique_coupler_item(
                    couplerfmrtc,
                    source_name,
                    target_name,
                    flow_to_rtc_seen_pairs,
                    gn_brackets,
                ):
                    flow_to_rtc_simple_count += 1

            coupler_exists = True
        # it could be that there are no PID/interval controllers, but there are complex controllers
        for block in complex_flow_to_rtc_blocks:
            if not coupler_exists:
                couplerfmrtc = ET.Element(gn_brackets + "coupler")
                couplerfmrtc.attrib = {"name": "flow_to_rtc"}
                couplerfmrtc.tail = "\n"

                sourcefmrtc = ET.SubElement(
                    couplerfmrtc, gn_brackets + "sourceComponent"
                )
                sourcefmrtc.text = "DFM"
                sourcefmrtc.tail = "\n"

                targetfmrtc = ET.SubElement(
                    couplerfmrtc, gn_brackets + "targetComponent"
                )
                targetfmrtc.text = "Real_Time_Control"
                targetfmrtc.tail = "\n"

                coupler_exists = True
            for iblock in block:
                iblock_tag = iblock.tag.split("}", 1)[-1] if iblock.tag.startswith("{") else iblock.tag
                if iblock_tag == "item":
                    source = iblock.find(".//{*}sourceName")
                    target = iblock.find(".//{*}targetName")
                    if self._add_unique_coupler_item(
                        couplerfmrtc,
                        source.text if source is not None else None,
                        target.text if target is not None else None,
                        flow_to_rtc_seen_pairs,
                        gn_brackets,
                    ):
                        flow_to_rtc_complex_count += 1
        # in any of those two cases, add the controllers and the logger
        if coupler_exists:
            logger.info(
                "dimr_config.xml: Wrote flow_to_rtc coupler with %d items "
                "(%d simple, %d complex).",
                len(flow_to_rtc_seen_pairs),
                flow_to_rtc_simple_count,
                flow_to_rtc_complex_count,
            )

            logrtc = ET.SubElement(couplerfmrtc, gn_brackets + "logger")
            logrtc.text = ""
            logrtc.tail = "\n"

            wdlogrtc = ET.SubElement(logrtc, gn_brackets + "workingDir")
            wdlogrtc.text = "."
            wdlogrtc.tail = "\n"

            flogrtc = ET.SubElement(logrtc, gn_brackets + "outputFile")
            flogrtc.text = "dflowfm_to_rtc.nc"
            flogrtc.tail = "\n"

            myroot.append(couplerfmrtc)

    if RR:
        # finally the RR to FM coupler
        couplerrrfm = ET.Element(gn_brackets + "coupler")
        couplerrrfm.attrib = {"name": "rr_to_flow"}
        couplerrrfm.tail = "\n"

        sourcerrfm = ET.SubElement(couplerrrfm, gn_brackets + "sourceComponent")
        sourcerrfm.text = "RR"
        sourcerrfm.tail = "\n"

        targetrrfm = ET.SubElement(couplerrrfm, gn_brackets + "targetComponent")
        targetrrfm.text = "DFM"
        targetrrfm.tail = "\n"

        for i in rr_model.external_forcings.boundary_nodes.keys():
            item = ET.SubElement(couplerrrfm, gn_brackets + "item")
            item.text = ""
            item.tail = "\n"

            source = ET.SubElement(item, gn_brackets + "sourceName")
            source.text = f"catchments/{i}/water_discharge"
            source.tail = "\n"

            target = ET.SubElement(item, gn_brackets + "targetName")
            target.text = f"laterals/{i}/water_discharge"
            target.tail = "\n"

        logrrfm = ET.SubElement(couplerrrfm, gn_brackets + "logger")
        logrrfm.text = ""
        logrrfm.tail = "\n"

        wdlogrrfm = ET.SubElement(logrrfm, gn_brackets + "workingDir")
        wdlogrrfm.text = "."
        wdlogrrfm.tail = "\n"

        flogrrfm = ET.SubElement(logrrfm, gn_brackets + "outputFile")
        flogrrfm.text = "rr_to_dflowfm.nc"
        flogrrfm.tail = "\n"

        myroot.append(couplerrrfm)
    # use the function from DRTC to finish the file with the correct namespace
    DRTCModel.finish_file(myroot, configfile, self.output_path / "dimr_config.xml")

write_runbat(debuglevel=6, runlog=None) -> None

Generates a run.bat to run the model in DIMR. The path to the executable is provided by the user in the class initialization, or set to the default D-Hydro installation folder.

Source code in hydrolib/dhydamo/io/dimrwriter.py

@validate_arguments
def write_runbat(self, debuglevel=6, runlog=None) -> None:
    """Generates a run.bat to run the model in DIMR. The path to the executable is provided by the user in the class initialization, or set to the default D-Hydro installation folder."""

    with open(os.path.join(self.output_path, "run.bat"), "w") as f:
        f.write("@ echo off\n")
        f.write("set OMP_NUM_THREADS=2\n")
        if runlog is None:
            f.write('call "' + str(self.run_dimr) + '" -d ' + str(debuglevel) + '\n')
        else:
            f.write('call "' + str(self.run_dimr) + '" -d ' + str(debuglevel) + ' > ' + str(runlog) + '\n')

drrreader

ExternalForcingsIO

Source code in hydrolib/dhydamo/io/drrreader.py

class ExternalForcingsIO:
    def __init__(self, external_forcings):
        self.external_forcings = external_forcings

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def seepage_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        seepage_folder: Path | str,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Perform zonal statistics to derive seepage time series per catchment. Time steps are derived from the data

        Args:
            catchments (ExtendedGeoDataFrame): catchment areas
            seepage_folder (str): folder where the seepage rasters are stored
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
        """
        warnings.filterwarnings("ignore")
        rasters = []
        times = []
        affine = raster_crs = nodata = None
        for file in sorted(os.listdir(seepage_folder)):
            path = os.path.join(seepage_folder, file)
            if file.endswith('.idf'):
                # Lazy: zonal_stats_per_raster reads it when its turn comes.
                dataset = idfreader.open(path).squeeze()
                if affine is None:
                    affine = idfreader.affine_from_idf(dataset)
                    raster_crs = dataset.attrs.get("crs")
                    nodata = dataset.attrs.get("nodata", np.nan)
                elif idfreader.affine_from_idf(dataset) != affine:
                    raise ValueError(f"IDF seepage file {file} differs in grid from the others.")
                rasters.append(dataset)
                times.append(pd.Timestamp(dataset["time"].values))
            else:
                rasters.append(path)
                times.append(self.external_forcings.drrmodel._time_from_filename(path))
        stats = zonal.zonal_stats_per_raster(
            gpd.GeoDataFrame(catchments),
            rasters,
            "mean",
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            affine=affine,
            raster_crs=raster_crs,
            nodata=nodata,
            progress=True,
        )
        result = pd.DataFrame(
            stats.to_numpy().T,
            index=times,
            columns=["sep_" + str(cat) for cat in catchments.code],
        )
        if affine is not None:
            # if an NHI model (IDF files) is used, convert units from m3 to mm/d
            result = (result / (1e-3 * (affine[0] * -affine[4]))) / (
                    (times[2] - times[1]).total_seconds() / 86400.0
            )
        [self.external_forcings.add_seepage(*sep) for sep in result.items()]


    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def precip_from_input(
        self,
        areas: ExtendedGeoDataFrame,
        precip_folder: Path | str | None = None,
        precip_file: Path | str | None = None,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Create time series of precipitation for every meteo_area, based on zonal statistics from rasters.

        Args:
            areas (ExtendedGeoDataFrame): meteo areas for which time series are created
            precip_folder (str, optional): folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.
            precip_file (str, optional): existing meteo-file, which is used if available.
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
        """
        if precip_file is None:
            warnings.filterwarnings("ignore")
            paths = [
                os.path.join(precip_folder, file)
                for file in sorted(os.listdir(precip_folder))
            ]
            stats = zonal.zonal_stats_per_raster(
                gpd.GeoDataFrame(areas),
                paths,
                "mean",
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                progress=True,
            )
            result = pd.DataFrame(
                stats.to_numpy().T,
                index=[self.external_forcings.drrmodel._time_from_filename(path) for path in paths],
                columns=["ms_" + str(area) for area in areas.code],
            )
            [self.external_forcings.add_precip(*prec) for prec in result.items()]
        else:
            self.external_forcings.precip = str(precip_file)

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def evap_from_input(
        self,
        areas: ExtendedGeoDataFrame,
        evap_folder: Path | str | None = None,
        evap_file: Path | str | None = None,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Create time series of evaporation for every meteo_area, based on zonal statistics from rasters.

        Args:
            areas (ExtendedGeoDataFrame): meteo areas for which time series are created
            evap_folder (str, optional): folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.
            evap_file (str, optional): existing meteo-file, which is used if available.
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
        """
        if evap_file is None:
            warnings.filterwarnings("ignore")
            # aggregated evap
            # areas['dissolve'] = 1
            # agg_areas = areas.iloc[0:len(areas),:].dissolve(by='dissolve',aggfunc='mean')
            paths = [
                os.path.join(evap_folder, file)
                for file in sorted(os.listdir(evap_folder))
            ]
            stats = zonal.zonal_stats_per_raster(
                gpd.GeoDataFrame(areas),
                paths,
                "mean",
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                progress=True,
            )
            result = pd.DataFrame(
                stats.to_numpy().T,
                index=[self.external_forcings.drrmodel._time_from_filename(path) for path in paths],
                columns=["ms_" + str(area) for area in areas.code],
            )
            [self.external_forcings.add_evap(*evap) for evap in result.items()]
        else:
            self.external_forcings.evap = str(evap_file)

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def boundary_from_input(
        self,
        boundary_nodes: ExtendedGeoDataFrame,
        catchments: ExtendedGeoDataFrame,
        drrmodel,
        overflows: ExtendedGeoDataFrame = None,
        greenhouse_laterals: ExtendedGeoDataFrame=None
    ) -> None:
        """Generate RR-boundary nodes to link to RR-nodes and to FM-laterals.

        Args:
            boundary_nodes (ExtendedGeoDataFrame): boundary nodes
            catchments (ExtendedGeoDataFrame): catchment areas associated with them
            drrmodel (_type_): drrmodel object
            overflows (ExtendedGeoDataFrame, optional): overflow locations, if applicable. Defaults to None.
            greenhouse_laterals (ExtendedGeoDataFrame, optional): overflow locations, if applicable. Defaults to None.

        """
        # find the catchments that have no area attached and no nodes that will be attached to the boundary
        not_occurring = []
        for cat in catchments.itertuples():
            occurs = False
            if cat.boundary_node in [
                val["boundary_node"]
                for val in drrmodel.unpaved.unp_nodes.values()
                if np.sum([float(d) for d in val["ar"].split(" ")]) > 0.0
            ]:
                occurs = True
            if cat.boundary_node in [
                val["boundary_node"]
                for val in drrmodel.paved.pav_nodes.values()
                if float(val["ar"]) > 0.0
            ]:
                occurs = True
            if cat.boundary_node in [
                val["boundary_node"]
                for val in drrmodel.greenhouse.gh_nodes.values()
                if float(val["ar"]) > 0.0
            ]:
                occurs = True
            if cat.boundary_node in [
                val["boundary_node"]
                for val in drrmodel.openwater.ow_nodes.values()
                if float(val["ar"]) > 0.0
            ]:
                occurs = True
            if not occurs:
                not_occurring.append(cat.boundary_node)


        drop_idx = catchments[catchments.boundary_node.isin(not_occurring)].index.to_list()
        if drop_idx:
            logger.warning(
                "%d catchments removed because of an area of 0 m2.",
                len(drop_idx),
            )
            # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
            catchments.drop(drop_idx, inplace=True)  # noqa: PD002

        catchment_lat = RelationsRR.CATCHMENT_LATERAL
        for cat in catchments.itertuples():
            lateral_id = getattr(cat, catchment_lat.foreign_key)
            lat = boundary_nodes[
                boundary_nodes[catchment_lat.parent_key] == lateral_id
            ]
            if lat.empty:
                logger.warning(
                    f"{cat.code} not connected to a boundary node. Skipping."
                )
                self.external_forcings.add_boundary_node(id="0.0", px="0.0", py="0.0")
                continue
            self.external_forcings.add_boundary_node(
                id=f"lat_{cat.code}",
                px=str(lat.geometry.x.iloc[0]).strip(),
                py=str(lat.geometry.y.iloc[0]).strip(),
            )
        if overflows is not None:
            logger.info("Adding overflows to the boundary nodes.")
            for ovf in overflows.itertuples():
                self.external_forcings.add_boundary_node(
                    id=str(ovf.code),
                    px=str(ovf.geometry.coords[0][0]),
                    py=str(ovf.geometry.coords[0][1]),
                )
        if greenhouse_laterals is not None:
            logger.info("Adding greenhouse_laterals to the boundary nodes.")
            for gh in greenhouse_laterals.itertuples():
                self.external_forcings.add_boundary_node(
                    id=str(gh.code),
                    px=str(gh.geometry.coords[0][0]),
                    py=str(gh.geometry.coords[0][1]),
                )

boundary_from_input(boundary_nodes: ExtendedGeoDataFrame, catchments: ExtendedGeoDataFrame, drrmodel, overflows: ExtendedGeoDataFrame = None, greenhouse_laterals: ExtendedGeoDataFrame = None) -> None

Generate RR-boundary nodes to link to RR-nodes and to FM-laterals.

Parameters:

Name	Type	Description	Default
boundary_nodes	ExtendedGeoDataFrame	boundary nodes	required
catchments	ExtendedGeoDataFrame	catchment areas associated with them	required
drrmodel	_type_	drrmodel object	required
overflows	ExtendedGeoDataFrame	overflow locations, if applicable. Defaults to None.	None
greenhouse_laterals	ExtendedGeoDataFrame	overflow locations, if applicable. Defaults to None.	None

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def boundary_from_input(
    self,
    boundary_nodes: ExtendedGeoDataFrame,
    catchments: ExtendedGeoDataFrame,
    drrmodel,
    overflows: ExtendedGeoDataFrame = None,
    greenhouse_laterals: ExtendedGeoDataFrame=None
) -> None:
    """Generate RR-boundary nodes to link to RR-nodes and to FM-laterals.

    Args:
        boundary_nodes (ExtendedGeoDataFrame): boundary nodes
        catchments (ExtendedGeoDataFrame): catchment areas associated with them
        drrmodel (_type_): drrmodel object
        overflows (ExtendedGeoDataFrame, optional): overflow locations, if applicable. Defaults to None.
        greenhouse_laterals (ExtendedGeoDataFrame, optional): overflow locations, if applicable. Defaults to None.

    """
    # find the catchments that have no area attached and no nodes that will be attached to the boundary
    not_occurring = []
    for cat in catchments.itertuples():
        occurs = False
        if cat.boundary_node in [
            val["boundary_node"]
            for val in drrmodel.unpaved.unp_nodes.values()
            if np.sum([float(d) for d in val["ar"].split(" ")]) > 0.0
        ]:
            occurs = True
        if cat.boundary_node in [
            val["boundary_node"]
            for val in drrmodel.paved.pav_nodes.values()
            if float(val["ar"]) > 0.0
        ]:
            occurs = True
        if cat.boundary_node in [
            val["boundary_node"]
            for val in drrmodel.greenhouse.gh_nodes.values()
            if float(val["ar"]) > 0.0
        ]:
            occurs = True
        if cat.boundary_node in [
            val["boundary_node"]
            for val in drrmodel.openwater.ow_nodes.values()
            if float(val["ar"]) > 0.0
        ]:
            occurs = True
        if not occurs:
            not_occurring.append(cat.boundary_node)


    drop_idx = catchments[catchments.boundary_node.isin(not_occurring)].index.to_list()
    if drop_idx:
        logger.warning(
            "%d catchments removed because of an area of 0 m2.",
            len(drop_idx),
        )
        # inplace: returning a copy of an ExtendedGeoDataFrame gives a plain GeoDataFrame
        catchments.drop(drop_idx, inplace=True)  # noqa: PD002

    catchment_lat = RelationsRR.CATCHMENT_LATERAL
    for cat in catchments.itertuples():
        lateral_id = getattr(cat, catchment_lat.foreign_key)
        lat = boundary_nodes[
            boundary_nodes[catchment_lat.parent_key] == lateral_id
        ]
        if lat.empty:
            logger.warning(
                f"{cat.code} not connected to a boundary node. Skipping."
            )
            self.external_forcings.add_boundary_node(id="0.0", px="0.0", py="0.0")
            continue
        self.external_forcings.add_boundary_node(
            id=f"lat_{cat.code}",
            px=str(lat.geometry.x.iloc[0]).strip(),
            py=str(lat.geometry.y.iloc[0]).strip(),
        )
    if overflows is not None:
        logger.info("Adding overflows to the boundary nodes.")
        for ovf in overflows.itertuples():
            self.external_forcings.add_boundary_node(
                id=str(ovf.code),
                px=str(ovf.geometry.coords[0][0]),
                py=str(ovf.geometry.coords[0][1]),
            )
    if greenhouse_laterals is not None:
        logger.info("Adding greenhouse_laterals to the boundary nodes.")
        for gh in greenhouse_laterals.itertuples():
            self.external_forcings.add_boundary_node(
                id=str(gh.code),
                px=str(gh.geometry.coords[0][0]),
                py=str(gh.geometry.coords[0][1]),
            )

evap_from_input(areas: ExtendedGeoDataFrame, evap_folder: Path | str | None = None, evap_file: Path | str | None = None, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Create time series of evaporation for every meteo_area, based on zonal statistics from rasters.

Parameters:

Name	Type	Description	Default
areas	ExtendedGeoDataFrame	meteo areas for which time series are created	required
evap_folder	str	folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.	None
evap_file	str	existing meteo-file, which is used if available.	None
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def evap_from_input(
    self,
    areas: ExtendedGeoDataFrame,
    evap_folder: Path | str | None = None,
    evap_file: Path | str | None = None,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Create time series of evaporation for every meteo_area, based on zonal statistics from rasters.

    Args:
        areas (ExtendedGeoDataFrame): meteo areas for which time series are created
        evap_folder (str, optional): folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.
        evap_file (str, optional): existing meteo-file, which is used if available.
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
    """
    if evap_file is None:
        warnings.filterwarnings("ignore")
        # aggregated evap
        # areas['dissolve'] = 1
        # agg_areas = areas.iloc[0:len(areas),:].dissolve(by='dissolve',aggfunc='mean')
        paths = [
            os.path.join(evap_folder, file)
            for file in sorted(os.listdir(evap_folder))
        ]
        stats = zonal.zonal_stats_per_raster(
            gpd.GeoDataFrame(areas),
            paths,
            "mean",
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            progress=True,
        )
        result = pd.DataFrame(
            stats.to_numpy().T,
            index=[self.external_forcings.drrmodel._time_from_filename(path) for path in paths],
            columns=["ms_" + str(area) for area in areas.code],
        )
        [self.external_forcings.add_evap(*evap) for evap in result.items()]
    else:
        self.external_forcings.evap = str(evap_file)

precip_from_input(areas: ExtendedGeoDataFrame, precip_folder: Path | str | None = None, precip_file: Path | str | None = None, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Create time series of precipitation for every meteo_area, based on zonal statistics from rasters.

Parameters:

Name	Type	Description	Default
areas	ExtendedGeoDataFrame	meteo areas for which time series are created	required
precip_folder	str	folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.	None
precip_file	str	existing meteo-file, which is used if available.	None
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def precip_from_input(
    self,
    areas: ExtendedGeoDataFrame,
    precip_folder: Path | str | None = None,
    precip_file: Path | str | None = None,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Create time series of precipitation for every meteo_area, based on zonal statistics from rasters.

    Args:
        areas (ExtendedGeoDataFrame): meteo areas for which time series are created
        precip_folder (str, optional): folder where precipitation rasters are stored. Only used if no precip_file is given. Defaults to None.
        precip_file (str, optional): existing meteo-file, which is used if available.
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
    """
    if precip_file is None:
        warnings.filterwarnings("ignore")
        paths = [
            os.path.join(precip_folder, file)
            for file in sorted(os.listdir(precip_folder))
        ]
        stats = zonal.zonal_stats_per_raster(
            gpd.GeoDataFrame(areas),
            paths,
            "mean",
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            progress=True,
        )
        result = pd.DataFrame(
            stats.to_numpy().T,
            index=[self.external_forcings.drrmodel._time_from_filename(path) for path in paths],
            columns=["ms_" + str(area) for area in areas.code],
        )
        [self.external_forcings.add_precip(*prec) for prec in result.items()]
    else:
        self.external_forcings.precip = str(precip_file)

seepage_from_input(catchments: ExtendedGeoDataFrame, seepage_folder: Path | str, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Perform zonal statistics to derive seepage time series per catchment. Time steps are derived from the data

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	catchment areas	required
seepage_folder	str	folder where the seepage rasters are stored	required
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def seepage_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    seepage_folder: Path | str,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Perform zonal statistics to derive seepage time series per catchment. Time steps are derived from the data

    Args:
        catchments (ExtendedGeoDataFrame): catchment areas
        seepage_folder (str): folder where the seepage rasters are stored
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
    """
    warnings.filterwarnings("ignore")
    rasters = []
    times = []
    affine = raster_crs = nodata = None
    for file in sorted(os.listdir(seepage_folder)):
        path = os.path.join(seepage_folder, file)
        if file.endswith('.idf'):
            # Lazy: zonal_stats_per_raster reads it when its turn comes.
            dataset = idfreader.open(path).squeeze()
            if affine is None:
                affine = idfreader.affine_from_idf(dataset)
                raster_crs = dataset.attrs.get("crs")
                nodata = dataset.attrs.get("nodata", np.nan)
            elif idfreader.affine_from_idf(dataset) != affine:
                raise ValueError(f"IDF seepage file {file} differs in grid from the others.")
            rasters.append(dataset)
            times.append(pd.Timestamp(dataset["time"].values))
        else:
            rasters.append(path)
            times.append(self.external_forcings.drrmodel._time_from_filename(path))
    stats = zonal.zonal_stats_per_raster(
        gpd.GeoDataFrame(catchments),
        rasters,
        "mean",
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
        affine=affine,
        raster_crs=raster_crs,
        nodata=nodata,
        progress=True,
    )
    result = pd.DataFrame(
        stats.to_numpy().T,
        index=times,
        columns=["sep_" + str(cat) for cat in catchments.code],
    )
    if affine is not None:
        # if an NHI model (IDF files) is used, convert units from m3 to mm/d
        result = (result / (1e-3 * (affine[0] * -affine[4]))) / (
                (times[2] - times[1]).total_seconds() / 86400.0
        )
    [self.external_forcings.add_seepage(*sep) for sep in result.items()]

GreenhouseIO

Source code in hydrolib/dhydamo/io/drrreader.py

class GreenhouseIO:
    def __init__(self, greenhouse):
        self.greenhouse = greenhouse

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def greenhouse_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        landuse: Path | str,
        surface_level: Path | str,
        roof_storage: StrictStr | float,
        meteo_areas: ExtendedGeoDataFrame,
        greenhouse_areas: ExtendedGeoDataFrame=None,
        greenhouse_laterals: ExtendedGeoDataFrame=None,
        basin_storage_class: int = 2,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Generate contents of an RR greenhouse node

        Args:
            catchments (ExtendedGeoDataFrame): catchment areas
            greenhouuse_areas (ExtendedGeoDataFrame): known set of greenhouse areas with attiributes
            greenhouse_laterals (ExtendedGeoDataFrame) : known set of outlet points for greenhouses
            landuse (str): filename of land use raster
            surface_level (str): filename of surface level raster
            roofstorage (Union): float for spatially uniform sewer storage (mm), or raster for distributed values
            meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
        """
        px_area, reference_crs = _raster_metadata(landuse)

        lu_counts = self.greenhouse.drrmodel.zonal_cache.zonal_category_counts(
            gpd.GeoDataFrame(catchments),
            landuse,
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )
        stats_elev = self.greenhouse.drrmodel.zonal_cache.zonal_stats(
            gpd.GeoDataFrame(catchments),
            surface_level,
            statistics=("median",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )
        if greenhouse_areas is not None:
            stats_elev_gh = zonal.zonal_stats(
                 gpd.GeoDataFrame(greenhouse_areas),
                 surface_level,
                 statistics=("median",),
                 min_coverage_frac=min_coverage_frac,
                 coverage_weight=coverage_weight,
                 strategy="feature-sequential",
            )

        # optional rasters
        if isinstance(roof_storage, (Path, str)):
            roofstors = zonal.zonal_stats(
                 gpd.GeoDataFrame(catchments),
                 roof_storage,
                 statistics=("mean",),
                 min_coverage_frac=min_coverage_frac,
                 coverage_weight=coverage_weight,
                 strategy="feature-sequential",
                 raster_crs=reference_crs,
            )
            if greenhouse_areas is not None:
                roofstors_gh = zonal.zonal_stats(
                     gpd.GeoDataFrame(greenhouse_areas),
                     roof_storage,
                     statistics=("mean",),
                     min_coverage_frac=min_coverage_frac,
                     coverage_weight=coverage_weight,
                     strategy="feature-sequential",
                     raster_crs=reference_crs,
                )

        gh_columns = [
            "id",
            "area",
            "surface_level",
            "roof_storage",
            "basin_storage_class",
            "meteo_area",
            "px",
            "py",
            "boundary_node",
        ]
        if greenhouse_areas is not None:
            for num, gh in enumerate(greenhouse_areas.itertuples()):
                # find corresponding meteo-station
                if pd.isna(stats_elev_gh.iloc[num]["median"]):
                    logger.warning(NO_RASTERDATA_WARNING, gh.code)
                    self.greenhouse.add_greenhouse(**dict.fromkeys(gh_columns, "0.0"))
                    continue
                tm = [
                    m
                    for m in meteo_areas.itertuples()
                    if m.geometry.contains(gh.geometry.centroid)
                ]
                ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

                elev = stats_elev_gh.iloc[num]["median"]
                if hasattr(gh, 'roof_storage_mm') and not np.isnan(gh.roof_storage_mm):
                    roof_storage_val = f"{gh.roof_storage_mm:.2f}"
                elif isinstance(roof_storage, float):
                    roof_storage_val = f"{roof_storage:.2f}"
                else:
                    roof_storage_val = f'{roofstors_gh.iloc[num]["mean"]:.2f}'
                if hasattr(gh, 'basin_storage_class') and not np.isnan(gh.basin_storage_class):
                    basin_storage_class_val = f"{gh.basin_storage_class:g}"
                else:
                    basin_storage_class_val = f'{basin_storage_class:g}'
                gh_area = RelationsRR.GREENHOUSE_LATERAL_AREA
                area_id = getattr(gh, gh_area.parent_key)
                related_laterals = greenhouse_laterals[
                    greenhouse_laterals[gh_area.foreign_key] == area_id
                ]
                latcode = related_laterals.code.to_numpy()[0]
                self.greenhouse.add_greenhouse(
                    id=str(gh.code),
                    area=gh.geometry.area,
                    surface_level=f"{elev:.2f}",
                    roof_storage=roof_storage_val,
                    basin_storage_class=basin_storage_class_val,
                    meteo_area=str(ms),
                    px=f"{gh.geometry.centroid.coords[0][0]:.0f}",
                    py=f"{gh.geometry.centroid.coords[0][1]:.0f}",
                    boundary_node=str(latcode),
                )

        for num, cat in enumerate(catchments.itertuples()):
            # if no rasterdata could be obtained for this catchment, skip it.
            if pd.isna(stats_elev.iloc[num]["median"]):
                logger.warning(NO_RASTERDATA_WARNING, cat.code)
                self.greenhouse.add_greenhouse(**dict.fromkeys(gh_columns, "0.0"))
                continue

            # find corresponding meteo-station
            tm = [
                m
                for m in meteo_areas.itertuples()
                if m.geometry.contains(cat.geometry.centroid)
            ]
            ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

            if greenhouse_areas is not None and cat.geometry.intersects(greenhouse_areas.geometry).any():
                intersection_area = cat.geometry.intersection(greenhouse_areas.geometry).area
                intersection_area = intersection_area[intersection_area > 0.].to_numpy()[0]                    
                landuse_counts = lu_counts.iloc[num]
                if landuse_counts.get(15, 0) > 0:
                    # divide area to subtract between greenhouses and the most occurring area                                                
                    logger.info(
                        "Catchment: %s: subtracting %s m2 from greenhouse area in landuse map.",
                        cat.code,
                        np.min([(landuse_counts[15] * px_area, intersection_area)]),
                    )
                    lu_counts.iloc[num, lu_counts.columns.get_loc(15)] = np.max(
                        [0.0, landuse_counts[15] - np.round(intersection_area / px_area)]
                    )

            elev = stats_elev.iloc[num]["median"]
            roof_storage_val = (
                f"{roof_storage:.2f}"
                if isinstance(roof_storage, float)
                else f'{roofstors.iloc[num]["mean"]:.2f}'
            )
            self.greenhouse.add_greenhouse(
                id=str(cat.code),
                area=(
                    str(lu_counts.iloc[num].get(15, 0) * px_area)
                ),
                surface_level=f"{elev:.2f}",
                roof_storage=roof_storage_val,
                basin_storage_class=f"{basin_storage_class:g}",
                meteo_area=str(ms),
                px=f"{cat.geometry.centroid.coords[0][0]+20:.0f}",
                py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
                boundary_node=cat.boundary_node,
            )

greenhouse_from_input(catchments: ExtendedGeoDataFrame, landuse: Path | str, surface_level: Path | str, roof_storage: StrictStr | float, meteo_areas: ExtendedGeoDataFrame, greenhouse_areas: ExtendedGeoDataFrame = None, greenhouse_laterals: ExtendedGeoDataFrame = None, basin_storage_class: int = 2, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Generate contents of an RR greenhouse node

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	catchment areas	required
greenhouuse_areas	ExtendedGeoDataFrame	known set of greenhouse areas with attiributes	required
greenhouse_laterals (ExtendedGeoDataFrame)		known set of outlet points for greenhouses	required
landuse	str	filename of land use raster	required
surface_level	str	filename of surface level raster	required
roofstorage	Union	float for spatially uniform sewer storage (mm), or raster for distributed values	required
meteo_areas	ExtendedGeoDataFrame	meteo areas, for each station a meteo time series is assigned	required
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def greenhouse_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    landuse: Path | str,
    surface_level: Path | str,
    roof_storage: StrictStr | float,
    meteo_areas: ExtendedGeoDataFrame,
    greenhouse_areas: ExtendedGeoDataFrame=None,
    greenhouse_laterals: ExtendedGeoDataFrame=None,
    basin_storage_class: int = 2,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Generate contents of an RR greenhouse node

    Args:
        catchments (ExtendedGeoDataFrame): catchment areas
        greenhouuse_areas (ExtendedGeoDataFrame): known set of greenhouse areas with attiributes
        greenhouse_laterals (ExtendedGeoDataFrame) : known set of outlet points for greenhouses
        landuse (str): filename of land use raster
        surface_level (str): filename of surface level raster
        roofstorage (Union): float for spatially uniform sewer storage (mm), or raster for distributed values
        meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
    """
    px_area, reference_crs = _raster_metadata(landuse)

    lu_counts = self.greenhouse.drrmodel.zonal_cache.zonal_category_counts(
        gpd.GeoDataFrame(catchments),
        landuse,
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )
    stats_elev = self.greenhouse.drrmodel.zonal_cache.zonal_stats(
        gpd.GeoDataFrame(catchments),
        surface_level,
        statistics=("median",),
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )
    if greenhouse_areas is not None:
        stats_elev_gh = zonal.zonal_stats(
             gpd.GeoDataFrame(greenhouse_areas),
             surface_level,
             statistics=("median",),
             min_coverage_frac=min_coverage_frac,
             coverage_weight=coverage_weight,
             strategy="feature-sequential",
        )

    # optional rasters
    if isinstance(roof_storage, (Path, str)):
        roofstors = zonal.zonal_stats(
             gpd.GeoDataFrame(catchments),
             roof_storage,
             statistics=("mean",),
             min_coverage_frac=min_coverage_frac,
             coverage_weight=coverage_weight,
             strategy="feature-sequential",
             raster_crs=reference_crs,
        )
        if greenhouse_areas is not None:
            roofstors_gh = zonal.zonal_stats(
                 gpd.GeoDataFrame(greenhouse_areas),
                 roof_storage,
                 statistics=("mean",),
                 min_coverage_frac=min_coverage_frac,
                 coverage_weight=coverage_weight,
                 strategy="feature-sequential",
                 raster_crs=reference_crs,
            )

    gh_columns = [
        "id",
        "area",
        "surface_level",
        "roof_storage",
        "basin_storage_class",
        "meteo_area",
        "px",
        "py",
        "boundary_node",
    ]
    if greenhouse_areas is not None:
        for num, gh in enumerate(greenhouse_areas.itertuples()):
            # find corresponding meteo-station
            if pd.isna(stats_elev_gh.iloc[num]["median"]):
                logger.warning(NO_RASTERDATA_WARNING, gh.code)
                self.greenhouse.add_greenhouse(**dict.fromkeys(gh_columns, "0.0"))
                continue
            tm = [
                m
                for m in meteo_areas.itertuples()
                if m.geometry.contains(gh.geometry.centroid)
            ]
            ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

            elev = stats_elev_gh.iloc[num]["median"]
            if hasattr(gh, 'roof_storage_mm') and not np.isnan(gh.roof_storage_mm):
                roof_storage_val = f"{gh.roof_storage_mm:.2f}"
            elif isinstance(roof_storage, float):
                roof_storage_val = f"{roof_storage:.2f}"
            else:
                roof_storage_val = f'{roofstors_gh.iloc[num]["mean"]:.2f}'
            if hasattr(gh, 'basin_storage_class') and not np.isnan(gh.basin_storage_class):
                basin_storage_class_val = f"{gh.basin_storage_class:g}"
            else:
                basin_storage_class_val = f'{basin_storage_class:g}'
            gh_area = RelationsRR.GREENHOUSE_LATERAL_AREA
            area_id = getattr(gh, gh_area.parent_key)
            related_laterals = greenhouse_laterals[
                greenhouse_laterals[gh_area.foreign_key] == area_id
            ]
            latcode = related_laterals.code.to_numpy()[0]
            self.greenhouse.add_greenhouse(
                id=str(gh.code),
                area=gh.geometry.area,
                surface_level=f"{elev:.2f}",
                roof_storage=roof_storage_val,
                basin_storage_class=basin_storage_class_val,
                meteo_area=str(ms),
                px=f"{gh.geometry.centroid.coords[0][0]:.0f}",
                py=f"{gh.geometry.centroid.coords[0][1]:.0f}",
                boundary_node=str(latcode),
            )

    for num, cat in enumerate(catchments.itertuples()):
        # if no rasterdata could be obtained for this catchment, skip it.
        if pd.isna(stats_elev.iloc[num]["median"]):
            logger.warning(NO_RASTERDATA_WARNING, cat.code)
            self.greenhouse.add_greenhouse(**dict.fromkeys(gh_columns, "0.0"))
            continue

        # find corresponding meteo-station
        tm = [
            m
            for m in meteo_areas.itertuples()
            if m.geometry.contains(cat.geometry.centroid)
        ]
        ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

        if greenhouse_areas is not None and cat.geometry.intersects(greenhouse_areas.geometry).any():
            intersection_area = cat.geometry.intersection(greenhouse_areas.geometry).area
            intersection_area = intersection_area[intersection_area > 0.].to_numpy()[0]                    
            landuse_counts = lu_counts.iloc[num]
            if landuse_counts.get(15, 0) > 0:
                # divide area to subtract between greenhouses and the most occurring area                                                
                logger.info(
                    "Catchment: %s: subtracting %s m2 from greenhouse area in landuse map.",
                    cat.code,
                    np.min([(landuse_counts[15] * px_area, intersection_area)]),
                )
                lu_counts.iloc[num, lu_counts.columns.get_loc(15)] = np.max(
                    [0.0, landuse_counts[15] - np.round(intersection_area / px_area)]
                )

        elev = stats_elev.iloc[num]["median"]
        roof_storage_val = (
            f"{roof_storage:.2f}"
            if isinstance(roof_storage, float)
            else f'{roofstors.iloc[num]["mean"]:.2f}'
        )
        self.greenhouse.add_greenhouse(
            id=str(cat.code),
            area=(
                str(lu_counts.iloc[num].get(15, 0) * px_area)
            ),
            surface_level=f"{elev:.2f}",
            roof_storage=roof_storage_val,
            basin_storage_class=f"{basin_storage_class:g}",
            meteo_area=str(ms),
            px=f"{cat.geometry.centroid.coords[0][0]+20:.0f}",
            py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
            boundary_node=cat.boundary_node,
        )

OpenwaterIO

Source code in hydrolib/dhydamo/io/drrreader.py

class OpenwaterIO:
    def __init__(self, openwater):
        self.openwater = openwater

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def openwater_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        landuse: Path | str,
        meteo_areas: ExtendedGeoDataFrame,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Generate contents of an RR open water node.

        Args:
            catchments (ExtendedGeoDataFrame): catchment areas
            landuse (str): filename of landuse raster
            meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".

        Returns:
            _type_: _description_
        """
        px_area, _ = _raster_metadata(landuse)

        lu_counts = self.openwater.drrmodel.zonal_cache.zonal_category_counts(
            gpd.GeoDataFrame(catchments),
            landuse,
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        for num, cat in enumerate(catchments.itertuples()):
            # find corresponding meteo-station
            tm = [
                m
                for m in meteo_areas.itertuples()
                if m.geometry.contains(cat.geometry.centroid)
            ]
            ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

            self.openwater.add_openwater(
                id=str(cat.code),
                area=str(lu_counts.iloc[num].get(13, 0) * px_area),
                meteo_area=str(ms),
                px=f"{cat.geometry.centroid.coords[0][0]-20:.0f}",
                py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
                boundary_node=cat.boundary_node,
            )

openwater_from_input(catchments: ExtendedGeoDataFrame, landuse: Path | str, meteo_areas: ExtendedGeoDataFrame, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Generate contents of an RR open water node.

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	catchment areas	required
landuse	str	filename of landuse raster	required
meteo_areas	ExtendedGeoDataFrame	meteo areas, for each station a meteo time series is assigned	required
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Returns:

Name	Type	Description
_type_	None	description

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def openwater_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    landuse: Path | str,
    meteo_areas: ExtendedGeoDataFrame,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Generate contents of an RR open water node.

    Args:
        catchments (ExtendedGeoDataFrame): catchment areas
        landuse (str): filename of landuse raster
        meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".

    Returns:
        _type_: _description_
    """
    px_area, _ = _raster_metadata(landuse)

    lu_counts = self.openwater.drrmodel.zonal_cache.zonal_category_counts(
        gpd.GeoDataFrame(catchments),
        landuse,
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )

    for num, cat in enumerate(catchments.itertuples()):
        # find corresponding meteo-station
        tm = [
            m
            for m in meteo_areas.itertuples()
            if m.geometry.contains(cat.geometry.centroid)
        ]
        ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

        self.openwater.add_openwater(
            id=str(cat.code),
            area=str(lu_counts.iloc[num].get(13, 0) * px_area),
            meteo_area=str(ms),
            px=f"{cat.geometry.centroid.coords[0][0]-20:.0f}",
            py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
            boundary_node=cat.boundary_node,
        )

PavedIO

Source code in hydrolib/dhydamo/io/drrreader.py

class PavedIO:
    def __init__(self, paved):
        self.paved = paved

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def paved_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        landuse: StrictStr | Path,
        surface_level: StrictStr | Path,
        street_storage:StrictStr | Path | float,
        sewer_storage:StrictStr | Path | float,
        pump_capacity:StrictStr | Path | float,
        meteo_areas: ExtendedGeoDataFrame,
        overflows: ExtendedGeoDataFrame = None,
        sewer_areas: ExtendedGeoDataFrame = None,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ) -> None:
        """Generate contents of RR paved nodes

        Args:
            catchments (ExtendedGeoDataFrame): catchment areas
            landuse (str): filename of landuse raster
            surface_level (str): file name of suface level raster
            street_storage (Union): numeric for spatially uniform street storage (mm), or raster for distributed values
            sewer_storage (Union): numeric for spatially uniform sewer storage (mm), or raster for distributed values
            pump_capacity (Union): numeric for spatially uniform pump capaities (mm), or raster for distributed values
            meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
            overflows (ExtendedGeoDataFrame, optional): overflow locations. Defaults to None.
            sewer_areas (ExtendedGeoDataFrame, optional): sewer area locations. Defaults to None.
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".

        Returns:
            _type_: _description_
        """
        px_area, reference_crs = _raster_metadata(landuse)

        lu_counts = self.paved.drrmodel.zonal_cache.zonal_category_counts(
            gpd.GeoDataFrame(catchments),
            landuse,
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )
        stats_elev = self.paved.drrmodel.zonal_cache.zonal_stats(
            gpd.GeoDataFrame(catchments),
            surface_level,
            statistics=("median",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        if isinstance(street_storage, (Path, str)):
            str_stors = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                street_storage,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )

        if isinstance(sewer_storage, (Path, str)):
            sew_stors = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                sewer_storage,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )

        if isinstance(pump_capacity,  (Path, str)):
            # raster of POC in mm/h
            pump_caps = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                pump_capacity,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        paved_columns = [
            "id",
            "area",
            "surface_level",
            "street_storage",
            "sewer_storage",
            "pump_capacity",
            "meteo_area",
            "px",
            "py",
            "boundary_node",
        ]
        if sewer_areas is not None:
            # if the parameters ara rasters, do the zonal statistics per sewage area as well.
            if isinstance(street_storage,(Path, str)):
                str_stors_sa = zonal.zonal_stats(
                    gpd.GeoDataFrame(sewer_areas),
                    street_storage,
                    statistics=("mean",),
                    min_coverage_frac=min_coverage_frac,
                    coverage_weight=coverage_weight,
                    strategy="feature-sequential",
                    raster_crs=reference_crs,
                )
            if isinstance(sewer_storage, (Path, str)):
                sew_stors_sa = zonal.zonal_stats(
                    gpd.GeoDataFrame(sewer_areas),
                    sewer_storage,
                    statistics=("mean",),
                    min_coverage_frac=min_coverage_frac,
                    coverage_weight=coverage_weight,
                    strategy="feature-sequential",
                    raster_crs=reference_crs,
                )
            if isinstance(pump_capacity, (Path, str)):
                pump_caps_sa = zonal.zonal_stats(
                    gpd.GeoDataFrame(sewer_areas),
                    pump_capacity,
                    statistics=("mean",),
                    min_coverage_frac=min_coverage_frac,
                    coverage_weight=coverage_weight,
                    strategy="feature-sequential",
                    raster_crs=reference_crs,
                )
            stats_sa_elev = zonal.zonal_stats(
                gpd.GeoDataFrame(sewer_areas),
                surface_level,
                statistics=("median",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
            )
            sewer_lu_counts = zonal.zonal_category_counts(
                gpd.GeoDataFrame(sewer_areas),
                landuse,
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
            )

            # find the paved area in the sewer areas
            for isew, sew in enumerate(sewer_areas.itertuples()):
                pav_area = 0
                pixels = sewer_lu_counts.iloc[isew]
                pav_pixels = pixels.get(14, 0)
                if pav_pixels == 0:
                    logger.warning("No paved area in sewer area %s.", sew.code)
                    self.paved.add_paved(**dict.fromkeys(paved_columns, "0.0"))
                    continue
                pav_area += pav_pixels * px_area

                # subtract it fromthe total paved area in this catchment, make sure at least 0 remains
                # lu_counts[cat_ind][14.0] -=  pav_pixels
                # if lu_counts[cat_ind][14.0] < 0: lu_counts[cat_ind][14.0]  = 0

                elev = stats_sa_elev.iloc[isew]["median"]
                # find overflows related to this sewer area
                overflow_sa = RelationsRR.OVERFLOW_SEWER_AREA
                sewer_id = getattr(sew, overflow_sa.parent_key)
                ovf = overflows[overflows[overflow_sa.foreign_key] == sewer_id]
                for ov in ovf.itertuples():
                    # find corresponding meteo-station
                    tm = [
                        m
                        for m in meteo_areas.itertuples()
                        if m.geometry.contains(sew.geometry.centroid)
                    ]
                    ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

                    # if a float is given, a standard value is passed. If a string is given, a rastername is assumed to zonal statistics are applied.
                    if isinstance(street_storage, float):
                        street_storage_val = f"{street_storage:.2f}"
                    elif isinstance(street_storage, (Path, str)):
                        street_storage_val = f'{str_stors_sa.iloc[isew]["mean"]:.2f}'
                    else:
                        raise TypeError('Street_storage has the wrong datatype. It should be a filename (Path or string) or number (float or int).')

                    # three options: it can be an attribute of a sewer area, a uniform value or a raster
                    riool_berging_mm = getattr(sew, "riool_berging_mm", None)
                    if riool_berging_mm is None or np.isnan(riool_berging_mm) or not isinstance(riool_berging_mm, float):
                        if isinstance(sewer_storage, float):
                            sewer_storage_val = f"{sewer_storage:.2f}"
                        elif isinstance(sewer_storage, (Path, str)):
                            sewer_storage_val = f'{sew_stors_sa.iloc[isew]["mean"]:.2f}'
                        else:
                            raise TypeError('Sewer_storage has the wrong datatype. It should be a filename (Path or string) or number (float or int).')
                    else:
                        sewer_storage_val = f'{riool_berging_mm:.2f}'

                    # three options: it can be an attribute of a sewer area, a uniform value or a raster
                    riool_poc_m3s = getattr(sew, "riool_poc_m3s", None)
                    if riool_poc_m3s is None or np.isnan(riool_poc_m3s) or not isinstance(riool_poc_m3s, float):
                        if isinstance(pump_capacity, float):
                            # convert the value from mm/h to m3/s
                            pump_capacity_val = f"{pump_capacity * (float(pav_area) * ov.fractie) / (1000. * 3600.):.8f}"
                        elif isinstance(pump_capacity, (Path, str)):
                            # convert the value (extracted from the raster) from mm/h to m3/s
                            pump_capacity_val = f'{pump_caps_sa.iloc[isew]["mean"] * (float(pav_area) * ov.fractie) / (1000. * 3600.):.8f}'
                        else:
                            raise TypeError('Pump_capacity has the wrong datatype. It should be a filename (Path or string) or number (float or int).')
                    else:
                        # use the attribute value
                        pump_capacity_val = f'{riool_poc_m3s * ov.fractie:.8f}'

                    # add prefix to the overflow id to create the paved-node id
                    self.paved.add_paved(
                        id=str(ov.code),
                        area=str(pav_area * ov.fractie),
                        surface_level=f"{elev:.2f}",
                        street_storage=street_storage_val,
                        sewer_storage=sewer_storage_val,
                        pump_capacity=pump_capacity_val,
                        meteo_area=str(ms),
                        px=f"{ov.geometry.coords[0][0]+10:.0f}",
                        py=f"{ov.geometry.coords[0][1]:.0f}",
                        boundary_node=ov.code,
                    )

        if sewer_areas is not None:
            # Land use of the catchment parts outside the sewer areas, by
            # catchment position. Catchments that overlap no sewer area use
            # lu_counts.
            sewer_union = sewer_areas.union_all()
            geometries = catchments.geometry.reset_index(drop=True)
            overlaps = geometries.intersects(sewer_union)
            outside = geometries[overlaps].difference(sewer_union)
            outside = outside[outside.area > 0]
            outside_counts = pd.DataFrame()
            if not outside.empty:
                outside_counts = zonal.zonal_category_counts(
                    outside,
                    landuse,
                    min_coverage_frac=min_coverage_frac,
                    coverage_weight=coverage_weight,
                    strategy="feature-sequential",
                )

        for num, cat in enumerate(catchments.itertuples()):
            # if no rasterdata could be obtained for this catchment, skip it.
            if pd.isna(stats_elev.iloc[num]["median"]):
                logger.warning(NO_RASTERDATA_WARNING, cat.code)
                self.paved.add_paved(**dict.fromkeys(paved_columns, "0.0"))
                continue
            if sewer_areas is None or not overlaps[num]:
                pav_area = str(lu_counts.iloc[num].get(14, 0) * px_area)
            elif num in outside_counts.index:
                pav_area = str(outside_counts.loc[num].get(14, 0) * px_area)
            else:
                logger.info(
                    f"No paved area outside sewer area in catchments {cat.code}."
                )
                pav_area = 0.0

            # find corresponding meteo-station
            tm = [
                m
                for m in meteo_areas.itertuples()
                if m.geometry.contains(cat.geometry.centroid)
            ]
            ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

            elev = stats_elev.iloc[num]["median"]
            # if a float is given, a standard value is passed. If a string is given, a rastername is assumed to zonal statistics are applied.
            street_storage_val = (
                f"{street_storage:.2f}"
                if isinstance(street_storage, float)
                else f'{str_stors.iloc[num]["mean"]:.2f}'
            )
            sewer_storage_val = (
                f"{sewer_storage:.2f}"
                if isinstance(sewer_storage, float)
                else f'{sew_stors.iloc[num]["mean"]:.2f}'
            )
            pump_capacity_val = (
                f'{(pump_capacity * float(pav_area)) / (1000. * 3600.):.8f}'
                if isinstance(pump_capacity, float)
                else f'{pump_caps.iloc[num]["mean"] * (float(pav_area)) / (1000. * 3600.):.8f}'
            )
            self.paved.add_paved(
                id=str(cat.code),
                area=str(pav_area),
                surface_level=f"{elev:.2f}",
                street_storage=street_storage_val,
                sewer_storage=sewer_storage_val,
                pump_capacity=pump_capacity_val,
                meteo_area=str(ms),
                px=f"{cat.geometry.centroid.coords[0][0]+10:.0f}",
                py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
                boundary_node=cat.boundary_node,
            )

paved_from_input(catchments: ExtendedGeoDataFrame, landuse: StrictStr | Path, surface_level: StrictStr | Path, street_storage: StrictStr | Path | float, sewer_storage: StrictStr | Path | float, pump_capacity: StrictStr | Path | float, meteo_areas: ExtendedGeoDataFrame, overflows: ExtendedGeoDataFrame = None, sewer_areas: ExtendedGeoDataFrame = None, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction') -> None

Generate contents of RR paved nodes

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	catchment areas	required
landuse	str	filename of landuse raster	required
surface_level	str	file name of suface level raster	required
street_storage	Union	numeric for spatially uniform street storage (mm), or raster for distributed values	required
sewer_storage	Union	numeric for spatially uniform sewer storage (mm), or raster for distributed values	required
pump_capacity	Union	numeric for spatially uniform pump capaities (mm), or raster for distributed values	required
meteo_areas	ExtendedGeoDataFrame	meteo areas, for each station a meteo time series is assigned	required
overflows	ExtendedGeoDataFrame	overflow locations. Defaults to None.	None
sewer_areas	ExtendedGeoDataFrame	sewer area locations. Defaults to None.	None
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Returns:

Name	Type	Description
_type_	None	description

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def paved_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    landuse: StrictStr | Path,
    surface_level: StrictStr | Path,
    street_storage:StrictStr | Path | float,
    sewer_storage:StrictStr | Path | float,
    pump_capacity:StrictStr | Path | float,
    meteo_areas: ExtendedGeoDataFrame,
    overflows: ExtendedGeoDataFrame = None,
    sewer_areas: ExtendedGeoDataFrame = None,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
) -> None:
    """Generate contents of RR paved nodes

    Args:
        catchments (ExtendedGeoDataFrame): catchment areas
        landuse (str): filename of landuse raster
        surface_level (str): file name of suface level raster
        street_storage (Union): numeric for spatially uniform street storage (mm), or raster for distributed values
        sewer_storage (Union): numeric for spatially uniform sewer storage (mm), or raster for distributed values
        pump_capacity (Union): numeric for spatially uniform pump capaities (mm), or raster for distributed values
        meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
        overflows (ExtendedGeoDataFrame, optional): overflow locations. Defaults to None.
        sewer_areas (ExtendedGeoDataFrame, optional): sewer area locations. Defaults to None.
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".

    Returns:
        _type_: _description_
    """
    px_area, reference_crs = _raster_metadata(landuse)

    lu_counts = self.paved.drrmodel.zonal_cache.zonal_category_counts(
        gpd.GeoDataFrame(catchments),
        landuse,
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )
    stats_elev = self.paved.drrmodel.zonal_cache.zonal_stats(
        gpd.GeoDataFrame(catchments),
        surface_level,
        statistics=("median",),
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )

    if isinstance(street_storage, (Path, str)):
        str_stors = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            street_storage,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )

    if isinstance(sewer_storage, (Path, str)):
        sew_stors = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            sewer_storage,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )

    if isinstance(pump_capacity,  (Path, str)):
        # raster of POC in mm/h
        pump_caps = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            pump_capacity,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )
    paved_columns = [
        "id",
        "area",
        "surface_level",
        "street_storage",
        "sewer_storage",
        "pump_capacity",
        "meteo_area",
        "px",
        "py",
        "boundary_node",
    ]
    if sewer_areas is not None:
        # if the parameters ara rasters, do the zonal statistics per sewage area as well.
        if isinstance(street_storage,(Path, str)):
            str_stors_sa = zonal.zonal_stats(
                gpd.GeoDataFrame(sewer_areas),
                street_storage,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        if isinstance(sewer_storage, (Path, str)):
            sew_stors_sa = zonal.zonal_stats(
                gpd.GeoDataFrame(sewer_areas),
                sewer_storage,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        if isinstance(pump_capacity, (Path, str)):
            pump_caps_sa = zonal.zonal_stats(
                gpd.GeoDataFrame(sewer_areas),
                pump_capacity,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        stats_sa_elev = zonal.zonal_stats(
            gpd.GeoDataFrame(sewer_areas),
            surface_level,
            statistics=("median",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )
        sewer_lu_counts = zonal.zonal_category_counts(
            gpd.GeoDataFrame(sewer_areas),
            landuse,
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        # find the paved area in the sewer areas
        for isew, sew in enumerate(sewer_areas.itertuples()):
            pav_area = 0
            pixels = sewer_lu_counts.iloc[isew]
            pav_pixels = pixels.get(14, 0)
            if pav_pixels == 0:
                logger.warning("No paved area in sewer area %s.", sew.code)
                self.paved.add_paved(**dict.fromkeys(paved_columns, "0.0"))
                continue
            pav_area += pav_pixels * px_area

            # subtract it fromthe total paved area in this catchment, make sure at least 0 remains
            # lu_counts[cat_ind][14.0] -=  pav_pixels
            # if lu_counts[cat_ind][14.0] < 0: lu_counts[cat_ind][14.0]  = 0

            elev = stats_sa_elev.iloc[isew]["median"]
            # find overflows related to this sewer area
            overflow_sa = RelationsRR.OVERFLOW_SEWER_AREA
            sewer_id = getattr(sew, overflow_sa.parent_key)
            ovf = overflows[overflows[overflow_sa.foreign_key] == sewer_id]
            for ov in ovf.itertuples():
                # find corresponding meteo-station
                tm = [
                    m
                    for m in meteo_areas.itertuples()
                    if m.geometry.contains(sew.geometry.centroid)
                ]
                ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

                # if a float is given, a standard value is passed. If a string is given, a rastername is assumed to zonal statistics are applied.
                if isinstance(street_storage, float):
                    street_storage_val = f"{street_storage:.2f}"
                elif isinstance(street_storage, (Path, str)):
                    street_storage_val = f'{str_stors_sa.iloc[isew]["mean"]:.2f}'
                else:
                    raise TypeError('Street_storage has the wrong datatype. It should be a filename (Path or string) or number (float or int).')

                # three options: it can be an attribute of a sewer area, a uniform value or a raster
                riool_berging_mm = getattr(sew, "riool_berging_mm", None)
                if riool_berging_mm is None or np.isnan(riool_berging_mm) or not isinstance(riool_berging_mm, float):
                    if isinstance(sewer_storage, float):
                        sewer_storage_val = f"{sewer_storage:.2f}"
                    elif isinstance(sewer_storage, (Path, str)):
                        sewer_storage_val = f'{sew_stors_sa.iloc[isew]["mean"]:.2f}'
                    else:
                        raise TypeError('Sewer_storage has the wrong datatype. It should be a filename (Path or string) or number (float or int).')
                else:
                    sewer_storage_val = f'{riool_berging_mm:.2f}'

                # three options: it can be an attribute of a sewer area, a uniform value or a raster
                riool_poc_m3s = getattr(sew, "riool_poc_m3s", None)
                if riool_poc_m3s is None or np.isnan(riool_poc_m3s) or not isinstance(riool_poc_m3s, float):
                    if isinstance(pump_capacity, float):
                        # convert the value from mm/h to m3/s
                        pump_capacity_val = f"{pump_capacity * (float(pav_area) * ov.fractie) / (1000. * 3600.):.8f}"
                    elif isinstance(pump_capacity, (Path, str)):
                        # convert the value (extracted from the raster) from mm/h to m3/s
                        pump_capacity_val = f'{pump_caps_sa.iloc[isew]["mean"] * (float(pav_area) * ov.fractie) / (1000. * 3600.):.8f}'
                    else:
                        raise TypeError('Pump_capacity has the wrong datatype. It should be a filename (Path or string) or number (float or int).')
                else:
                    # use the attribute value
                    pump_capacity_val = f'{riool_poc_m3s * ov.fractie:.8f}'

                # add prefix to the overflow id to create the paved-node id
                self.paved.add_paved(
                    id=str(ov.code),
                    area=str(pav_area * ov.fractie),
                    surface_level=f"{elev:.2f}",
                    street_storage=street_storage_val,
                    sewer_storage=sewer_storage_val,
                    pump_capacity=pump_capacity_val,
                    meteo_area=str(ms),
                    px=f"{ov.geometry.coords[0][0]+10:.0f}",
                    py=f"{ov.geometry.coords[0][1]:.0f}",
                    boundary_node=ov.code,
                )

    if sewer_areas is not None:
        # Land use of the catchment parts outside the sewer areas, by
        # catchment position. Catchments that overlap no sewer area use
        # lu_counts.
        sewer_union = sewer_areas.union_all()
        geometries = catchments.geometry.reset_index(drop=True)
        overlaps = geometries.intersects(sewer_union)
        outside = geometries[overlaps].difference(sewer_union)
        outside = outside[outside.area > 0]
        outside_counts = pd.DataFrame()
        if not outside.empty:
            outside_counts = zonal.zonal_category_counts(
                outside,
                landuse,
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
            )

    for num, cat in enumerate(catchments.itertuples()):
        # if no rasterdata could be obtained for this catchment, skip it.
        if pd.isna(stats_elev.iloc[num]["median"]):
            logger.warning(NO_RASTERDATA_WARNING, cat.code)
            self.paved.add_paved(**dict.fromkeys(paved_columns, "0.0"))
            continue
        if sewer_areas is None or not overlaps[num]:
            pav_area = str(lu_counts.iloc[num].get(14, 0) * px_area)
        elif num in outside_counts.index:
            pav_area = str(outside_counts.loc[num].get(14, 0) * px_area)
        else:
            logger.info(
                f"No paved area outside sewer area in catchments {cat.code}."
            )
            pav_area = 0.0

        # find corresponding meteo-station
        tm = [
            m
            for m in meteo_areas.itertuples()
            if m.geometry.contains(cat.geometry.centroid)
        ]
        ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code

        elev = stats_elev.iloc[num]["median"]
        # if a float is given, a standard value is passed. If a string is given, a rastername is assumed to zonal statistics are applied.
        street_storage_val = (
            f"{street_storage:.2f}"
            if isinstance(street_storage, float)
            else f'{str_stors.iloc[num]["mean"]:.2f}'
        )
        sewer_storage_val = (
            f"{sewer_storage:.2f}"
            if isinstance(sewer_storage, float)
            else f'{sew_stors.iloc[num]["mean"]:.2f}'
        )
        pump_capacity_val = (
            f'{(pump_capacity * float(pav_area)) / (1000. * 3600.):.8f}'
            if isinstance(pump_capacity, float)
            else f'{pump_caps.iloc[num]["mean"] * (float(pav_area)) / (1000. * 3600.):.8f}'
        )
        self.paved.add_paved(
            id=str(cat.code),
            area=str(pav_area),
            surface_level=f"{elev:.2f}",
            street_storage=street_storage_val,
            sewer_storage=sewer_storage_val,
            pump_capacity=pump_capacity_val,
            meteo_area=str(ms),
            px=f"{cat.geometry.centroid.coords[0][0]+10:.0f}",
            py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
            boundary_node=cat.boundary_node,
        )

UnpavedIO

Source code in hydrolib/dhydamo/io/drrreader.py

class UnpavedIO:
    def __init__(self, unpaved):
        self.unpaved = unpaved

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def unpaved_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        landuse: StrictStr | Path,
        surface_level: StrictStr | Path,
        soiltype: StrictStr | Path,
        surface_storage: StrictStr | Path | float,
        infiltration_capacity: StrictStr | Path | float,
        initial_gwd: StrictStr | Path | float,
        meteo_areas: ExtendedGeoDataFrame,
        greenhouse_areas: ExtendedGeoDataFrame = None,
        *,
        min_coverage_frac: float = 0.0,
        coverage_weight: zonal.CoverageWeight = "fraction",
    ):
        """Generate contents of an unpaved node from raster data

        Args:
            catchments (ExtendedGeoDataFrame): catchment areas
            landuse (str): filename of land use raster
            surface_level (str): file name of surface level raster
            soiltype (str): file name of soiltype raster
            surface_storage (Union): numeric for spatially uniform surface storage (mm), or raster for distributed values
            infiltration_capacity (Union): numeric for spatially uniform infiltration capacity (mm/d), or raster for distributed values
            initial_gwd (Union): numeric for spatially uniform initial groundwater levels (m), or raster for distributed values
            meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
            min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
            coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
        """
        px_area, reference_crs = _raster_metadata(landuse)

        # required rasters
        warnings.filterwarnings("ignore")
        lu_counts = self.unpaved.drrmodel.zonal_cache.zonal_category_counts(
            gpd.GeoDataFrame(catchments),
            landuse,
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        soiltypes = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            soiltype,
            statistics=("mode",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        stats_elev = self.unpaved.drrmodel.zonal_cache.zonal_stats(
            gpd.GeoDataFrame(catchments),
            surface_level,
            statistics=("median",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
        )

        # optional rasters
        if isinstance(surface_storage, (Path, str)):
            sstores = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                surface_storage,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        elif isinstance(surface_storage, int):
            surface_storage = float(surface_storage)
        if isinstance(infiltration_capacity, (Path, str)):
            infcaps = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                infiltration_capacity,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        elif isinstance(infiltration_capacity, int):
            infiltration_capacity = float(infiltration_capacity)
        if isinstance(initial_gwd, (Path, str)):
            ini_gwds = zonal.zonal_stats(
                gpd.GeoDataFrame(catchments),
                initial_gwd,
                statistics=("mean",),
                min_coverage_frac=min_coverage_frac,
                coverage_weight=coverage_weight,
                strategy="feature-sequential",
                raster_crs=reference_crs,
            )
        elif isinstance(initial_gwd, int):
            initial_gwd = float(initial_gwd)

        # HyDAMO Crop code; hydamo name, sobek index, sobek name:
        # 1 aardappelen   3 potatoes
        # 2 graan         5 grain
        # 3 suikerbiet    4 sugarbeet
        # 4 mais          2 corn
        # 5 overige gew. 15 vegetables
        # 6 bloembollen  10 bulbous plants
        # 7 boomgaard     9 orchard
        # 8 gras          1 grass
        # 9 loofbos      11 dediduous
        # 10 naaldbos    12 conferous
        # 11 natuuur     13 nature
        # 12 braak       14 fallow
        sobek_indices = [3, 5, 4, 2, 15, 10, 9, 1, 11, 12, 13, 14]
        for num, cat in enumerate(catchments.itertuples()):
            # if no rasterdata could be obtained for this catchment, skip it.
            if pd.isna(stats_elev.iloc[num]["median"]):
                logger.warning(NO_RASTERDATA_WARNING, cat.code)
                self.unpaved.add_unpaved(
                    id="0.0",
                    total_area="0.0",
                    lu_areas="0.0",
                    surface_level="0.0",
                    soiltype="0.0",
                    surface_storage="0.0",
                    infiltration_capacity="0.0",
                    initial_gwd="0.0",
                    meteo_area="0.0",
                    px="0.0",
                    py="0.0",
                    boundary_node="0.0",
                )
                continue
            tm = [
                m
                for m in meteo_areas.itertuples()
                if m.geometry.contains(cat.geometry.centroid)
            ]
            ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code
            landuse_counts = lu_counts.iloc[num].copy()
            mapping = np.zeros(16, dtype=int)

            # subtract greenhouse area from most occurring land use if no greenhouse area is in the landuse map
            if greenhouse_areas is not None and cat.geometry.intersects(greenhouse_areas.geometry).any():
                intersection_area = cat.geometry.intersection(greenhouse_areas.geometry).area
                intersection_area = intersection_area[intersection_area > 0.].to_numpy()[0]                    
                if landuse_counts.get(15, 0) > 0:
                    # divide area to subtract between greenhouses and the most occurring area
                    remainder = np.max([0., intersection_area - float(landuse_counts[15] * px_area)])
                else:    
                    remainder = intersection_area
                maxind = landuse_counts.idxmax()
                logger.info(
                    "Catchment %s: subtracting %s m2 from class %s for supplied greenhouse area.",
                    cat.code,
                    remainder,
                    maxind,
                )
                landuse_counts[maxind] = np.max(
                    [0.0, landuse_counts[maxind] - np.round(remainder / px_area)]
                )

            for i in range(1, 13):
                mapping[sobek_indices[i - 1] - 1] = landuse_counts.get(i, 0) * px_area
            lu_map = " ".join(map(str, mapping))
            elev = stats_elev.iloc[num]["median"]
            self.unpaved.add_unpaved(
                id=str(cat.code),
                total_area=f"{cat.geometry.area:.0f}",
                lu_areas=lu_map,
                surface_level=f"{elev:.2f}",
                soiltype=f'{soiltypes.iloc[num]["mode"] + 100.0:.0f}',
                surface_storage=(
                    f"{surface_storage:.3f}"
                    if isinstance(surface_storage, float)
                    else f'{sstores.iloc[num]["mean"]:.3f}'
                ),
                infiltration_capacity=(
                    f"{infiltration_capacity:.3f}"
                    if isinstance(infiltration_capacity, float)
                    else f'{infcaps.iloc[num]["mean"]:.3f}'
                ),
                initial_gwd=(
                    f"{initial_gwd:.2f}"
                    if isinstance(initial_gwd, float)
                    else f'{ini_gwds.iloc[num]["mean"]:.2f}'
                ),
                meteo_area=str(ms),
                px=f"{cat.geometry.centroid.coords[0][0]-10:.0f}",
                py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
                boundary_node=cat.boundary_node,
            )

    @validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
    def ernst_from_input(
        self,
        catchments: ExtendedGeoDataFrame,
        depths: list,
        resistance: list,
        infiltration_resistance: float | None = None,
        runoff_resistance: float | None = None,
    ) -> None:
        """Generate an Ernst definition for an unpaved node.

        Args:
            catchments (ExtendedGeoDataFrame): Cahchment areas
            depths (list): list of layer depths (m)
            resistance (list): list of layer Ernst resistances (d-1)
            infiltration_resistance (int or float, optional): resistance to infiltration. Defaults to 300 d-1.
            runoff_resistance (int or float, optional): resistance to suface runoff. Defaults to 1 d-1.
        """
        if infiltration_resistance is None:
            infiltration_resistance = 300.0
        if runoff_resistance is None:
            runoff_resistance = 1.0

        cvo = " ".join([str(res) for res in resistance])
        lv = " ".join([str(depth) for depth in depths])
        cvi = f"{infiltration_resistance:.2f}"
        cvs = f"{runoff_resistance:.2f}"
        for cat in catchments.itertuples():
            self.unpaved.add_ernst_def(
                id=str(cat.code), cvo=cvo, lv=lv, cvi=cvi, cvs=cvs
            )

ernst_from_input(catchments: ExtendedGeoDataFrame, depths: list, resistance: list, infiltration_resistance: float | None = None, runoff_resistance: float | None = None) -> None

Generate an Ernst definition for an unpaved node.

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	Cahchment areas	required
depths	list	list of layer depths (m)	required
resistance	list	list of layer Ernst resistances (d-1)	required
infiltration_resistance	int or float	resistance to infiltration. Defaults to 300 d-1.	None
runoff_resistance	int or float	resistance to suface runoff. Defaults to 1 d-1.	None

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def ernst_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    depths: list,
    resistance: list,
    infiltration_resistance: float | None = None,
    runoff_resistance: float | None = None,
) -> None:
    """Generate an Ernst definition for an unpaved node.

    Args:
        catchments (ExtendedGeoDataFrame): Cahchment areas
        depths (list): list of layer depths (m)
        resistance (list): list of layer Ernst resistances (d-1)
        infiltration_resistance (int or float, optional): resistance to infiltration. Defaults to 300 d-1.
        runoff_resistance (int or float, optional): resistance to suface runoff. Defaults to 1 d-1.
    """
    if infiltration_resistance is None:
        infiltration_resistance = 300.0
    if runoff_resistance is None:
        runoff_resistance = 1.0

    cvo = " ".join([str(res) for res in resistance])
    lv = " ".join([str(depth) for depth in depths])
    cvi = f"{infiltration_resistance:.2f}"
    cvs = f"{runoff_resistance:.2f}"
    for cat in catchments.itertuples():
        self.unpaved.add_ernst_def(
            id=str(cat.code), cvo=cvo, lv=lv, cvi=cvi, cvs=cvs
        )

unpaved_from_input(catchments: ExtendedGeoDataFrame, landuse: StrictStr | Path, surface_level: StrictStr | Path, soiltype: StrictStr | Path, surface_storage: StrictStr | Path | float, infiltration_capacity: StrictStr | Path | float, initial_gwd: StrictStr | Path | float, meteo_areas: ExtendedGeoDataFrame, greenhouse_areas: ExtendedGeoDataFrame = None, *, min_coverage_frac: float = 0.0, coverage_weight: zonal.CoverageWeight = 'fraction')

Generate contents of an unpaved node from raster data

Parameters:

Name	Type	Description	Default
catchments	ExtendedGeoDataFrame	catchment areas	required
landuse	str	filename of land use raster	required
surface_level	str	file name of surface level raster	required
soiltype	str	file name of soiltype raster	required
surface_storage	Union	numeric for spatially uniform surface storage (mm), or raster for distributed values	required
infiltration_capacity	Union	numeric for spatially uniform infiltration capacity (mm/d), or raster for distributed values	required
initial_gwd	Union	numeric for spatially uniform initial groundwater levels (m), or raster for distributed values	required
meteo_areas	ExtendedGeoDataFrame	meteo areas, for each station a meteo time series is assigned	required
min_coverage_frac	float	smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.	0.0
coverage_weight	str	"fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".	'fraction'

Source code in hydrolib/dhydamo/io/drrreader.py

@validate_arguments(config=ConfigDict(arbitrary_types_allowed=True))
def unpaved_from_input(
    self,
    catchments: ExtendedGeoDataFrame,
    landuse: StrictStr | Path,
    surface_level: StrictStr | Path,
    soiltype: StrictStr | Path,
    surface_storage: StrictStr | Path | float,
    infiltration_capacity: StrictStr | Path | float,
    initial_gwd: StrictStr | Path | float,
    meteo_areas: ExtendedGeoDataFrame,
    greenhouse_areas: ExtendedGeoDataFrame = None,
    *,
    min_coverage_frac: float = 0.0,
    coverage_weight: zonal.CoverageWeight = "fraction",
):
    """Generate contents of an unpaved node from raster data

    Args:
        catchments (ExtendedGeoDataFrame): catchment areas
        landuse (str): filename of land use raster
        surface_level (str): file name of surface level raster
        soiltype (str): file name of soiltype raster
        surface_storage (Union): numeric for spatially uniform surface storage (mm), or raster for distributed values
        infiltration_capacity (Union): numeric for spatially uniform infiltration capacity (mm/d), or raster for distributed values
        initial_gwd (Union): numeric for spatially uniform initial groundwater levels (m), or raster for distributed values
        meteo_areas (ExtendedGeoDataFrame): meteo areas, for each station a meteo time series is assigned
        min_coverage_frac (float, optional): smallest fraction of a pixel a zone must cover for the pixel to count in the zonal statistics. Defaults to 0, every touched pixel.
        coverage_weight (str, optional): "fraction" counts pixels by the fraction of them the zone covers, "none" counts them fully. Defaults to "fraction".
    """
    px_area, reference_crs = _raster_metadata(landuse)

    # required rasters
    warnings.filterwarnings("ignore")
    lu_counts = self.unpaved.drrmodel.zonal_cache.zonal_category_counts(
        gpd.GeoDataFrame(catchments),
        landuse,
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )

    soiltypes = zonal.zonal_stats(
        gpd.GeoDataFrame(catchments),
        soiltype,
        statistics=("mode",),
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )

    stats_elev = self.unpaved.drrmodel.zonal_cache.zonal_stats(
        gpd.GeoDataFrame(catchments),
        surface_level,
        statistics=("median",),
        min_coverage_frac=min_coverage_frac,
        coverage_weight=coverage_weight,
        strategy="feature-sequential",
    )

    # optional rasters
    if isinstance(surface_storage, (Path, str)):
        sstores = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            surface_storage,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )
    elif isinstance(surface_storage, int):
        surface_storage = float(surface_storage)
    if isinstance(infiltration_capacity, (Path, str)):
        infcaps = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            infiltration_capacity,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )
    elif isinstance(infiltration_capacity, int):
        infiltration_capacity = float(infiltration_capacity)
    if isinstance(initial_gwd, (Path, str)):
        ini_gwds = zonal.zonal_stats(
            gpd.GeoDataFrame(catchments),
            initial_gwd,
            statistics=("mean",),
            min_coverage_frac=min_coverage_frac,
            coverage_weight=coverage_weight,
            strategy="feature-sequential",
            raster_crs=reference_crs,
        )
    elif isinstance(initial_gwd, int):
        initial_gwd = float(initial_gwd)

    # HyDAMO Crop code; hydamo name, sobek index, sobek name:
    # 1 aardappelen   3 potatoes
    # 2 graan         5 grain
    # 3 suikerbiet    4 sugarbeet
    # 4 mais          2 corn
    # 5 overige gew. 15 vegetables
    # 6 bloembollen  10 bulbous plants
    # 7 boomgaard     9 orchard
    # 8 gras          1 grass
    # 9 loofbos      11 dediduous
    # 10 naaldbos    12 conferous
    # 11 natuuur     13 nature
    # 12 braak       14 fallow
    sobek_indices = [3, 5, 4, 2, 15, 10, 9, 1, 11, 12, 13, 14]
    for num, cat in enumerate(catchments.itertuples()):
        # if no rasterdata could be obtained for this catchment, skip it.
        if pd.isna(stats_elev.iloc[num]["median"]):
            logger.warning(NO_RASTERDATA_WARNING, cat.code)
            self.unpaved.add_unpaved(
                id="0.0",
                total_area="0.0",
                lu_areas="0.0",
                surface_level="0.0",
                soiltype="0.0",
                surface_storage="0.0",
                infiltration_capacity="0.0",
                initial_gwd="0.0",
                meteo_area="0.0",
                px="0.0",
                py="0.0",
                boundary_node="0.0",
            )
            continue
        tm = [
            m
            for m in meteo_areas.itertuples()
            if m.geometry.contains(cat.geometry.centroid)
        ]
        ms = meteo_areas.iloc[0, 0] if not tm else tm[0].code
        landuse_counts = lu_counts.iloc[num].copy()
        mapping = np.zeros(16, dtype=int)

        # subtract greenhouse area from most occurring land use if no greenhouse area is in the landuse map
        if greenhouse_areas is not None and cat.geometry.intersects(greenhouse_areas.geometry).any():
            intersection_area = cat.geometry.intersection(greenhouse_areas.geometry).area
            intersection_area = intersection_area[intersection_area > 0.].to_numpy()[0]                    
            if landuse_counts.get(15, 0) > 0:
                # divide area to subtract between greenhouses and the most occurring area
                remainder = np.max([0., intersection_area - float(landuse_counts[15] * px_area)])
            else:    
                remainder = intersection_area
            maxind = landuse_counts.idxmax()
            logger.info(
                "Catchment %s: subtracting %s m2 from class %s for supplied greenhouse area.",
                cat.code,
                remainder,
                maxind,
            )
            landuse_counts[maxind] = np.max(
                [0.0, landuse_counts[maxind] - np.round(remainder / px_area)]
            )

        for i in range(1, 13):
            mapping[sobek_indices[i - 1] - 1] = landuse_counts.get(i, 0) * px_area
        lu_map = " ".join(map(str, mapping))
        elev = stats_elev.iloc[num]["median"]
        self.unpaved.add_unpaved(
            id=str(cat.code),
            total_area=f"{cat.geometry.area:.0f}",
            lu_areas=lu_map,
            surface_level=f"{elev:.2f}",
            soiltype=f'{soiltypes.iloc[num]["mode"] + 100.0:.0f}',
            surface_storage=(
                f"{surface_storage:.3f}"
                if isinstance(surface_storage, float)
                else f'{sstores.iloc[num]["mean"]:.3f}'
            ),
            infiltration_capacity=(
                f"{infiltration_capacity:.3f}"
                if isinstance(infiltration_capacity, float)
                else f'{infcaps.iloc[num]["mean"]:.3f}'
            ),
            initial_gwd=(
                f"{initial_gwd:.2f}"
                if isinstance(initial_gwd, float)
                else f'{ini_gwds.iloc[num]["mean"]:.2f}'
            ),
            meteo_area=str(ms),
            px=f"{cat.geometry.centroid.coords[0][0]-10:.0f}",
            py=f"{cat.geometry.centroid.coords[0][1]:.0f}",
            boundary_node=cat.boundary_node,
        )

drrwriter

DRRWriter

Writer for RR files

Source code in hydrolib/dhydamo/io/drrwriter.py

class DRRWriter:
    """Writer for RR files"""

    def __init__(self, rrmodel, output_dir, name=None, wwtp=None):
        self.rrmodel = rrmodel
        # self.geometries = parent.geometries
        # self.boundary_conditions = parent.boundary_conditions
        self.output_dir = os.path.join(output_dir, "rr")

        if name is not None:
            logger.warning(
                "The 'name' parameter is deprecated and will be removed in a future version.",                
        )

        hydamo = HyDAMO()
        self.version = hydamo.version

        if os.path.exists(self.output_dir):
            shutil.rmtree(self.output_dir)

        self.d3b_parameters = rrmodel.d3b_parameters

        self.precip_df = {}
        self.run_dimrpad = rrmodel.dimr_path

        if wwtp is None:
            self.wwtp = Point((1e5, 5e5))
        else:
            self.wwtp = Point((wwtp[0], wwtp[1]))

    def write_all(self):  # write all RR files
        """
        Wrapper method to write all components
        """
        # copy all 'standard' files from the resources folder
        self.copyRRFiles()

        self.write_topology()

        self.write_unpaved()

        self.write_paved()

        self.write_greenhouse()

        self.write_openwater()

        self.write_meteo()

        # Change d3b parameters
        for parameter, value in self.d3b_parameters.items():
            self.change_d3b_parameter(parameter, value)

    def copyRRFiles(self):
        """
        Method to copy the RR-files from resources.

        Returns
        -------
        bool
            Succesfull or not.

        """
        srcRR = os.path.join(os.path.dirname(__file__), "..", "resources", "RR")
        targetRR = os.path.join(self.output_dir)
        shutil.copytree(srcRR, targetRR)
        return True

    def write_topology(self):
        """
        Wrapper to write the topolgy files for RR. The following files are written:
            3B_NOD.TP: IDs and coordinates of all nodes: unpaved, paved, greenhouse, open water and boundary. Only nodes with a total area of >0 are written. Nodes are placed up to 20 m to the east and west of the catchment centroid.
            3B_LINK.TP: Links between all nodes and the correct boundaries.
            BOUND3B.3B: coupling information for every RR boundary.
            BoundaryConditions.bc: see bound3b.3b
        Returns
        -------
        None.

        """
        filepath = os.path.join(self.output_dir, "3B_NOD.TP")
        with open(filepath, "w") as f:

            # Unpaved nodes
            if any(self.rrmodel.unpaved.unp_nodes):
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        f.write(
                            "NODE id '"
                            + dct["id"]
                            + "' nm '"
                            + dct["id"]
                            + "' ri '-1' mt 1 '2' nt 44 ObID '3B_UNPAVED' px "
                            + dct["px"]
                            + " py "
                            + dct["py"]
                            + " node\n"
                        )

            # Paved nodes
            if any(self.rrmodel.paved.pav_nodes):
                for _, dct in self.rrmodel.paved.pav_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "NODE id '"
                            + dct["id"]
                            + "' nm '"
                            + dct["id"]
                            + "' ri '-1' mt 1 '1' nt 43 ObID '3B_PAVED' px "
                            + dct["px"]
                            + " py "
                            + dct["py"]
                            + " node\n"
                        )

                f.write(
                    "NODE id 'WWTP' nm 'WWTP' ri '-1' mt 1 '14' nt 56 ObID '3B_WWTP' px "
                    + str(self.wwtp.coords[0][0])
                    + " py "
                    + str(self.wwtp.coords[0][1])
                    + " node\n"
                )
                f.write(
                    "NODE id 'WWTP_BND' nm 'WWTP_BND' ri '-1' mt 1 '6' nt 47 ObID '3B_BOUNDARY' px "
                    + str(self.wwtp.coords[0][0] + 50.0)
                    + " py "
                    + str(self.wwtp.coords[0][1] + 50.0)
                    + " node\n"
                )

            # Greenhouse nodes
            if any(self.rrmodel.greenhouse.gh_nodes):
                for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "NODE id '"
                            + dct["id"]
                            + "' nm '"
                            + dct["id"]
                            + "' ri '-1' mt 1 '3' nt 45 ObID '3B_GREENHOUSE' px "
                            + dct["px"]
                            + " py "
                            + dct["py"]
                            + " node\n"
                        )

            # Openwater nodes
            if any(self.rrmodel.openwater.ow_nodes):
                for _, dct in self.rrmodel.openwater.ow_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "NODE id '"
                            + dct["id"]
                            + "' nm '"
                            + dct["id"]
                            + "' ri '-1' mt 1 '21' nt 46 ObID 'OW_PRECIP' px "
                            + dct["px"]
                            + " py "
                            + dct["py"]
                            + " node\n"
                        )

            # Boundary_nodes
            if any(self.rrmodel.external_forcings.boundary_nodes):
                for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
                    f.write(
                        "NODE id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' ri '-1' mt 1 '6' nt 78 ObID 'SBK_SBK-3B-NODE' px "
                        + dct["px"]
                        + " py "
                        + dct["py"]
                        + " node\n"
                    )

        filepath = os.path.join(self.output_dir, "3B_LINK.TP")
        with open(filepath, "w") as f:

            cnt = 0
            if any(self.rrmodel.unpaved.unp_nodes):
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        cnt += 1
                        f.write(
                            "BRCH id 'link_"
                            + str(cnt)
                            + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                            + dct["id"]
                            + "' en '"
                            + dct["boundary_node"]
                            + "' brch\n"
                        )
            if any(self.rrmodel.paved.pav_nodes):
                for _, dct in self.rrmodel.paved.pav_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        cnt += 1
                        f.write(
                            "BRCH id 'link_"
                            + str(cnt)
                            + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                            + dct["id"]
                            + "' en '"
                            + dct["boundary_node"]
                            + "' brch\n"
                        )
                        cnt += 1
                        f.write(
                            "BRCH id 'link_"
                            + str(cnt)
                            + "' ri '-1' mt 1 '1' bt 18 ObID '3B_LINK_RWZI' bn '"
                            + dct["id"]
                            + "' en 'WWTP' brch\n"
                        )
                cnt += 1
                f.write(
                    "BRCH id 'link_"
                    + str(cnt)
                    + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn 'WWTP' en 'WWTP_BND' brch\n"
                )

            if any(self.rrmodel.greenhouse.gh_nodes):
                for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        cnt += 1
                        f.write(
                            "BRCH id 'link_"
                            + str(cnt)
                            + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                            + dct["id"]
                            + "' en '"
                            + dct["boundary_node"]
                            + "' brch\n"
                        )
            if any(self.rrmodel.openwater.ow_nodes):
                for _, dct in self.rrmodel.openwater.ow_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        cnt += 1
                        f.write(
                            "BRCH id 'link_"
                            + str(cnt)
                            + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                            + dct["id"]
                            + "' en '"
                            + dct["boundary_node"]
                            + "' brch\n"
                        )

        # bound3b.3b
        filepath = os.path.join(self.output_dir, "Bound3B.3B")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
                f.write("BOUN id '" + dct["id"] + "' bl 2 '0' is 0 boun\n")
            if any(self.rrmodel.paved.pav_nodes):
                f.write("BOUN id 'WWTP_BND' bl 0 -3 is 0 boun\n")

        # BoundaryConditions.bc
        filepath = os.path.join(self.output_dir, "BoundaryConditions.bc")
        header = {"fileVersion": "1.01", "fileType": "boundConds"}
        with open(filepath, "w") as f:
            self._write_dict(f, header, "General", "\n")
            for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
                temp = {
                    "name": "" + dct["id"],
                    "function": "constant",
                    "quantity": "water_level",
                    "unit": "m",
                }
                self._write_dict(f, temp, "Boundary", "    0\n\n")
            if any(self.rrmodel.paved.pav_nodes):
                temp = {
                    "name": "WWTP_BND",
                    "function": "constant",
                    "quantity": "water_level",
                    "unit": "m",
                }
                self._write_dict(f, temp, "Boundary", "    0\n\n")

    def write_unpaved(self):
        """
        Method to write all files associated with unpaved nodes: UNPAVED.3B, UNPAVED.ALF, UNPAVED.STO, UNPAVED.INF and UNPAVED.SEP. All files contain a definition for every node  because they may or may not be spatially distributed.

        Returns
        -------
        None.

        """
        if any(self.rrmodel.unpaved.unp_nodes):
            filepath = os.path.join(self.output_dir, "UNPAVED.3B")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        f.write(
                            "UNPV id '"
                            + dct["id"]
                            + "' na "
                            + dct["na"]
                            + " ga "
                            + dct["ga"]
                            + " ar "
                            + dct["ar"]
                            + " lv "
                            + dct["lv"]
                            + " co "
                            + dct["co"]
                            + " su "
                            + dct["su"]
                            + " sd 'sto_"
                            + dct["id"]
                            + "' ic 'inf_"
                            + dct["id"]
                            + "' bt "
                            + dct["bt"]
                            + " ed '"
                            + dct["ed"]
                            + "' sp '"
                            + dct["sp"]
                            + "' ig 0 "
                            + dct["ig"]
                            + " mg "
                            + dct["mg"]
                            + " gl "
                            + dct["gl"]
                            + " is "
                            + dct["is"]
                            + " ms '"
                            + dct["ms"]
                            + "' unpv\n"
                        )

            filepath = os.path.join(self.output_dir, "UNPAVED.STO")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        f.write(
                            "STDF id 'sto_"
                            + dct["id"]
                            + "' nm 'sto_"
                            + dct["id"]
                            + "' ml "
                            + dct["sd"]
                            + " il 0 stdf\n"
                        )

            filepath = os.path.join(self.output_dir, "UNPAVED.INF")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        f.write(
                            "INFC id 'inf_"
                            + dct["id"]
                            + "' nm 'inf_"
                            + dct["id"]
                            + "' ic "
                            + dct["ic"]
                            + " infc\n"
                        )

            filepath = os.path.join(self.output_dir, "UNPAVED.ALF")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.unpaved.ernst_defs.items():
                    if (
                        np.sum(
                            [
                                float(d)
                                for d in self.rrmodel.unpaved.unp_nodes[_]["ar"].split(
                                    " "
                                )
                            ]
                        )
                        > 0.0
                    ):
                        f.write(
                            "ERNS id '"
                            + dct["id"]
                            + "' nm '"
                            + dct["id"]
                            + "' cvi "
                            + dct["cvi"]
                            + " cvo 0 "
                            + dct["cvo"]
                            + " cvs "
                            + dct["cvs"]
                            + " lv "
                            + dct["lv"]
                            + " erns\n"
                        )

            filepath = os.path.join(self.output_dir, "UNPAVED.SEP")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                    if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                        f.write(
                            "SEEP id '"
                            + dct["sp"]
                            + "' nm '"
                            + dct["sp"]
                            + "' co 4 PDIN 0 0 pdin ss 0\n"
                        )
                        f.write("TBLE\n")
                        for row in self.rrmodel.external_forcings.seepage[dct["sp"]][
                            "seepage"
                        ].items():
                            f.write(
                                "'"
                                + row[0].strftime("%Y/%m/%d;%H:%M:%S")
                                + "' "
                                + f"{row[1]:.5f} <\n"
                            )
                        f.write("tble\nseep\n")

    def write_paved(self):
        """
        Method to write all files associated with paved nodes: PAVED.3B, PAVED.STO and PAVED.DWA. The latter contains only one definition.

        Returns
        -------
        None.

        """
        if any(self.rrmodel.paved.pav_nodes):
            filepath = os.path.join(self.output_dir, "PAVED.3B")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.paved.pav_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "PAVE id '"
                            + dct["id"]
                            + "' ar "
                            + dct["ar"]
                            + " lv "
                            + dct["lv"]
                            + " sd 'sto_"
                            + dct["id"]
                            + "' ss 0 qc 0 "
                            + dct["qc"]
                            + " 0 qo 2 2 ms '"
                            + dct["ms"]
                            + "' is "
                            + dct["is"]
                            + " np "
                            + dct["np"]
                            + " dw 'Def_DWA' ro "
                            + dct["ro"]
                            + " ru "
                            + dct["ru"]
                            + " qh '' pave\n"
                        )

            filepath = os.path.join(self.output_dir, "PAVED.STO")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.paved.pav_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "STDF id 'sto_"
                            + dct["id"]
                            + "' nm 'sto_"
                            + dct["id"]
                            + "' ms "
                            + dct["strs"]
                            + " is 0 mr "
                            + dct["sews"]
                            + " 0.0 ir 0.0 0.0 stdf\n"
                        )

            filepath = os.path.join(self.output_dir, "PAVED.DWA")
            with open(filepath, "w") as f:
                f.write(
                    "DWA id 'Def_DWA' nm 'Def_DWA' do 1 wc 0 wd 0 wh 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 dwa\n"
                )

            filepath = os.path.join(self.output_dir, "WWTP.3B")
            with open(filepath, "w") as f:
                f.write("WWTP id 'WWTP' tb 0  wwtp\n")

            filepath = os.path.join(self.output_dir, "WWTP.tbl")
            with open(filepath, "w") as f:
                f.write("\n")

    def write_greenhouse(self):
        """
        Method to write all files associated with greenhouse nodes: GREENHSE.3B, GREENHSE.RF and GREENHSE.SIL. The latter contains only one definition.

        Returns
        -------
        None.

        """
        if any(self.rrmodel.greenhouse.gh_nodes):
            filepath = os.path.join(self.output_dir, "GREENHSE.3B")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.greenhouse.gh_nodes.items():                    
                    if float(dct["ar"]) > 0.0:
                        area_string = ' '.join(['0' for _ in range(int(dct['basin_storage_class'])-1)]+ [dct['ar']]+['0' for _ in range(10-(int(dct['basin_storage_class'])))])
                        f.write(
                            "GRHS id '"
                            + dct["id"]
                            + "' na 10 ar "
                            + area_string
                            + " sl "
                            + dct["sl"]
                            + " as 0 si 'Def_silo' sd 'sto_"
                            + dct["id"]
                            + "' ms '"
                            + dct["ms"]
                            + "' is "
                            + dct["is"]
                            + " grhs\n"
                        )

            filepath = os.path.join(self.output_dir, "GREENHSE.SIL")
            with open(filepath, "w") as f:
                f.write("SILO id 'Def_silo' nm 'Def_silo' sc 0.0 pc 0.0 silo\n")

            filepath = os.path.join(self.output_dir, "GREENHSE.RF")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "STDF id 'sto_"
                            + dct["id"]
                            + "' nm 'sto_"
                            + dct["id"]
                            + "' mk "
                            + dct["sd"]
                            + " ik 0 stdf\n"
                        )

    def write_openwater(self):
        """
         Method to write OPENWATE.3B file.

        Returns
        -------
        None.

        """
        # write openwater.3b
        if any(self.rrmodel.openwater.ow_nodes):
            filepath = os.path.join(self.output_dir, "OPENWATE.3B")
            with open(filepath, "w") as f:
                for _, dct in self.rrmodel.openwater.ow_nodes.items():
                    if float(dct["ar"]) > 0.0:
                        f.write(
                            "OWRR id '"
                            + dct["id"]
                            + "' ar "
                            + dct["ar"]
                            + " ms '"
                            + dct["ms"]
                            + "' owrr\n"
                        )

    def write_meteo(self):
        """
        Method to write meteofiles (DEFAULT.BUI and DEFAULT.EVP) based on values per catchment.

        Returns
        -------
        None.

        """
        # precip
        filepath = os.path.join(self.output_dir, "DEFAULT.BUI")
        if isinstance(self.rrmodel.external_forcings.precip, str):
            shutil.copy(self.rrmodel.external_forcings.precip, filepath)
        else:
            df = {}
            for ii, ms in enumerate(self.rrmodel.external_forcings.precip.items()):
                df[ms[0]] = pd.Series(ms[1]["precip"])

            f = open(filepath, "w")
            f.write("*Name of this file: c:\\Result\\1058\\DEFAULT.BUI\n")
            f.close()
            f = open(filepath, "a")
            f.write("*Date and time of construction: 00/00/2000 00:00:00.\n")
            f.write("1\n")
            f.write("*Aantal stations\n")
            f.write(str(len(self.rrmodel.external_forcings.precip)) + "\n")
            f.write("*Namen van stations\n")
            for ms in self.rrmodel.external_forcings.precip.items():
                f.write("'" + ms[0] + "'" + "\n")
            f.write(
                "*Aantal gebeurtenissen (omdat het 1 bui betreft is dit altijd 1)\n"
            )
            f.write("*en het aantal seconden per waarnemingstijdstap\n")
            f.write(
                f'1 {(ms[1]["precip"].index[1]-ms[1]["precip"].index[0]).total_seconds():.0f}\n'
            )
            f.write("*Elke commentaarregel wordt begonnen met een * (asterisk).\n")
            f.write(
                "*Eerste record bevat startdatum en -tijd, lengte van de gebeurtenis in dd hh mm ss\n"
            )
            f.write("*Het format is: yyyymmdd:hhmmss:ddhhmmss")
            f.write("*Daarna voor elk station de neerslag in mm per tijdstap.\n")
            duration = (
                ms[1]["precip"].index[-1] - ms[1]["precip"].index[0]
            ).total_seconds()
            f.write(
                ms[1]["precip"].index[0].strftime("%Y %#m %#d %#H %#M %#S ")
                + str(int(duration / 86400.0))
                + " "
                + str(int(np.remainder(duration, 86400.0) / 3600.0))
                + " 0 0\n"
            )
            f.close()
            self._dict_to_df()
            self.precip_df.to_csv(
                filepath,
                header=False,
                index=False,
                sep=" ",
                float_format="%.3f",
                mode="a",
            )
        #%%

        #%%%
        # evaporation
        filepath = os.path.join(self.output_dir, "DEFAULT.EVP")
        if isinstance(self.rrmodel.external_forcings.evap, str):
            shutil.copy(self.rrmodel.external_forcings.evap, filepath)
        else:
            f = open(filepath, "w")
            f.write("*Verdampingsfile\n")
            f.write("*Meteo data: evaporation intensity in mm/day\n")
            f.write("*First record: start date, data in mm/day\n")
            f.write(
                "*Datum (year month day), verdamping (mm/dag) voor elk weerstation\n"
            )
            f.write("*jaar maand dag verdamping[mm]\n")
            f.close()
            table = list(self.rrmodel.external_forcings.evap.values())[0]["evap"]
            table = table.sort_index()
            date_format = "%#Y  %#m  %#d "
            float_format = "%.3f"
            sep = " "
            with open(filepath, "a") as f:
                for ts, val in table.items():
                    date = ts.strftime(date_format)
                    f.write(f"{date}{sep}{float_format % val}\n")

    def _dict_to_df(self):
        """
        Method to convert the precip dictionary to a dataframe
        Returns
        -------
        None.

        """
        df = {}
        for ii, stat in enumerate(self.rrmodel.external_forcings.precip.items()):
            df[stat[0]] = pd.Series(stat[1]["precip"])
        self.precip_df = pd.DataFrame(df)

    def _write_dict(self, f, dct, header, endline):
        """
        Write a dictionary to file, without transforming with mapping
        """
        # Write header
        f.write(f"[{header}]\n")
        for key, value in dct.items():
            f.write(f"    {key:20s} = {value:20s}\n")
        f.write(endline)

    def change_d3b_parameter(self, parameter, value):
        """
        Change parameter value in delft3b.ini file
        """
        with open(os.path.join(self.output_dir, "DELFT_3B.INI")) as f:
            lines = f.readlines()

        for i, line in enumerate(lines):
            if line.lower().strip().startswith(parameter.lower()):
                items = re.split("=", line)
                key, oldvalue = items[0], items[1]
                if not key.strip().lower() == parameter.lower():
                    continue
                lines[i] = line.replace(
                    "=" + oldvalue, f"={value}\n"
                )
                break

        with open(os.path.join(self.output_dir, "DELFT_3B.INI"), "w") as f:
            f.write("".join(lines))

change_d3b_parameter(parameter, value)

Change parameter value in delft3b.ini file

Source code in hydrolib/dhydamo/io/drrwriter.py

def change_d3b_parameter(self, parameter, value):
    """
    Change parameter value in delft3b.ini file
    """
    with open(os.path.join(self.output_dir, "DELFT_3B.INI")) as f:
        lines = f.readlines()

    for i, line in enumerate(lines):
        if line.lower().strip().startswith(parameter.lower()):
            items = re.split("=", line)
            key, oldvalue = items[0], items[1]
            if not key.strip().lower() == parameter.lower():
                continue
            lines[i] = line.replace(
                "=" + oldvalue, f"={value}\n"
            )
            break

    with open(os.path.join(self.output_dir, "DELFT_3B.INI"), "w") as f:
        f.write("".join(lines))

copyRRFiles()

Method to copy the RR-files from resources.

Returns

bool Succesfull or not.

Source code in hydrolib/dhydamo/io/drrwriter.py

def copyRRFiles(self):
    """
    Method to copy the RR-files from resources.

    Returns
    -------
    bool
        Succesfull or not.

    """
    srcRR = os.path.join(os.path.dirname(__file__), "..", "resources", "RR")
    targetRR = os.path.join(self.output_dir)
    shutil.copytree(srcRR, targetRR)
    return True

write_all()

Wrapper method to write all components

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_all(self):  # write all RR files
    """
    Wrapper method to write all components
    """
    # copy all 'standard' files from the resources folder
    self.copyRRFiles()

    self.write_topology()

    self.write_unpaved()

    self.write_paved()

    self.write_greenhouse()

    self.write_openwater()

    self.write_meteo()

    # Change d3b parameters
    for parameter, value in self.d3b_parameters.items():
        self.change_d3b_parameter(parameter, value)

write_greenhouse()

Method to write all files associated with greenhouse nodes: GREENHSE.3B, GREENHSE.RF and GREENHSE.SIL. The latter contains only one definition.

Returns

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_greenhouse(self):
    """
    Method to write all files associated with greenhouse nodes: GREENHSE.3B, GREENHSE.RF and GREENHSE.SIL. The latter contains only one definition.

    Returns
    -------
    None.

    """
    if any(self.rrmodel.greenhouse.gh_nodes):
        filepath = os.path.join(self.output_dir, "GREENHSE.3B")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.greenhouse.gh_nodes.items():                    
                if float(dct["ar"]) > 0.0:
                    area_string = ' '.join(['0' for _ in range(int(dct['basin_storage_class'])-1)]+ [dct['ar']]+['0' for _ in range(10-(int(dct['basin_storage_class'])))])
                    f.write(
                        "GRHS id '"
                        + dct["id"]
                        + "' na 10 ar "
                        + area_string
                        + " sl "
                        + dct["sl"]
                        + " as 0 si 'Def_silo' sd 'sto_"
                        + dct["id"]
                        + "' ms '"
                        + dct["ms"]
                        + "' is "
                        + dct["is"]
                        + " grhs\n"
                    )

        filepath = os.path.join(self.output_dir, "GREENHSE.SIL")
        with open(filepath, "w") as f:
            f.write("SILO id 'Def_silo' nm 'Def_silo' sc 0.0 pc 0.0 silo\n")

        filepath = os.path.join(self.output_dir, "GREENHSE.RF")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "STDF id 'sto_"
                        + dct["id"]
                        + "' nm 'sto_"
                        + dct["id"]
                        + "' mk "
                        + dct["sd"]
                        + " ik 0 stdf\n"
                    )

write_meteo()

Method to write meteofiles (DEFAULT.BUI and DEFAULT.EVP) based on values per catchment.

Returns

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_meteo(self):
    """
    Method to write meteofiles (DEFAULT.BUI and DEFAULT.EVP) based on values per catchment.

    Returns
    -------
    None.

    """
    # precip
    filepath = os.path.join(self.output_dir, "DEFAULT.BUI")
    if isinstance(self.rrmodel.external_forcings.precip, str):
        shutil.copy(self.rrmodel.external_forcings.precip, filepath)
    else:
        df = {}
        for ii, ms in enumerate(self.rrmodel.external_forcings.precip.items()):
            df[ms[0]] = pd.Series(ms[1]["precip"])

        f = open(filepath, "w")
        f.write("*Name of this file: c:\\Result\\1058\\DEFAULT.BUI\n")
        f.close()
        f = open(filepath, "a")
        f.write("*Date and time of construction: 00/00/2000 00:00:00.\n")
        f.write("1\n")
        f.write("*Aantal stations\n")
        f.write(str(len(self.rrmodel.external_forcings.precip)) + "\n")
        f.write("*Namen van stations\n")
        for ms in self.rrmodel.external_forcings.precip.items():
            f.write("'" + ms[0] + "'" + "\n")
        f.write(
            "*Aantal gebeurtenissen (omdat het 1 bui betreft is dit altijd 1)\n"
        )
        f.write("*en het aantal seconden per waarnemingstijdstap\n")
        f.write(
            f'1 {(ms[1]["precip"].index[1]-ms[1]["precip"].index[0]).total_seconds():.0f}\n'
        )
        f.write("*Elke commentaarregel wordt begonnen met een * (asterisk).\n")
        f.write(
            "*Eerste record bevat startdatum en -tijd, lengte van de gebeurtenis in dd hh mm ss\n"
        )
        f.write("*Het format is: yyyymmdd:hhmmss:ddhhmmss")
        f.write("*Daarna voor elk station de neerslag in mm per tijdstap.\n")
        duration = (
            ms[1]["precip"].index[-1] - ms[1]["precip"].index[0]
        ).total_seconds()
        f.write(
            ms[1]["precip"].index[0].strftime("%Y %#m %#d %#H %#M %#S ")
            + str(int(duration / 86400.0))
            + " "
            + str(int(np.remainder(duration, 86400.0) / 3600.0))
            + " 0 0\n"
        )
        f.close()
        self._dict_to_df()
        self.precip_df.to_csv(
            filepath,
            header=False,
            index=False,
            sep=" ",
            float_format="%.3f",
            mode="a",
        )
    #%%

    #%%%
    # evaporation
    filepath = os.path.join(self.output_dir, "DEFAULT.EVP")
    if isinstance(self.rrmodel.external_forcings.evap, str):
        shutil.copy(self.rrmodel.external_forcings.evap, filepath)
    else:
        f = open(filepath, "w")
        f.write("*Verdampingsfile\n")
        f.write("*Meteo data: evaporation intensity in mm/day\n")
        f.write("*First record: start date, data in mm/day\n")
        f.write(
            "*Datum (year month day), verdamping (mm/dag) voor elk weerstation\n"
        )
        f.write("*jaar maand dag verdamping[mm]\n")
        f.close()
        table = list(self.rrmodel.external_forcings.evap.values())[0]["evap"]
        table = table.sort_index()
        date_format = "%#Y  %#m  %#d "
        float_format = "%.3f"
        sep = " "
        with open(filepath, "a") as f:
            for ts, val in table.items():
                date = ts.strftime(date_format)
                f.write(f"{date}{sep}{float_format % val}\n")

write_openwater()

Method to write OPENWATE.3B file.

Returns

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_openwater(self):
    """
     Method to write OPENWATE.3B file.

    Returns
    -------
    None.

    """
    # write openwater.3b
    if any(self.rrmodel.openwater.ow_nodes):
        filepath = os.path.join(self.output_dir, "OPENWATE.3B")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.openwater.ow_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "OWRR id '"
                        + dct["id"]
                        + "' ar "
                        + dct["ar"]
                        + " ms '"
                        + dct["ms"]
                        + "' owrr\n"
                    )

write_paved()

Method to write all files associated with paved nodes: PAVED.3B, PAVED.STO and PAVED.DWA. The latter contains only one definition.

Returns

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_paved(self):
    """
    Method to write all files associated with paved nodes: PAVED.3B, PAVED.STO and PAVED.DWA. The latter contains only one definition.

    Returns
    -------
    None.

    """
    if any(self.rrmodel.paved.pav_nodes):
        filepath = os.path.join(self.output_dir, "PAVED.3B")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.paved.pav_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "PAVE id '"
                        + dct["id"]
                        + "' ar "
                        + dct["ar"]
                        + " lv "
                        + dct["lv"]
                        + " sd 'sto_"
                        + dct["id"]
                        + "' ss 0 qc 0 "
                        + dct["qc"]
                        + " 0 qo 2 2 ms '"
                        + dct["ms"]
                        + "' is "
                        + dct["is"]
                        + " np "
                        + dct["np"]
                        + " dw 'Def_DWA' ro "
                        + dct["ro"]
                        + " ru "
                        + dct["ru"]
                        + " qh '' pave\n"
                    )

        filepath = os.path.join(self.output_dir, "PAVED.STO")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.paved.pav_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "STDF id 'sto_"
                        + dct["id"]
                        + "' nm 'sto_"
                        + dct["id"]
                        + "' ms "
                        + dct["strs"]
                        + " is 0 mr "
                        + dct["sews"]
                        + " 0.0 ir 0.0 0.0 stdf\n"
                    )

        filepath = os.path.join(self.output_dir, "PAVED.DWA")
        with open(filepath, "w") as f:
            f.write(
                "DWA id 'Def_DWA' nm 'Def_DWA' do 1 wc 0 wd 0 wh 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 dwa\n"
            )

        filepath = os.path.join(self.output_dir, "WWTP.3B")
        with open(filepath, "w") as f:
            f.write("WWTP id 'WWTP' tb 0  wwtp\n")

        filepath = os.path.join(self.output_dir, "WWTP.tbl")
        with open(filepath, "w") as f:
            f.write("\n")

write_topology()

Wrapper to write the topolgy files for RR. The following files are written: 3B_NOD.TP: IDs and coordinates of all nodes: unpaved, paved, greenhouse, open water and boundary. Only nodes with a total area of >0 are written. Nodes are placed up to 20 m to the east and west of the catchment centroid. 3B_LINK.TP: Links between all nodes and the correct boundaries. BOUND3B.3B: coupling information for every RR boundary. BoundaryConditions.bc: see bound3b.3b Returns

---

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_topology(self):
    """
    Wrapper to write the topolgy files for RR. The following files are written:
        3B_NOD.TP: IDs and coordinates of all nodes: unpaved, paved, greenhouse, open water and boundary. Only nodes with a total area of >0 are written. Nodes are placed up to 20 m to the east and west of the catchment centroid.
        3B_LINK.TP: Links between all nodes and the correct boundaries.
        BOUND3B.3B: coupling information for every RR boundary.
        BoundaryConditions.bc: see bound3b.3b
    Returns
    -------
    None.

    """
    filepath = os.path.join(self.output_dir, "3B_NOD.TP")
    with open(filepath, "w") as f:

        # Unpaved nodes
        if any(self.rrmodel.unpaved.unp_nodes):
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    f.write(
                        "NODE id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' ri '-1' mt 1 '2' nt 44 ObID '3B_UNPAVED' px "
                        + dct["px"]
                        + " py "
                        + dct["py"]
                        + " node\n"
                    )

        # Paved nodes
        if any(self.rrmodel.paved.pav_nodes):
            for _, dct in self.rrmodel.paved.pav_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "NODE id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' ri '-1' mt 1 '1' nt 43 ObID '3B_PAVED' px "
                        + dct["px"]
                        + " py "
                        + dct["py"]
                        + " node\n"
                    )

            f.write(
                "NODE id 'WWTP' nm 'WWTP' ri '-1' mt 1 '14' nt 56 ObID '3B_WWTP' px "
                + str(self.wwtp.coords[0][0])
                + " py "
                + str(self.wwtp.coords[0][1])
                + " node\n"
            )
            f.write(
                "NODE id 'WWTP_BND' nm 'WWTP_BND' ri '-1' mt 1 '6' nt 47 ObID '3B_BOUNDARY' px "
                + str(self.wwtp.coords[0][0] + 50.0)
                + " py "
                + str(self.wwtp.coords[0][1] + 50.0)
                + " node\n"
            )

        # Greenhouse nodes
        if any(self.rrmodel.greenhouse.gh_nodes):
            for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "NODE id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' ri '-1' mt 1 '3' nt 45 ObID '3B_GREENHOUSE' px "
                        + dct["px"]
                        + " py "
                        + dct["py"]
                        + " node\n"
                    )

        # Openwater nodes
        if any(self.rrmodel.openwater.ow_nodes):
            for _, dct in self.rrmodel.openwater.ow_nodes.items():
                if float(dct["ar"]) > 0.0:
                    f.write(
                        "NODE id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' ri '-1' mt 1 '21' nt 46 ObID 'OW_PRECIP' px "
                        + dct["px"]
                        + " py "
                        + dct["py"]
                        + " node\n"
                    )

        # Boundary_nodes
        if any(self.rrmodel.external_forcings.boundary_nodes):
            for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
                f.write(
                    "NODE id '"
                    + dct["id"]
                    + "' nm '"
                    + dct["id"]
                    + "' ri '-1' mt 1 '6' nt 78 ObID 'SBK_SBK-3B-NODE' px "
                    + dct["px"]
                    + " py "
                    + dct["py"]
                    + " node\n"
                )

    filepath = os.path.join(self.output_dir, "3B_LINK.TP")
    with open(filepath, "w") as f:

        cnt = 0
        if any(self.rrmodel.unpaved.unp_nodes):
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    cnt += 1
                    f.write(
                        "BRCH id 'link_"
                        + str(cnt)
                        + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                        + dct["id"]
                        + "' en '"
                        + dct["boundary_node"]
                        + "' brch\n"
                    )
        if any(self.rrmodel.paved.pav_nodes):
            for _, dct in self.rrmodel.paved.pav_nodes.items():
                if float(dct["ar"]) > 0.0:
                    cnt += 1
                    f.write(
                        "BRCH id 'link_"
                        + str(cnt)
                        + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                        + dct["id"]
                        + "' en '"
                        + dct["boundary_node"]
                        + "' brch\n"
                    )
                    cnt += 1
                    f.write(
                        "BRCH id 'link_"
                        + str(cnt)
                        + "' ri '-1' mt 1 '1' bt 18 ObID '3B_LINK_RWZI' bn '"
                        + dct["id"]
                        + "' en 'WWTP' brch\n"
                    )
            cnt += 1
            f.write(
                "BRCH id 'link_"
                + str(cnt)
                + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn 'WWTP' en 'WWTP_BND' brch\n"
            )

        if any(self.rrmodel.greenhouse.gh_nodes):
            for _, dct in self.rrmodel.greenhouse.gh_nodes.items():
                if float(dct["ar"]) > 0.0:
                    cnt += 1
                    f.write(
                        "BRCH id 'link_"
                        + str(cnt)
                        + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                        + dct["id"]
                        + "' en '"
                        + dct["boundary_node"]
                        + "' brch\n"
                    )
        if any(self.rrmodel.openwater.ow_nodes):
            for _, dct in self.rrmodel.openwater.ow_nodes.items():
                if float(dct["ar"]) > 0.0:
                    cnt += 1
                    f.write(
                        "BRCH id 'link_"
                        + str(cnt)
                        + "' ri '-1' mt 1 '0' bt 17 ObID '3B_LINK' bn '"
                        + dct["id"]
                        + "' en '"
                        + dct["boundary_node"]
                        + "' brch\n"
                    )

    # bound3b.3b
    filepath = os.path.join(self.output_dir, "Bound3B.3B")
    with open(filepath, "w") as f:
        for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
            f.write("BOUN id '" + dct["id"] + "' bl 2 '0' is 0 boun\n")
        if any(self.rrmodel.paved.pav_nodes):
            f.write("BOUN id 'WWTP_BND' bl 0 -3 is 0 boun\n")

    # BoundaryConditions.bc
    filepath = os.path.join(self.output_dir, "BoundaryConditions.bc")
    header = {"fileVersion": "1.01", "fileType": "boundConds"}
    with open(filepath, "w") as f:
        self._write_dict(f, header, "General", "\n")
        for _, dct in self.rrmodel.external_forcings.boundary_nodes.items():
            temp = {
                "name": "" + dct["id"],
                "function": "constant",
                "quantity": "water_level",
                "unit": "m",
            }
            self._write_dict(f, temp, "Boundary", "    0\n\n")
        if any(self.rrmodel.paved.pav_nodes):
            temp = {
                "name": "WWTP_BND",
                "function": "constant",
                "quantity": "water_level",
                "unit": "m",
            }
            self._write_dict(f, temp, "Boundary", "    0\n\n")

write_unpaved()

Method to write all files associated with unpaved nodes: UNPAVED.3B, UNPAVED.ALF, UNPAVED.STO, UNPAVED.INF and UNPAVED.SEP. All files contain a definition for every node because they may or may not be spatially distributed.

Returns

None.

Source code in hydrolib/dhydamo/io/drrwriter.py

def write_unpaved(self):
    """
    Method to write all files associated with unpaved nodes: UNPAVED.3B, UNPAVED.ALF, UNPAVED.STO, UNPAVED.INF and UNPAVED.SEP. All files contain a definition for every node  because they may or may not be spatially distributed.

    Returns
    -------
    None.

    """
    if any(self.rrmodel.unpaved.unp_nodes):
        filepath = os.path.join(self.output_dir, "UNPAVED.3B")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    f.write(
                        "UNPV id '"
                        + dct["id"]
                        + "' na "
                        + dct["na"]
                        + " ga "
                        + dct["ga"]
                        + " ar "
                        + dct["ar"]
                        + " lv "
                        + dct["lv"]
                        + " co "
                        + dct["co"]
                        + " su "
                        + dct["su"]
                        + " sd 'sto_"
                        + dct["id"]
                        + "' ic 'inf_"
                        + dct["id"]
                        + "' bt "
                        + dct["bt"]
                        + " ed '"
                        + dct["ed"]
                        + "' sp '"
                        + dct["sp"]
                        + "' ig 0 "
                        + dct["ig"]
                        + " mg "
                        + dct["mg"]
                        + " gl "
                        + dct["gl"]
                        + " is "
                        + dct["is"]
                        + " ms '"
                        + dct["ms"]
                        + "' unpv\n"
                    )

        filepath = os.path.join(self.output_dir, "UNPAVED.STO")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    f.write(
                        "STDF id 'sto_"
                        + dct["id"]
                        + "' nm 'sto_"
                        + dct["id"]
                        + "' ml "
                        + dct["sd"]
                        + " il 0 stdf\n"
                    )

        filepath = os.path.join(self.output_dir, "UNPAVED.INF")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    f.write(
                        "INFC id 'inf_"
                        + dct["id"]
                        + "' nm 'inf_"
                        + dct["id"]
                        + "' ic "
                        + dct["ic"]
                        + " infc\n"
                    )

        filepath = os.path.join(self.output_dir, "UNPAVED.ALF")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.unpaved.ernst_defs.items():
                if (
                    np.sum(
                        [
                            float(d)
                            for d in self.rrmodel.unpaved.unp_nodes[_]["ar"].split(
                                " "
                            )
                        ]
                    )
                    > 0.0
                ):
                    f.write(
                        "ERNS id '"
                        + dct["id"]
                        + "' nm '"
                        + dct["id"]
                        + "' cvi "
                        + dct["cvi"]
                        + " cvo 0 "
                        + dct["cvo"]
                        + " cvs "
                        + dct["cvs"]
                        + " lv "
                        + dct["lv"]
                        + " erns\n"
                    )

        filepath = os.path.join(self.output_dir, "UNPAVED.SEP")
        with open(filepath, "w") as f:
            for _, dct in self.rrmodel.unpaved.unp_nodes.items():
                if np.sum([float(d) for d in dct["ar"].split(" ")]) > 0.0:
                    f.write(
                        "SEEP id '"
                        + dct["sp"]
                        + "' nm '"
                        + dct["sp"]
                        + "' co 4 PDIN 0 0 pdin ss 0\n"
                    )
                    f.write("TBLE\n")
                    for row in self.rrmodel.external_forcings.seepage[dct["sp"]][
                        "seepage"
                    ].items():
                        f.write(
                            "'"
                            + row[0].strftime("%Y/%m/%d;%H:%M:%S")
                            + "' "
                            + f"{row[1]:.5f} <\n"
                        )
                    f.write("tble\nseep\n")

idfreader

Functions for reading and writing iMOD Data Files (IDFs) to xarray objects.

The primary functions to use are :func:imod.idf.open and :func:imod.idf.save, though lower level functions are also available.

The function in this file are taken from the imod-python package version 0.10.1.

The only functionality needed in dhydamo is the capability to read IDF files.

affine_from_idf(dataset)

Build a rasterio affine transform from an opened IDF dataset.

Parameters

dataset : xarray.DataArray Single-time, single-layer array as returned by idfreader.open(...).squeeze(), with scalar dx/dy coordinates and pixel-centre x/y coordinates.

Returns

affine.Affine Geotransform for the raster, derived from its pixel-centre coordinates and cell size.

Source code in hydrolib/dhydamo/io/idfreader.py

def affine_from_idf(dataset):
    """Build a rasterio affine transform from an opened IDF dataset.

    Parameters
    ----------
    dataset : xarray.DataArray
        Single-time, single-layer array as returned by
        ``idfreader.open(...).squeeze()``, with scalar ``dx``/``dy``
        coordinates and pixel-centre ``x``/``y`` coordinates.

    Returns
    -------
    affine.Affine
        Geotransform for the raster, derived from its pixel-centre
        coordinates and cell size.
    """
    dx = float(dataset["dx"])
    dy = float(dataset["dy"])
    west = float(dataset["x"].min()) - dx / 2.0
    north = float(dataset["y"].max()) - dy / 2.0
    return from_origin(west, north, dx, -dy)

decompose(path, pattern=None)

Parse a path, returning a dict of the parts, following the iMOD conventions.

Parameters

path : str or pathlib.Path Path to the file. Upper case is ignored. pattern : str, regex pattern, optional If the path is not made up of standard paths, and the default decompose does not produce the right result, specify the used pattern here. See the examples below.

Returns

d : dict Dictionary with name of variable and dimensions

Examples

Decompose a path, relying on default conventions:

decompose("head_20010101_l1.idf")

Do the same, by specifying a format string pattern, excluding extension:

decompose("head_20010101_l1.idf", pattern="{name}_{time}_l{layer}")

This supports an arbitrary number of variables:

decompose("head_slr_20010101_l1.idf", pattern="{name}{scenario}{time}_l{layer}")

The format string pattern will only work on tidy paths, where variables are separated by underscores. You can also pass a compiled regex pattern. Make sure to include the re.IGNORECASE flag since all paths are lowered.

import re pattern = re.compile(r"(?P[\w]+)L(?P[\d+]*)") decompose("headL11", pattern=pattern)

However, this requires constructing regular expressions, which is generally a fiddly process. The website https://regex101.com is a nice help. Alternatively, the most pragmatic solution may be to just rename your files.

Source code in hydrolib/dhydamo/io/idfreader.py

def decompose(path, pattern=None):
    r"""
    Parse a path, returning a dict of the parts, following the iMOD conventions.

    Parameters
    ----------
    path : str or pathlib.Path
        Path to the file. Upper case is ignored.
    pattern : str, regex pattern, optional
        If the path is not made up of standard paths, and the default decompose
        does not produce the right result, specify the used pattern here. See
        the examples below.

    Returns
    -------
    d : dict
        Dictionary with name of variable and dimensions

    Examples
    --------
    Decompose a path, relying on default conventions:

    >>> decompose("head_20010101_l1.idf")

    Do the same, by specifying a format string pattern, excluding extension:

    >>> decompose("head_20010101_l1.idf", pattern="{name}_{time}_l{layer}")

    This supports an arbitrary number of variables:

    >>> decompose("head_slr_20010101_l1.idf", pattern="{name}_{scenario}_{time}_l{layer}")

    The format string pattern will only work on tidy paths, where variables are
    separated by underscores. You can also pass a compiled regex pattern.
    Make sure to include the ``re.IGNORECASE`` flag since all paths are lowered.

    >>> import re
    >>> pattern = re.compile(r"(?P<name>[\w]+)L(?P<layer>[\d+]*)")
    >>> decompose("headL11", pattern=pattern)

    However, this requires constructing regular expressions, which is generally
    a fiddly process. The website https://regex101.com is a nice help.
    Alternatively, the most pragmatic solution may be to just rename your files.
    """
    path = pathlib.Path(path)
    # We'll ignore upper case
    stem = path.stem.lower()

    d = _groupdict(stem, pattern)
    dims = list(d.keys())
    # If name is not provided, generate one from other fields
    if "name" not in d.keys():
        d["name"] = "_".join(d.values())
    else:
        dims.remove("name")

    # String -> type conversion
    if "layer" in d.keys():
        d["layer"] = int(d["layer"])
    if "species" in d.keys():
        d["species"] = int(d["species"])
    if "time" in d.keys():
        d["time"] = to_datetime(d["time"])
    if "steady-state" in d["name"]:
        # steady-state as time identifier isn't picked up by <time>[0-9] regex
        d["name"] = d["name"].replace("_steady-state", "")
        d["time"] = "steady-state"
        dims.append("time")

    d["extension"] = path.suffix
    d["directory"] = path.parent
    d["dims"] = dims
    return d

header(path, pattern)

Read the IDF header information into a dictionary

Source code in hydrolib/dhydamo/io/idfreader.py

def header(path, pattern):
    """Read the IDF header information into a dictionary"""
    attrs = decompose(path, pattern)
    with f_open(path, "rb") as f:
        reclen_id = struct.unpack("i", f.read(4))[0]  # Lahey RecordLength Ident.
        if reclen_id == 1271:
            floatsize = intsize = 4
            floatformat = "f"
            intformat = "i"
            dtype = "float32"
            doubleprecision = False
        # 2296 was a typo in the iMOD manual. Keep 2296 around in case some IDFs
        # were written with this identifier to avoid possible incompatibility
        # issues.
        elif reclen_id == 2295 or reclen_id == 2296:
            floatsize = intsize = 8
            floatformat = "d"
            intformat = "q"
            dtype = "float64"
            doubleprecision = True
        else:
            raise ValueError(
                f"Not a supported IDF file: {path}\n"
                "Record length identifier should be 1271 or 2295, "
                f"received {reclen_id} instead."
            )

        # Header is fully doubled in size in case of double precision ...
        # This means integers are also turned into 8 bytes
        # and requires padding with some additional bytes
        if doubleprecision:
            f.read(4)  # not used

        ncol = struct.unpack(intformat, f.read(intsize))[0]
        nrow = struct.unpack(intformat, f.read(intsize))[0]
        attrs["xmin"] = struct.unpack(floatformat, f.read(floatsize))[0]
        attrs["xmax"] = struct.unpack(floatformat, f.read(floatsize))[0]
        attrs["ymin"] = struct.unpack(floatformat, f.read(floatsize))[0]
        attrs["ymax"] = struct.unpack(floatformat, f.read(floatsize))[0]
        # dmin and dmax are recomputed during writing
        f.read(floatsize)  # dmin, minimum data value present
        f.read(floatsize)  # dmax, maximum data value present
        nodata = struct.unpack(floatformat, f.read(floatsize))[0]
        attrs["nodata"] = nodata
        # flip definition here such that True means equidistant
        # equidistant IDFs
        ieq = not struct.unpack("?", f.read(1))[0]
        itb = struct.unpack("?", f.read(1))[0]

        f.read(2)  # not used
        if doubleprecision:
            f.read(4)  # not used

        if ieq:
            # dx and dy are stored positively in the IDF
            # dy is made negative here to be consistent with the nonequidistant case
            attrs["dx"] = struct.unpack(floatformat, f.read(floatsize))[0]
            attrs["dy"] = -struct.unpack(floatformat, f.read(floatsize))[0]

        if itb:
            attrs["top"] = struct.unpack(floatformat, f.read(floatsize))[0]
            attrs["bot"] = struct.unpack(floatformat, f.read(floatsize))[0]

        if not ieq:
            # dx and dy are stored positive in the IDF, but since the difference between
            # successive y coordinates is negative, it is made negative here
            attrs["dx"] = np.fromfile(f, dtype, ncol)
            attrs["dy"] = -np.fromfile(f, dtype, nrow)

        # These are derived, remove after using them downstream
        attrs["headersize"] = f.tell()
        attrs["ncol"] = ncol
        attrs["nrow"] = nrow
        attrs["dtype"] = dtype

    return attrs

open(path, use_cftime=False, pattern=None)

Open one or more IDF files as an xarray.DataArray.

In accordance with xarray's design, open loads the data of IDF files lazily. This means the data of the IDFs are not loaded into memory until the data is needed. This allows for easier handling of large datasets, and more efficient computations.

Parameters

path : str, Path or list This can be a single file, 'head_l1.idf', a glob pattern expansion, 'head_l*.idf', or a list of files, ['head_l1.idf', 'head_l2.idf']. Note that each file needs to be of the same name (part before the first underscore) but have a different layer and/or timestamp, such that they can be combined in a single xarray.DataArray. use_cftime : bool, optional Use cftime.DatetimeProlepticGregorian instead of np.datetime64[ns] for the time axis.

Dates are normally encoded as ``np.datetime64[ns]``; however, if dates
fall before 1678 or after 2261, they are automatically encoded as
``cftime.DatetimeProlepticGregorian`` objects rather than
``np.datetime64[ns]``.

pattern : str, regex pattern, optional If the filenames do match default naming conventions of {name}_{time}_l{layer}, a custom pattern can be defined here either as a string, or as a compiled regular expression pattern. See the examples below.

Returns

xarray.DataArray A float xarray.DataArray of the values in the IDF file(s). All metadata needed for writing the file to IDF or other formats using imod.rasterio are included in the xarray.DataArray.attrs.

Examples

Open an IDF file:

da = imod.idf.open("example.idf")

Open an IDF file, relying on default naming conventions to identify layer:

da = imod.idf.open("example_l1.idf")

Open an IDF file, relying on default naming conventions to identify layer and time:

head = imod.idf.open("head_20010101_l1.idf")

Open multiple IDF files, in this case files for the year 2001 for all layers, again relying on default conventions for naming:

head = imod.idf.open("head_2001_l.idf")

The same, this time explicitly specifying name, time, and layer:

head = imod.idf.open("head_2001_l.idf", pattern="{name}_{time}_l{layer}")

The format string pattern will only work on tidy paths, where variables are separated by underscores. You can also pass a compiled regex pattern. Make sure to include the re.IGNORECASE flag since all paths are lowered.

import re pattern = re.compile(r"(?P[\w]+)L(?P[\d+]*)", re.IGNORECASE) head = imod.idf.open("headL11", pattern=pattern)

However, this requires constructing regular expressions, which is generally a fiddly process. Regex notation is also impossible to remember. The website https://regex101.com is a nice help. Alternatively, the most pragmatic solution may be to just rename your files.

Source code in hydrolib/dhydamo/io/idfreader.py

def open(path, use_cftime=False, pattern=None):
    r"""
    Open one or more IDF files as an xarray.DataArray.

    In accordance with xarray's design, ``open`` loads the data of IDF files
    lazily. This means the data of the IDFs are not loaded into memory until the
    data is needed. This allows for easier handling of large datasets, and
    more efficient computations.

    Parameters
    ----------
    path : str, Path or list
        This can be a single file, 'head_l1.idf', a glob pattern expansion,
        'head_l*.idf', or a list of files, ['head_l1.idf', 'head_l2.idf'].
        Note that each file needs to be of the same name (part before the
        first underscore) but have a different layer and/or timestamp,
        such that they can be combined in a single xarray.DataArray.
    use_cftime : bool, optional
        Use ``cftime.DatetimeProlepticGregorian`` instead of `np.datetime64[ns]`
        for the time axis.

        Dates are normally encoded as ``np.datetime64[ns]``; however, if dates
        fall before 1678 or after 2261, they are automatically encoded as
        ``cftime.DatetimeProlepticGregorian`` objects rather than
        ``np.datetime64[ns]``.
    pattern : str, regex pattern, optional
        If the filenames do match default naming conventions of
        {name}_{time}_l{layer}, a custom pattern can be defined here either
        as a string, or as a compiled regular expression pattern. See the
        examples below.

    Returns
    -------
    xarray.DataArray
        A float xarray.DataArray of the values in the IDF file(s).
        All metadata needed for writing the file to IDF or other formats
        using imod.rasterio are included in the xarray.DataArray.attrs.

    Examples
    --------
    Open an IDF file:

    >>> da = imod.idf.open("example.idf")

    Open an IDF file, relying on default naming conventions to identify
    layer:

    >>> da = imod.idf.open("example_l1.idf")

    Open an IDF file, relying on default naming conventions to identify layer
    and time:

    >>> head = imod.idf.open("head_20010101_l1.idf")

    Open multiple IDF files, in this case files for the year 2001 for all
    layers, again relying on default conventions for naming:

    >>> head = imod.idf.open("head_2001*_l*.idf")

    The same, this time explicitly specifying ``name``, ``time``, and ``layer``:

    >>> head = imod.idf.open("head_2001*_l*.idf", pattern="{name}_{time}_l{layer}")

    The format string pattern will only work on tidy paths, where variables are
    separated by underscores. You can also pass a compiled regex pattern.
    Make sure to include the ``re.IGNORECASE`` flag since all paths are lowered.

    >>> import re
    >>> pattern = re.compile(r"(?P<name>[\w]+)L(?P<layer>[\d+]*)", re.IGNORECASE)
    >>> head = imod.idf.open("headL11", pattern=pattern)

    However, this requires constructing regular expressions, which is
    generally a fiddly process. Regex notation is also impossible to
    remember. The website https://regex101.com is a nice help. Alternatively,
    the most pragmatic solution may be to just rename your files.
    """

    if isinstance(path, list):
        return _load(path, use_cftime, pattern, _read, header)
    elif isinstance(path, pathlib.Path):
        path = str(path)

    paths = [pathlib.Path(p) for p in glob.glob(path)]
    n = len(paths)
    if n == 0:
        raise FileNotFoundError(f"Could not find any files matching {path}")
    return _load(paths, use_cftime, pattern, _read, header)

read(path, pattern=None)

Read a single IDF file to a numpy.ndarray

Parameters

path : str or Path Path to the IDF file to be read pattern : str, regex pattern, optional If the filenames do match default naming conventions of {name}_{time}_l{layer}, a custom pattern can be defined here either as a string, or as a compiled regular expression pattern. Please refer to the examples for imod.idf.open.

Returns

numpy.ndarray A float numpy.ndarray with shape (nrow, ncol) of the values in the IDF file. On opening all nodata values are changed to NaN in the numpy.ndarray. dict A dict with all metadata.

Source code in hydrolib/dhydamo/io/idfreader.py

def read(path, pattern=None):
    """
    Read a single IDF file to a numpy.ndarray

    Parameters
    ----------
    path : str or Path
        Path to the IDF file to be read
    pattern : str, regex pattern, optional
        If the filenames do match default naming conventions of
        {name}_{time}_l{layer}, a custom pattern can be defined here either
        as a string, or as a compiled regular expression pattern. Please refer
        to the examples for ``imod.idf.open``.

    Returns
    -------
    numpy.ndarray
        A float numpy.ndarray with shape (nrow, ncol) of the values
        in the IDF file. On opening all nodata values are changed
        to NaN in the numpy.ndarray.
    dict
        A dict with all metadata.
    """
    warnings.warn(
        "The idf.read() function is deprecated. To get a numpy array of an IDF, "
        "use instead: imod.idf.open(path).values",
        FutureWarning,
    )
    attrs = header(path, pattern)
    headersize = attrs.pop("headersize")
    nrow = attrs.pop("nrow")
    ncol = attrs.pop("ncol")
    nodata = attrs.pop("nodata")
    dtype = attrs.pop("dtype")
    return _read(path, headersize, nrow, ncol, nodata, dtype), attrs