Parameters

Parameters

hapi.rrm.parameters.Parameters

Parameter distribution class for hydrological model calibration.

The Parameters class distributes values from a parameter vector during the calibration process into a 3D array, handling lumped parameters and hydrologic response units (HRUs).

Source code in src/hapi/rrm/parameters.py

class Parameters:
    """Parameter distribution class for hydrological model calibration.

    The Parameters class distributes values from a parameter vector during
    the calibration process into a 3D array, handling lumped parameters
    and hydrologic response units (HRUs).
    """

    def __init__(
        self,
        raster: Dataset,
        no_parameters: int,
        no_lumped_par: int = 0,
        lumped_par_pos: list[int] | None = None,
        lake: bool = False,
        snow: bool = False,
        hru: bool = False,
        function: int = 1,
        k_upper_bound: int = 1,
        k_lower_bound: int = 50,
        muskingum: bool = False,
    ):
        """Initialize the Parameters class.

        To initiate the Parameters class, you have to provide the Flow Acc
        raster.

        Args:
            raster: A pyramids ``Dataset`` to get the spatial information
                of the catchment (DEM, flow accumulation or flow direction
                raster), read it using ``Dataset.read_file``.
            no_parameters: Number of parameters in the HBV model.
            no_lumped_par: Number of lumped parameters. You have to enter
                the value of the lumped parameter at the end of the list.
                Defaults to 0 (no lumped parameters).
            lumped_par_pos: List of the order or position of lumped
                parameters among all the parameters of the lumped model
                (order starts from 0 to the length of the model
                parameters). Defaults to None (empty). The following
                order of parameters is used for the lumped HBV model:
                [ltt, utt, rfcf, sfcf, ttm, cfmax, cwh, cfr, fc, beta,
                e_corr, etf, lp, c_flux, k, k1, alpha, perc, pcorr,
                Kmuskingum, Xmuskingum].
            lake: True if there is a lake, False otherwise.
                Defaults to False.
            snow: True to run the snow-related processes, False otherwise.
                When True, parameters related to snow simulation have to
                be provided. Defaults to False.
            hru: True if the parameters will consider using HRUs.
                Defaults to False.
            function: Which parameter-distribution strategy to bind to
                :attr:`Function`. One of ``1`` (:meth:`par3d_lumped`), ``2``
                (:meth:`par3d`), ``3`` (:meth:`par2d_lumped_k1_lake`) or ``4``
                (:meth:`hydrologic_response_units`). Defaults to 1. Any other value
                raises :class:`ValueError`. Note that ``hru=True`` overrides the choice
                with :meth:`hydrologic_response_units` regardless, but the selector is
                still validated so a typo is not masked.
            k_upper_bound: Upper bound of K value (traveling time in
                muskingum routing method). Defaults to 1 hour.
            k_lower_bound: Lower bound of K value (traveling time in
                muskingum routing method). Defaults to 50.
            muskingum: True if the routing function is muskingum.
                Defaults to False.

        Raises:
            TypeError: If `raster` is not a pyramids Dataset.
            ValueError: If `function` is not one of the ints 1, 2, 3 or 4. A `bool`,
                a `float` such as ``2.0``, and an unhashable value are all rejected
                rather than coerced or allowed to raise `TypeError`.
            AssertionError: If `no_parameters` is not an integer, if
                `no_lumped_par` is not an integer, or if the length of
                `lumped_par_pos` does not match `no_lumped_par`.
            ValueError: If `lumped_par_pos` is not a list when
                `no_lumped_par` >= 1.

        Note:
            Cells outside the catchment are identified by pyramids via
            ``read_array(masked=True)`` and stored as ``NaN`` in
            :attr:`raster_array`, which is promoted to floating point so it can hold
            them. ``no_elem``, ``celli``/``cellj`` and the width of :attr:`Par2d` all
            derive from that mask, so the parameter vector length follows the raster's
            real domain.

        Examples:
            - Build the distributor from a small raster and inspect the domain it
              derived. The bottom-right cell is no-data, leaving three cells to
              parameterise:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from hapi.rrm.parameters import Parameters
                >>> raster = Dataset.create_from_array(
                ...     np.array([[1, 2], [3, -9999]], dtype="int32"),
                ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
                ...     no_data_value=-9999,
                ... )
                >>> distributor = Parameters(raster, 12)
                >>> distributor.no_elem
                3
                >>> distributor.Par2d.shape
                (12, 3)
                >>> list(zip(distributor.celli, distributor.cellj))
                [(0, 0), (0, 1), (1, 0)]

                ```
            - A real value within 0.1% of the sentinel is kept, so it is treated as a
              catchment cell and widens the parameter array:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from hapi.rrm.parameters import Parameters
                >>> raster = Dataset.create_from_array(
                ...     np.array([[1, 2], [3, -9990]], dtype="int32"),
                ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
                ...     no_data_value=-9999,
                ... )
                >>> distributor = Parameters(raster, 12)
                >>> distributor.no_elem
                4
                >>> float(distributor.raster_array[1, 1])
                -9990.0

                ```
        """
        if lumped_par_pos is None:
            lumped_par_pos = []

        if not isinstance(raster, Dataset):
            raise TypeError(
                "raster should be a pyramids Dataset, read it using pyramids.dataset.Dataset.read_file"
            )
        assert isinstance(no_parameters, int), " no_parameters should be integer number"
        assert isinstance(no_lumped_par, int), (
            "no of lumped parameters should be integer"
        )

        if no_lumped_par >= 1:
            if isinstance(lumped_par_pos, list):
                assert no_lumped_par == len(lumped_par_pos), (
                    f"you have to entered {no_lumped_par} no of lumped parameters but only {len(lumped_par_pos)} "
                    f"position "
                )
            else:  # if not int or list
                raise ValueError(
                    "you have one or more lumped parameters, so the position has to be entered as a list"
                )

        # Reject an unrecognised selector here rather than leaving `Function` unbound:
        # it is invoked on every calibration iteration, so a silent miss surfaces far
        # from the mistake as a bare AttributeError.
        strategies = {
            1: self.par3d_lumped,
            2: self.par3d,
            3: self.par2d_lumped_k1_lake,
            4: self.hydrologic_response_units,
        }
        # Test the type before the membership lookup: an unhashable selector (a list,
        # a set) raises TypeError from `in` rather than the documented ValueError, and
        # bool is a subclass of int, so True would otherwise silently select strategy 1.
        # A float is rejected outright rather than coerced -- 2.0 is a caller error, not
        # a spelling of 2.
        if isinstance(function, bool) or not isinstance(function, int):
            raise ValueError(
                f"function must be one of {sorted(strategies)}; got {function!r} of type "
                f"{type(function).__name__}. 1 = par3d_lumped, 2 = par3d, "
                "3 = par2d_lumped_k1_lake, 4 = hydrologic_response_units."
            )
        if function not in strategies:
            raise ValueError(
                f"function must be one of {sorted(strategies)}; got {function!r}. "
                "1 = par3d_lumped, 2 = par3d, 3 = par2d_lumped_k1_lake, "
                "4 = hydrologic_response_units."
            )
        self.Lake = lake
        self.Snow = snow
        self.no_lumped_par = no_lumped_par
        self.lumped_par_pos = lumped_par_pos
        self.HRUs = hru
        self.Kub = k_upper_bound
        self.Klb = k_lower_bound
        self.Maskingum = muskingum
        # read the raster
        self.raster = raster
        # No-data masking is delegated to pyramids: vectorised, dtype-aware, and it
        # also honours the band's GDAL mask band. Filling with NaN keeps the
        # float-array-with-NaN contract the rest of this class relies on.
        self.raster_array = np.ma.filled(
            raster.read_array(band=0, masked=True).astype(float), np.nan
        )
        # get the shape of the raster
        self.rows = raster.rows
        self.cols = raster.columns
        # get the no_value of in the raster
        self.noval = raster.no_data_value[0]
        if self.noval is None:
            warnings.warn(
                "the raster declares no no-data value, so every cell is treated as "
                "inside the catchment and the parameter vector is sized for the whole "
                "grid. If it has a sentinel, set it on the band.",
                UserWarning,
                stacklevel=2,
            )

        # count the number of non-empty cells
        if self.HRUs:
            self.values = list(
                set(
                    [
                        int(self.raster_array[i, j])
                        for i in range(self.rows)
                        for j in range(self.cols)
                        if not np.isnan(self.raster_array[i, j])
                    ]
                )
            )
            self.no_elem = len(self.values)
        else:
            # Count the cells the pyramids mask left intact (see
            # Catchment.read_flow_acc for why not Dataset.count_domain_cells).
            self.no_elem = int(np.count_nonzero(~np.isnan(self.raster_array)))

        self.no_parameters = no_parameters

        # store the indexes of the non-empty cells
        self.celli = []
        self.cellj = []
        for i in range(self.rows):
            for j in range(self.cols):
                if not np.isnan(self.raster_array[i, j]):
                    self.celli.append(i)
                    self.cellj.append(j)

        # create an empty 3D array [[raster dimension], no_parameters]
        self.Par3d = np.zeros([self.rows, self.cols, self.no_parameters]) * np.nan

        if no_lumped_par >= 1:
            # parameters in an array
            # remove a place for the lumped parameter (k1) lower zone coefficient
            self.no_parameters = self.no_parameters - no_lumped_par

        # all parameters lumped and distributed
        self.totnumberpar = self.no_parameters * self.no_elem + no_lumped_par
        # parameters in array
        # create a 2d array [no_parameters, no_cells]
        self.Par2d = np.zeros(
            shape=(self.no_parameters, self.no_elem), dtype=np.float32
        )

        self.Function = strategies[function]
        # to overwrite any choice user choose if the is HRUs
        if self.HRUs == 1:
            self.Function = self.hydrologic_response_units

        self.parameters_number()

        pass

    def par3d(self, par_g: list | np.ndarray):  # , kub=1,klb=0.5, Maskingum=True
        """Distribute parameters horizontally across grid cells.

        Takes a list of parameters (saved as one column or generated as a
        1D list from an optimization algorithm) and distributes them
        horizontally on the number of cells given by a raster.

        Args:
            par_g: 1D list or numpy array of parameters. For totally
                distributed parameters, the length should be
                ``no_elem * no_parameters``. For lumped parameters, the
                lumped parameter values should be appended at the end.

        Raises:
            AssertionError: If the length of `par_g` does not match the
                expected number of parameters based on the number of
                elements and lumped parameters.
        """
        # input data validation
        # data type
        # assert type(par_g)==np.ndarray or type(par_g)==list, "par_g should be of type 1d array or list"
        # assert isinstance(kub,numbers.Number) , " kub should be a number"
        # assert isinstance(klb,numbers.Number) , " klb should be a number"

        # input values
        if self.no_lumped_par > 0:
            par_no = (self.no_elem * self.no_parameters) + self.no_lumped_par

            assert len(par_g) == par_no, (
                f"As there is {self.no_lumped_par} lumped parameters, length of input parameters should be "
                f"{self.no_elem}"
                + f"*({self.no_parameters + self.no_lumped_par} - {self.no_lumped_par}) + {self.no_lumped_par} = "
                + f"{self.no_elem * (self.no_parameters - self.no_lumped_par) + self.no_lumped_par} not {len(par_g)}"
                + " probably you have to add the value of the lumped parameter at the end of the list"
            )
        else:
            # if there are no lumped parameters
            par_no = self.no_elem * self.no_parameters
            assert len(par_g) == par_no, (
                f"As there is no lumped parameters length of input parameters should be {self.no_elem} * "
                + f"{self.no_parameters} = {self.no_elem * self.no_parameters}"
            )

        # parameters in array
        # create a 2d array [no_parameters, no_cells]
        self.Par2d = np.ones((self.no_parameters, self.no_elem))  # type: ignore[assignment]
        # take the parameters from the generated parameters or the 1D list and
        # assign them to each cell
        for i in range(self.no_elem):
            self.Par2d[:, i] = par_g[
                i * self.no_parameters : (i * self.no_parameters) + self.no_parameters
            ]

        # lumped parameters
        if self.no_lumped_par > 0:
            for i in range(self.no_lumped_par):
                # create a list with the value of the lumped parameter(k1)
                # (stored at the end of the list of the parameters)
                pk1 = (
                    np.ones((1, self.no_elem))
                    * par_g[(self.no_parameters * np.shape(self.Par2d)[1]) + i]
                )
                # put the list of parameter k1 at the 6th row.
                self.Par2d = np.vstack(
                    [
                        self.Par2d[: self.lumped_par_pos[i], :],
                        pk1,
                        self.Par2d[self.lumped_par_pos[i] :, :],
                    ]
                )

        # assign the parameters from the array (no_parameters, no_cells) to
        # the spatially corrected location in par2d
        for i in range(self.no_elem):
            self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

        # calculate the value of k(travelling time in muskingum based on value of
        # x and the position and upper, lower bound of k value

        # if Maskingum:
        #     for i in range(self.no_elem):
        #         self.Par3d[self.celli[i],self.cellj[i],-2]=
        #         Parameters.calculateK(
        #               self.Par3d[self.celli[i], self.cellj[i],-1], self.Par3d[self.celli[i], self.cellj[i],-2], kub,
        #               klb
        #              )

    def par3d_lumped(
        self, par_g: list | np.ndarray
    ):  # , kub=1, klb=0.5, Maskingum = True
        r"""Distribute lumped parameters horizontally across grid cells.

        Takes a list of parameters (saved as one column or generated as a
        1D list from an optimization algorithm) and distributes them
        horizontally on the number of cells given by a raster, where all
        parameters are lumped (same value for every cell).

        Args:
            par_g: 1D list or numpy array of lumped parameters.
                The length should equal ``no_parameters``.

        Raises:
            ValueError: If `par_g` is not a numpy ndarray or a list.
        """
        # input data validation
        # data type
        if not (isinstance(par_g, np.ndarray) or isinstance(par_g, list)):
            raise ValueError("par_g should be of type 1d array or list")
        # assert isinstance(kub,numbers.Number) , " kub should be a number"
        # assert isinstance(klb,numbers.Number) , " klb should be a number"

        # take the parameters from the generated parameters or the 1D list and
        # assign them to each cell
        for i in range(self.no_elem):
            self.Par2d[:, i] = par_g

        # assign the parameters from the array (no_parameters, no_cells) to
        # the spatially corrected location in par2d
        for i in range(self.no_elem):
            self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

        # calculate the value of k(travelling time in muskingum based on value of
        # x and the position and upper, lower bound of k value
        # if Maskingum == True:
        #     for i in range(self.no_elem):
        #         self.Par3d[self.celli[i],self.cellj[i],-2] = Parameters.calculateK(
        #         self.Par3d[self.celli[i],self.cellj[i],-1], self.Par3d[self.celli[i],self.cellj[i],-2],kub,klb)

    @staticmethod
    def calculate_k(
        x: float, position: int, upper_bound: float, lower_bound: float
    ) -> float:
        """Calculate K parameter for Muskingum routing.

        Takes the value of x parameter and generates 100 random values of
        the K parameter between the upper and lower constraints, then
        returns the value corresponding to the given position.

        Args:
            x: Weighting coefficient to determine the linearity of the
                water surface (one of the parameters of the Muskingum
                routing method).
            position: Random position between upper and lower bounds of
                the K parameter.
            upper_bound: Upper bound for the K parameter.
            lower_bound: Lower bound for the K parameter.

        Returns:
            The K parameter value corresponding to the given position
                within the constrained range.
        """
        # k has to be smaller than this constraint
        constraint1 = 0.5 * 1 / (1 - x)
        # k has to be greater than this constraint
        constraint2 = 0.5 * 1 / x
        # if constraint is higher than UB take UB
        if constraint2 >= upper_bound:
            constraint2 = upper_bound
        # if constraint is lower than LB take UB
        if constraint1 <= lower_bound:
            constraint1 = lower_bound

        generated_k = np.linspace(constraint1, constraint2, 50)
        k = generated_k[int(round(position, 0))]
        return float(k)

    def par2d_lumped_k1_lake(
        self, par_g: list | np.ndarray, no_parameters_lake: int
    ):  # ,kub,klb
        """Distribute parameters with a lumped K1 and lake parameters.

        Takes a list of parameters and distributes them horizontally on
        the number of cells given by a raster. All parameters are
        distributed except the lower zone coefficient (K1), which is
        lumped and appended at the end of the parameter list. Lake
        parameters are extracted from the end of the parameter list.

        Args:
            par_g: 1D list or numpy array of parameters. Each cell's
                distributed parameters are listed sequentially, followed
                by the lumped K1 value, followed by lake parameters at
                the end. For example, with 14 cells and 11 distributed
                parameters: ``14 * 11 = 154 + 1 (K1) = 155``.
            no_parameters_lake: Number of lake parameters to extract
                from the end of `par_g`.
        """
        # parameters in array
        # remove a place for the lumped parameter (k1) lower zone coefficient
        no_parameters = self.no_parameters - 1

        # create a 2d array [no_parameters, no_cells]
        self.Par2d = np.ones((no_parameters, self.no_elem))  # type: ignore[assignment]

        # take the parameters from the generated parameters or the 1D list and
        # assign them to each cell
        for i in range(self.no_elem):
            self.Par2d[:, i] = par_g[
                i * no_parameters : (i * no_parameters) + no_parameters
            ]

        # create a list with the value of the lumped parameter(k1)
        # (stored at the end of the list of the parameters)
        pk1 = (
            np.ones((1, self.no_elem))
            * par_g[(np.shape(self.Par2d)[0] * np.shape(self.Par2d)[1])]
        )

        # put the list of parameter k1 at the 6 row
        self.Par2d = np.vstack([self.Par2d[:6, :], pk1, self.Par2d[6:, :]])

        # assign the parameters from the array (no_parameters, no_cells) to
        # the spatially corrected location in par2d
        for i in range(self.no_elem):
            self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

        # calculate the value of k(travelling time in muskingum based on value of
        # x and the position and upper, lower bound of k value
        # for i in range(self.no_elem):
        #     self.Par3d[self.celli[i],self.cellj[i],-2] = Parameters.calculateK(
        #     self.Par3d[self.celli[i],self.cellj[i],-1],self.Par3d[self.celli[i],self.cellj[i],-2],kub,klb)

        # lake parameters
        self.lake_par = par_g[len(par_g) - no_parameters_lake :]
        # self.lake_par[-2] = Parameters.calculateK(self.lake_par[-1],self.lake_par[-2],kub,klb)

        # return self.Par3d, lake_par

    def hydrologic_response_units(self, par_g: list | np.ndarray):  # ,kub=1,klb=0.5
        """Distribute parameters using Hydrologic Response Units (HRUs).

        Takes a list of parameters (saved as one column or generated as a
        1D list from an optimization algorithm) and distributes them
        horizontally on the number of cells given by a raster. The input
        raster should be a classified raster (by numbers) into classes to
        define the HRUs. Each HRU receives the same set of generated
        parameters.

        Args:
            par_g: 1D list or numpy array of parameters. For HRU without
                lumped parameters, the length should be
                ``no_elem * no_parameters``. For HRU with lumped
                parameters, the lumped parameter values should be
                appended at the end.

        Raises:
            ValueError: If `par_g` is not a numpy ndarray or a list, or
                if the length of `par_g` does not match the expected
                number of parameters.
            AssertionError: If there are lumped parameters and the length
                of `par_g` does not match the expected total.
        """
        # input data validation
        # data type
        if not (isinstance(par_g, np.ndarray) or isinstance(par_g, list)):
            raise ValueError("par_g should be of type 1d array or list")
        # assert isinstance(kub,numbers.Number) , " kub should be a number"
        # assert isinstance(klb,numbers.Number) , " klb should be a number"

        # input values
        if self.no_lumped_par > 0:
            par_no = (self.no_elem * self.no_parameters) + self.no_lumped_par
            assert len(par_g) == par_no, (
                f"As there is {self.no_lumped_par} lumped parameters, length of input parameters should be "
                f"{self.no_elem}*({self.no_parameters}-{self.no_lumped_par})+{self.no_lumped_par}="
                + str(
                    self.no_elem * (self.no_parameters - self.no_lumped_par)
                    + self.no_lumped_par
                )
                + f" not {len(par_g)} probably you have to add the value of the lumped parameter at the end of the list"
            )
        else:
            # if there is no lumped parameters
            if not len(par_g) == self.no_elem * self.no_parameters:
                raise ValueError(
                    f"As there is no lumped parameters length of input parameters should be {self.no_elem}*"
                    f"{self.no_parameters}={self.no_elem * self.no_parameters}"
                )

        # take the parameters from the generated parameters or the 1D list and
        # assign them to each cell
        self.Par2d = np.zeros(
            shape=(self.no_parameters, self.no_elem), dtype=np.float64
        )
        for i in range(self.no_elem):
            self.Par2d[:, i] = par_g[
                i * self.no_parameters : (i * self.no_parameters) + self.no_parameters
            ]

        # lumped parameters
        if self.no_lumped_par > 0:
            for i in range(self.no_lumped_par):
                # create a list with the value of the lumped parameter(k1)
                # (stored at the end of the list of the parameters)
                pk1 = (
                    np.ones((1, self.no_elem))
                    * par_g[(self.no_parameters * np.shape(self.Par2d)[1]) + i]
                )
                # put the list of parameter k1 at the 6 row
                self.Par2d = np.vstack(
                    [
                        self.Par2d[: self.lumped_par_pos[i], :],
                        pk1,
                        self.Par2d[self.lumped_par_pos[i] :, :],
                    ]
                )

        # calculate the value of k(travelling time in muskingum based on value of
        # x and the position and upper, lower bound of k value
        # for i in range(self.no_elem):
        #     self.Par2d[-2,i] = Parameters.calculateK(self.Par2d[-1,i],self.Par2d[-2,i],kub,klb)

        # assign the parameters from the array (no_parameters, no_cells) to
        # the spatially corrected location in par2d each soil type will have the same
        # generated parameters
        for i in range(self.no_elem):
            self.Par3d[self.raster_array == self.values[i]] = self.Par2d[:, i]

    @staticmethod
    def hru_hand(
        dem: Dataset,
        flow_direction: Dataset,
        flow_path_length: Dataset,
        river: Dataset,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Calculate Height Above Nearest Drainage (HAND) for HRU classification.

        Calculates inputs for the HAND method for land use
        classification by tracing flow direction from each cell to the
        nearest river reach, then computing the elevation difference
        and the flow path distance.

        Args:
            dem: A pyramids `Dataset` of the DEM raster.
            flow_direction: A pyramids `Dataset` of the flow direction
                raster.
            flow_path_length: A pyramids `Dataset` of the flow path
                length raster.
            river: A pyramids `Dataset` of the river location raster,
                where cells with value 1 indicate river presence.

        Returns:
            A tuple of two numpy ndarrays:
                - hand: Height above nearest drainage for each cell.
                - dist_to_nearest_drain: Distance to nearest drainage
                  for each cell.

        Raises:
            ValueError: If the catchment boundaries contain anomalies
                (e.g., after cropping with a polygon).
        """
        # Use DEM raster information to run all loops
        dem_a = dem.read_array(band=0)
        no_val = np.float32(dem.no_data_value[0])
        rows = dem.rows
        cols = dem.columns

        # get the indices of the flow direction path
        fd_index = DEM(flow_direction.raster).flow_direction_index()

        # read the river location raster
        river_a = river.read_array(band=0)

        # read the flow path length raster
        fpl_a = flow_path_length.read_array(band=0)

        # trace the flow direction to the nearest river reach and store the location
        # of that nearest reach
        nearest_network = Parameters._trace_nearest_drainage(
            dem_a, no_val, river_a, fd_index, rows, cols
        )

        # the elevation difference to the nearest drainage cell is the height above
        # nearest drainage; the same difference over the flow path length raster is
        # the distance to that drainage
        hand = Parameters._difference_to_nearest_drainage(
            dem_a, dem_a, no_val, nearest_network, rows, cols
        )
        dist_to_nearest_drain = Parameters._difference_to_nearest_drainage(
            fpl_a, dem_a, no_val, nearest_network, rows, cols
        )

        return hand, dist_to_nearest_drain

    @staticmethod
    def _trace_nearest_drainage(
        dem_a: np.ndarray,
        no_val: float | np.floating,
        river_a: np.ndarray,
        fd_index: np.ndarray,
        rows: int,
        cols: int,
    ) -> np.ndarray:
        """Locate the nearest downstream river cell for every domain cell.

        Args:
            dem_a (np.ndarray): DEM values, used to identify domain cells.
            no_val (float): No-data value of the DEM.
            river_a (np.ndarray): River raster; a value of 1 marks a river cell.
            fd_index (np.ndarray): Downstream cell indices, shape (rows, cols, 2).
            rows (int): Number of raster rows.
            cols (int): Number of raster columns.

        Returns:
            np.ndarray: Array of shape (rows, cols, 2) holding the row and column
                of the nearest drainage cell, NaN outside the domain.

        Raises:
            ValueError: If a flow path runs off the grid, which happens when the
                catchment boundary has anomalies.
        """
        nearest_network = np.ones((rows, cols, 2)) * np.nan
        try:
            for i in range(rows):
                for j in range(cols):
                    if dem_a[i, j] == no_val:
                        continue
                    # a river cell is its own nearest drainage
                    row, col = i, j
                    while river_a[row, col] != 1:
                        # not at the river yet, step to the downstream cell
                        row, col = (
                            int(fd_index[row, col, 0]),
                            int(fd_index[row, col, 1]),
                        )
                    nearest_network[i, j, 0] = row
                    nearest_network[i, j, 1] = col
        except (IndexError, ValueError) as e:
            raise ValueError(
                "please check the boundaries of your catchment. After cropping the catchment using a polygon, it "
                "creates anomalies at the boundary"
            ) from e

        return nearest_network

    @staticmethod
    def _difference_to_nearest_drainage(
        values: np.ndarray,
        dem_a: np.ndarray,
        no_val: float | np.floating,
        nearest_network: np.ndarray,
        rows: int,
        cols: int,
    ) -> np.ndarray:
        """Subtract each cell's nearest-drainage value from the cell's own value.

        Args:
            values (np.ndarray): Raster to difference (elevation or flow path length).
            dem_a (np.ndarray): DEM values, used to identify domain cells.
            no_val (float): No-data value of the DEM.
            nearest_network (np.ndarray): Nearest drainage indices from
                `_trace_nearest_drainage`.
            rows (int): Number of raster rows.
            cols (int): Number of raster columns.

        Returns:
            np.ndarray: The difference for every domain cell, NaN elsewhere.
        """
        difference = np.ones((rows, cols)) * np.nan
        for i in range(rows):
            for j in range(cols):
                if dem_a[i, j] == no_val:
                    continue
                drain_row = int(nearest_network[i, j, 0])
                drain_col = int(nearest_network[i, j, 1])
                difference[i, j] = values[i, j] - values[drain_row, drain_col]

        return difference

    def parameters_number(self):
        """Calculate the total number of parameters for the optimization.

        Calculates the number of parameters that the optimization
        algorithm will search for. Use this only in case of totally
        distributed catchment parameters. In case of lumped parameters,
        the number of parameters is the same as the number of parameters
        of the conceptual model.

        The result is stored in the ``ParametersNO`` attribute.

        Note:
            The Parameters object should have the following attributes
            before calling this method: ``raster``, ``no_parameters``,
            ``no_lumped_par``, and ``HRUs``.
        """
        if not self.HRUs:
            if self.no_lumped_par > 0:
                # self.ParametersNO = (self.no_elem *( self.no_parameters - self.no_lumped_par)) + self.no_lumped_par
                self.ParametersNO = (
                    self.no_elem * self.no_parameters
                ) + self.no_lumped_par
            else:
                # if there is no lumped parameters
                self.ParametersNO = self.no_elem * self.no_parameters
        else:
            if self.no_lumped_par > 0:
                # self.ParametersNO = (self.no_elem * (self.no_parameters - self.no_lumped_par)) + self.no_lumped_par
                self.ParametersNO = (
                    self.no_elem * self.no_parameters
                ) + self.no_lumped_par
            else:
                # if there is no lumped parameters
                self.ParametersNO = self.no_elem * self.no_parameters

    def save_parameters(self, path: str | None):
        """Save distributed parameters as raster files.

        Takes the generated 3D parameter array and saves each parameter
        layer as a separate GeoTIFF raster file.

        Args:
            path: Path to the folder where the parameter rasters will
                be saved.

        Raises:
            TypeError: `path` is not a ``str``. Output names are built by concatenation
                (``path + name``), which a ``Path`` does not support.
            FileNotFoundError: The output directory does not exist. Checked up front so
                the failure does not surface midway through writing.

        Note:
            The Parameters object should have the following attributes
            set before calling this method: ``DistParFn``, ``raster``,
            ``Par``, ``no_parameters``, ``snow``, ``kub``, and ``klb``.
        """
        # Not delegated to pyramids: the output names are built by string concatenation
        # below (`path + name`), so this genuinely needs a str, and the directory must
        # exist before the first raster is written. Raised rather than asserted so the
        # checks survive `python -O`.
        if not isinstance(path, str):
            raise TypeError(f"path should be of type string, given: {type(path)}")
        if not os.path.exists(path):
            raise FileNotFoundError(f"{path} you have provided does not exist")

        # save
        if self.Snow == 0:  # now snow subroutine
            pnme = [
                "01_rfcf",
                "02_FC",
                "03_BETA",
                "04_ETF",
                "05_LP",
                "06_K0",
                "07_K1",
                "08_K2",
                "09_UZL",
                "10_PERC",
                "11_Kmuskingum",
                "12_Xmuskingum",
            ]
        else:  # there is snow subtoutine
            pnme = [
                "01_ltt",
                "02_rfcf",
                "03_sfcf",
                "04_cfmax",
                "05_cwh",
                "06_cfr",
                "07_fc",
                "08_beta",
                "09_etf",
                "10_lp",
                "11_k0",
                "12_k1",
                "13_k2",
                "14_uzl",
                "18_perc",
            ]

        if path is not None:
            pnme = [
                path + i + "_" + str(dt.datetime.now())[0:10] + ".tif" for i in pnme
            ]

        for i in range(np.shape(self.Par3d)[2]):
            Dataset.dataset_like(self.raster, self.Par3d[:, :, i], path=pnme[i])

__init__(raster: Dataset, no_parameters: int, no_lumped_par: int = 0, lumped_par_pos: list[int] | None = None, lake: bool = False, snow: bool = False, hru: bool = False, function: int = 1, k_upper_bound: int = 1, k_lower_bound: int = 50, muskingum: bool = False)

Initialize the Parameters class.

To initiate the Parameters class, you have to provide the Flow Acc raster.

Parameters:

Name	Type	Description	Default
raster	Dataset	A pyramids Dataset to get the spatial information of the catchment (DEM, flow accumulation or flow direction raster), read it using Dataset.read_file.	required
no_parameters	int	Number of parameters in the HBV model.	required
no_lumped_par	int	Number of lumped parameters. You have to enter the value of the lumped parameter at the end of the list. Defaults to 0 (no lumped parameters).	0
lumped_par_pos	list[int] | None	List of the order or position of lumped parameters among all the parameters of the lumped model (order starts from 0 to the length of the model parameters). Defaults to None (empty). The following order of parameters is used for the lumped HBV model: [ltt, utt, rfcf, sfcf, ttm, cfmax, cwh, cfr, fc, beta, e_corr, etf, lp, c_flux, k, k1, alpha, perc, pcorr, Kmuskingum, Xmuskingum].	None
lake	bool	True if there is a lake, False otherwise. Defaults to False.	False
snow	bool	True to run the snow-related processes, False otherwise. When True, parameters related to snow simulation have to be provided. Defaults to False.	False
hru	bool	True if the parameters will consider using HRUs. Defaults to False.	False
function	int	Which parameter-distribution strategy to bind to :attr:Function. One of 1 (:meth:par3d_lumped), 2 (:meth:par3d), 3 (:meth:par2d_lumped_k1_lake) or 4 (:meth:hydrologic_response_units). Defaults to 1. Any other value raises :class:ValueError. Note that hru=True overrides the choice with :meth:hydrologic_response_units regardless, but the selector is still validated so a typo is not masked.	1
k_upper_bound	int	Upper bound of K value (traveling time in muskingum routing method). Defaults to 1 hour.	1
k_lower_bound	int	Lower bound of K value (traveling time in muskingum routing method). Defaults to 50.	50
muskingum	bool	True if the routing function is muskingum. Defaults to False.	False

Raises:

Type	Description
TypeError	If raster is not a pyramids Dataset.
ValueError	If function is not one of the ints 1, 2, 3 or 4. A bool, a float such as 2.0, and an unhashable value are all rejected rather than coerced or allowed to raise TypeError.
AssertionError	If no_parameters is not an integer, if no_lumped_par is not an integer, or if the length of lumped_par_pos does not match no_lumped_par.
ValueError	If lumped_par_pos is not a list when no_lumped_par >= 1.

Note

Cells outside the catchment are identified by pyramids via read_array(masked=True) and stored as NaN in :attr:raster_array, which is promoted to floating point so it can hold them. no_elem, celli/cellj and the width of :attr:Par2d all derive from that mask, so the parameter vector length follows the raster's real domain.

Examples:

• Build the distributor from a small raster and inspect the domain it derived. The bottom-right cell is no-data, leaving three cells to parameterise:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from hapi.rrm.parameters import Parameters
  >>> raster = Dataset.create_from_array(
  ...     np.array([[1, 2], [3, -9999]], dtype="int32"),
  ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
  ...     no_data_value=-9999,
  ... )
  >>> distributor = Parameters(raster, 12)
  >>> distributor.no_elem
  3
  >>> distributor.Par2d.shape
  (12, 3)
  >>> list(zip(distributor.celli, distributor.cellj))
  [(0, 0), (0, 1), (1, 0)]
  

• A real value within 0.1% of the sentinel is kept, so it is treated as a catchment cell and widens the parameter array:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from hapi.rrm.parameters import Parameters
  >>> raster = Dataset.create_from_array(
  ...     np.array([[1, 2], [3, -9990]], dtype="int32"),
  ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
  ...     no_data_value=-9999,
  ... )
  >>> distributor = Parameters(raster, 12)
  >>> distributor.no_elem
  4
  >>> float(distributor.raster_array[1, 1])
  -9990.0
  

Source code in src/hapi/rrm/parameters.py

def __init__(
    self,
    raster: Dataset,
    no_parameters: int,
    no_lumped_par: int = 0,
    lumped_par_pos: list[int] | None = None,
    lake: bool = False,
    snow: bool = False,
    hru: bool = False,
    function: int = 1,
    k_upper_bound: int = 1,
    k_lower_bound: int = 50,
    muskingum: bool = False,
):
    """Initialize the Parameters class.

    To initiate the Parameters class, you have to provide the Flow Acc
    raster.

    Args:
        raster: A pyramids ``Dataset`` to get the spatial information
            of the catchment (DEM, flow accumulation or flow direction
            raster), read it using ``Dataset.read_file``.
        no_parameters: Number of parameters in the HBV model.
        no_lumped_par: Number of lumped parameters. You have to enter
            the value of the lumped parameter at the end of the list.
            Defaults to 0 (no lumped parameters).
        lumped_par_pos: List of the order or position of lumped
            parameters among all the parameters of the lumped model
            (order starts from 0 to the length of the model
            parameters). Defaults to None (empty). The following
            order of parameters is used for the lumped HBV model:
            [ltt, utt, rfcf, sfcf, ttm, cfmax, cwh, cfr, fc, beta,
            e_corr, etf, lp, c_flux, k, k1, alpha, perc, pcorr,
            Kmuskingum, Xmuskingum].
        lake: True if there is a lake, False otherwise.
            Defaults to False.
        snow: True to run the snow-related processes, False otherwise.
            When True, parameters related to snow simulation have to
            be provided. Defaults to False.
        hru: True if the parameters will consider using HRUs.
            Defaults to False.
        function: Which parameter-distribution strategy to bind to
            :attr:`Function`. One of ``1`` (:meth:`par3d_lumped`), ``2``
            (:meth:`par3d`), ``3`` (:meth:`par2d_lumped_k1_lake`) or ``4``
            (:meth:`hydrologic_response_units`). Defaults to 1. Any other value
            raises :class:`ValueError`. Note that ``hru=True`` overrides the choice
            with :meth:`hydrologic_response_units` regardless, but the selector is
            still validated so a typo is not masked.
        k_upper_bound: Upper bound of K value (traveling time in
            muskingum routing method). Defaults to 1 hour.
        k_lower_bound: Lower bound of K value (traveling time in
            muskingum routing method). Defaults to 50.
        muskingum: True if the routing function is muskingum.
            Defaults to False.

    Raises:
        TypeError: If `raster` is not a pyramids Dataset.
        ValueError: If `function` is not one of the ints 1, 2, 3 or 4. A `bool`,
            a `float` such as ``2.0``, and an unhashable value are all rejected
            rather than coerced or allowed to raise `TypeError`.
        AssertionError: If `no_parameters` is not an integer, if
            `no_lumped_par` is not an integer, or if the length of
            `lumped_par_pos` does not match `no_lumped_par`.
        ValueError: If `lumped_par_pos` is not a list when
            `no_lumped_par` >= 1.

    Note:
        Cells outside the catchment are identified by pyramids via
        ``read_array(masked=True)`` and stored as ``NaN`` in
        :attr:`raster_array`, which is promoted to floating point so it can hold
        them. ``no_elem``, ``celli``/``cellj`` and the width of :attr:`Par2d` all
        derive from that mask, so the parameter vector length follows the raster's
        real domain.

    Examples:
        - Build the distributor from a small raster and inspect the domain it
          derived. The bottom-right cell is no-data, leaving three cells to
          parameterise:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from hapi.rrm.parameters import Parameters
            >>> raster = Dataset.create_from_array(
            ...     np.array([[1, 2], [3, -9999]], dtype="int32"),
            ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
            ...     no_data_value=-9999,
            ... )
            >>> distributor = Parameters(raster, 12)
            >>> distributor.no_elem
            3
            >>> distributor.Par2d.shape
            (12, 3)
            >>> list(zip(distributor.celli, distributor.cellj))
            [(0, 0), (0, 1), (1, 0)]

            ```
        - A real value within 0.1% of the sentinel is kept, so it is treated as a
          catchment cell and widens the parameter array:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from hapi.rrm.parameters import Parameters
            >>> raster = Dataset.create_from_array(
            ...     np.array([[1, 2], [3, -9990]], dtype="int32"),
            ...     top_left_corner=(0.0, 8000.0), cell_size=4000.0, epsg=32618,
            ...     no_data_value=-9999,
            ... )
            >>> distributor = Parameters(raster, 12)
            >>> distributor.no_elem
            4
            >>> float(distributor.raster_array[1, 1])
            -9990.0

            ```
    """
    if lumped_par_pos is None:
        lumped_par_pos = []

    if not isinstance(raster, Dataset):
        raise TypeError(
            "raster should be a pyramids Dataset, read it using pyramids.dataset.Dataset.read_file"
        )
    assert isinstance(no_parameters, int), " no_parameters should be integer number"
    assert isinstance(no_lumped_par, int), (
        "no of lumped parameters should be integer"
    )

    if no_lumped_par >= 1:
        if isinstance(lumped_par_pos, list):
            assert no_lumped_par == len(lumped_par_pos), (
                f"you have to entered {no_lumped_par} no of lumped parameters but only {len(lumped_par_pos)} "
                f"position "
            )
        else:  # if not int or list
            raise ValueError(
                "you have one or more lumped parameters, so the position has to be entered as a list"
            )

    # Reject an unrecognised selector here rather than leaving `Function` unbound:
    # it is invoked on every calibration iteration, so a silent miss surfaces far
    # from the mistake as a bare AttributeError.
    strategies = {
        1: self.par3d_lumped,
        2: self.par3d,
        3: self.par2d_lumped_k1_lake,
        4: self.hydrologic_response_units,
    }
    # Test the type before the membership lookup: an unhashable selector (a list,
    # a set) raises TypeError from `in` rather than the documented ValueError, and
    # bool is a subclass of int, so True would otherwise silently select strategy 1.
    # A float is rejected outright rather than coerced -- 2.0 is a caller error, not
    # a spelling of 2.
    if isinstance(function, bool) or not isinstance(function, int):
        raise ValueError(
            f"function must be one of {sorted(strategies)}; got {function!r} of type "
            f"{type(function).__name__}. 1 = par3d_lumped, 2 = par3d, "
            "3 = par2d_lumped_k1_lake, 4 = hydrologic_response_units."
        )
    if function not in strategies:
        raise ValueError(
            f"function must be one of {sorted(strategies)}; got {function!r}. "
            "1 = par3d_lumped, 2 = par3d, 3 = par2d_lumped_k1_lake, "
            "4 = hydrologic_response_units."
        )
    self.Lake = lake
    self.Snow = snow
    self.no_lumped_par = no_lumped_par
    self.lumped_par_pos = lumped_par_pos
    self.HRUs = hru
    self.Kub = k_upper_bound
    self.Klb = k_lower_bound
    self.Maskingum = muskingum
    # read the raster
    self.raster = raster
    # No-data masking is delegated to pyramids: vectorised, dtype-aware, and it
    # also honours the band's GDAL mask band. Filling with NaN keeps the
    # float-array-with-NaN contract the rest of this class relies on.
    self.raster_array = np.ma.filled(
        raster.read_array(band=0, masked=True).astype(float), np.nan
    )
    # get the shape of the raster
    self.rows = raster.rows
    self.cols = raster.columns
    # get the no_value of in the raster
    self.noval = raster.no_data_value[0]
    if self.noval is None:
        warnings.warn(
            "the raster declares no no-data value, so every cell is treated as "
            "inside the catchment and the parameter vector is sized for the whole "
            "grid. If it has a sentinel, set it on the band.",
            UserWarning,
            stacklevel=2,
        )

    # count the number of non-empty cells
    if self.HRUs:
        self.values = list(
            set(
                [
                    int(self.raster_array[i, j])
                    for i in range(self.rows)
                    for j in range(self.cols)
                    if not np.isnan(self.raster_array[i, j])
                ]
            )
        )
        self.no_elem = len(self.values)
    else:
        # Count the cells the pyramids mask left intact (see
        # Catchment.read_flow_acc for why not Dataset.count_domain_cells).
        self.no_elem = int(np.count_nonzero(~np.isnan(self.raster_array)))

    self.no_parameters = no_parameters

    # store the indexes of the non-empty cells
    self.celli = []
    self.cellj = []
    for i in range(self.rows):
        for j in range(self.cols):
            if not np.isnan(self.raster_array[i, j]):
                self.celli.append(i)
                self.cellj.append(j)

    # create an empty 3D array [[raster dimension], no_parameters]
    self.Par3d = np.zeros([self.rows, self.cols, self.no_parameters]) * np.nan

    if no_lumped_par >= 1:
        # parameters in an array
        # remove a place for the lumped parameter (k1) lower zone coefficient
        self.no_parameters = self.no_parameters - no_lumped_par

    # all parameters lumped and distributed
    self.totnumberpar = self.no_parameters * self.no_elem + no_lumped_par
    # parameters in array
    # create a 2d array [no_parameters, no_cells]
    self.Par2d = np.zeros(
        shape=(self.no_parameters, self.no_elem), dtype=np.float32
    )

    self.Function = strategies[function]
    # to overwrite any choice user choose if the is HRUs
    if self.HRUs == 1:
        self.Function = self.hydrologic_response_units

    self.parameters_number()

    pass

calculate_k(x: float, position: int, upper_bound: float, lower_bound: float) -> float staticmethod

Calculate K parameter for Muskingum routing.

Takes the value of x parameter and generates 100 random values of the K parameter between the upper and lower constraints, then returns the value corresponding to the given position.

Parameters:

Name	Type	Description	Default
x	float	Weighting coefficient to determine the linearity of the water surface (one of the parameters of the Muskingum routing method).	required
position	int	Random position between upper and lower bounds of the K parameter.	required
upper_bound	float	Upper bound for the K parameter.	required
lower_bound	float	Lower bound for the K parameter.	required

Returns:

Type	Description
float	The K parameter value corresponding to the given position within the constrained range.

Source code in src/hapi/rrm/parameters.py

@staticmethod
def calculate_k(
    x: float, position: int, upper_bound: float, lower_bound: float
) -> float:
    """Calculate K parameter for Muskingum routing.

    Takes the value of x parameter and generates 100 random values of
    the K parameter between the upper and lower constraints, then
    returns the value corresponding to the given position.

    Args:
        x: Weighting coefficient to determine the linearity of the
            water surface (one of the parameters of the Muskingum
            routing method).
        position: Random position between upper and lower bounds of
            the K parameter.
        upper_bound: Upper bound for the K parameter.
        lower_bound: Lower bound for the K parameter.

    Returns:
        The K parameter value corresponding to the given position
            within the constrained range.
    """
    # k has to be smaller than this constraint
    constraint1 = 0.5 * 1 / (1 - x)
    # k has to be greater than this constraint
    constraint2 = 0.5 * 1 / x
    # if constraint is higher than UB take UB
    if constraint2 >= upper_bound:
        constraint2 = upper_bound
    # if constraint is lower than LB take UB
    if constraint1 <= lower_bound:
        constraint1 = lower_bound

    generated_k = np.linspace(constraint1, constraint2, 50)
    k = generated_k[int(round(position, 0))]
    return float(k)

hru_hand(dem: Dataset, flow_direction: Dataset, flow_path_length: Dataset, river: Dataset) -> tuple[np.ndarray, np.ndarray] staticmethod

Calculate Height Above Nearest Drainage (HAND) for HRU classification.

Calculates inputs for the HAND method for land use classification by tracing flow direction from each cell to the nearest river reach, then computing the elevation difference and the flow path distance.

Parameters:

Name	Type	Description	Default
dem	Dataset	A pyramids Dataset of the DEM raster.	required
flow_direction	Dataset	A pyramids Dataset of the flow direction raster.	required
flow_path_length	Dataset	A pyramids Dataset of the flow path length raster.	required
river	Dataset	A pyramids Dataset of the river location raster, where cells with value 1 indicate river presence.	required

Returns:

Type	Description
A tuple of two numpy ndarrays	hand: Height above nearest drainage for each cell. dist_to_nearest_drain: Distance to nearest drainage for each cell.

Raises:

Type	Description
ValueError	If the catchment boundaries contain anomalies (e.g., after cropping with a polygon).

Source code in src/hapi/rrm/parameters.py

@staticmethod
def hru_hand(
    dem: Dataset,
    flow_direction: Dataset,
    flow_path_length: Dataset,
    river: Dataset,
) -> tuple[np.ndarray, np.ndarray]:
    """Calculate Height Above Nearest Drainage (HAND) for HRU classification.

    Calculates inputs for the HAND method for land use
    classification by tracing flow direction from each cell to the
    nearest river reach, then computing the elevation difference
    and the flow path distance.

    Args:
        dem: A pyramids `Dataset` of the DEM raster.
        flow_direction: A pyramids `Dataset` of the flow direction
            raster.
        flow_path_length: A pyramids `Dataset` of the flow path
            length raster.
        river: A pyramids `Dataset` of the river location raster,
            where cells with value 1 indicate river presence.

    Returns:
        A tuple of two numpy ndarrays:
            - hand: Height above nearest drainage for each cell.
            - dist_to_nearest_drain: Distance to nearest drainage
              for each cell.

    Raises:
        ValueError: If the catchment boundaries contain anomalies
            (e.g., after cropping with a polygon).
    """
    # Use DEM raster information to run all loops
    dem_a = dem.read_array(band=0)
    no_val = np.float32(dem.no_data_value[0])
    rows = dem.rows
    cols = dem.columns

    # get the indices of the flow direction path
    fd_index = DEM(flow_direction.raster).flow_direction_index()

    # read the river location raster
    river_a = river.read_array(band=0)

    # read the flow path length raster
    fpl_a = flow_path_length.read_array(band=0)

    # trace the flow direction to the nearest river reach and store the location
    # of that nearest reach
    nearest_network = Parameters._trace_nearest_drainage(
        dem_a, no_val, river_a, fd_index, rows, cols
    )

    # the elevation difference to the nearest drainage cell is the height above
    # nearest drainage; the same difference over the flow path length raster is
    # the distance to that drainage
    hand = Parameters._difference_to_nearest_drainage(
        dem_a, dem_a, no_val, nearest_network, rows, cols
    )
    dist_to_nearest_drain = Parameters._difference_to_nearest_drainage(
        fpl_a, dem_a, no_val, nearest_network, rows, cols
    )

    return hand, dist_to_nearest_drain

hydrologic_response_units(par_g: list | np.ndarray)

Distribute parameters using Hydrologic Response Units (HRUs).

Takes a list of parameters (saved as one column or generated as a 1D list from an optimization algorithm) and distributes them horizontally on the number of cells given by a raster. The input raster should be a classified raster (by numbers) into classes to define the HRUs. Each HRU receives the same set of generated parameters.

Parameters:

Name	Type	Description	Default
par_g	list | ndarray	1D list or numpy array of parameters. For HRU without lumped parameters, the length should be no_elem * no_parameters. For HRU with lumped parameters, the lumped parameter values should be appended at the end.	required

Raises:

Type	Description
ValueError	If par_g is not a numpy ndarray or a list, or if the length of par_g does not match the expected number of parameters.
AssertionError	If there are lumped parameters and the length of par_g does not match the expected total.

Source code in src/hapi/rrm/parameters.py

def hydrologic_response_units(self, par_g: list | np.ndarray):  # ,kub=1,klb=0.5
    """Distribute parameters using Hydrologic Response Units (HRUs).

    Takes a list of parameters (saved as one column or generated as a
    1D list from an optimization algorithm) and distributes them
    horizontally on the number of cells given by a raster. The input
    raster should be a classified raster (by numbers) into classes to
    define the HRUs. Each HRU receives the same set of generated
    parameters.

    Args:
        par_g: 1D list or numpy array of parameters. For HRU without
            lumped parameters, the length should be
            ``no_elem * no_parameters``. For HRU with lumped
            parameters, the lumped parameter values should be
            appended at the end.

    Raises:
        ValueError: If `par_g` is not a numpy ndarray or a list, or
            if the length of `par_g` does not match the expected
            number of parameters.
        AssertionError: If there are lumped parameters and the length
            of `par_g` does not match the expected total.
    """
    # input data validation
    # data type
    if not (isinstance(par_g, np.ndarray) or isinstance(par_g, list)):
        raise ValueError("par_g should be of type 1d array or list")
    # assert isinstance(kub,numbers.Number) , " kub should be a number"
    # assert isinstance(klb,numbers.Number) , " klb should be a number"

    # input values
    if self.no_lumped_par > 0:
        par_no = (self.no_elem * self.no_parameters) + self.no_lumped_par
        assert len(par_g) == par_no, (
            f"As there is {self.no_lumped_par} lumped parameters, length of input parameters should be "
            f"{self.no_elem}*({self.no_parameters}-{self.no_lumped_par})+{self.no_lumped_par}="
            + str(
                self.no_elem * (self.no_parameters - self.no_lumped_par)
                + self.no_lumped_par
            )
            + f" not {len(par_g)} probably you have to add the value of the lumped parameter at the end of the list"
        )
    else:
        # if there is no lumped parameters
        if not len(par_g) == self.no_elem * self.no_parameters:
            raise ValueError(
                f"As there is no lumped parameters length of input parameters should be {self.no_elem}*"
                f"{self.no_parameters}={self.no_elem * self.no_parameters}"
            )

    # take the parameters from the generated parameters or the 1D list and
    # assign them to each cell
    self.Par2d = np.zeros(
        shape=(self.no_parameters, self.no_elem), dtype=np.float64
    )
    for i in range(self.no_elem):
        self.Par2d[:, i] = par_g[
            i * self.no_parameters : (i * self.no_parameters) + self.no_parameters
        ]

    # lumped parameters
    if self.no_lumped_par > 0:
        for i in range(self.no_lumped_par):
            # create a list with the value of the lumped parameter(k1)
            # (stored at the end of the list of the parameters)
            pk1 = (
                np.ones((1, self.no_elem))
                * par_g[(self.no_parameters * np.shape(self.Par2d)[1]) + i]
            )
            # put the list of parameter k1 at the 6 row
            self.Par2d = np.vstack(
                [
                    self.Par2d[: self.lumped_par_pos[i], :],
                    pk1,
                    self.Par2d[self.lumped_par_pos[i] :, :],
                ]
            )

    # calculate the value of k(travelling time in muskingum based on value of
    # x and the position and upper, lower bound of k value
    # for i in range(self.no_elem):
    #     self.Par2d[-2,i] = Parameters.calculateK(self.Par2d[-1,i],self.Par2d[-2,i],kub,klb)

    # assign the parameters from the array (no_parameters, no_cells) to
    # the spatially corrected location in par2d each soil type will have the same
    # generated parameters
    for i in range(self.no_elem):
        self.Par3d[self.raster_array == self.values[i]] = self.Par2d[:, i]

par2d_lumped_k1_lake(par_g: list | np.ndarray, no_parameters_lake: int)

Distribute parameters with a lumped K1 and lake parameters.

Takes a list of parameters and distributes them horizontally on the number of cells given by a raster. All parameters are distributed except the lower zone coefficient (K1), which is lumped and appended at the end of the parameter list. Lake parameters are extracted from the end of the parameter list.

Parameters:

Name	Type	Description	Default
par_g	list | ndarray	1D list or numpy array of parameters. Each cell's distributed parameters are listed sequentially, followed by the lumped K1 value, followed by lake parameters at the end. For example, with 14 cells and 11 distributed parameters: 14 * 11 = 154 + 1 (K1) = 155.	required
no_parameters_lake	int	Number of lake parameters to extract from the end of par_g.	required

Source code in src/hapi/rrm/parameters.py

def par2d_lumped_k1_lake(
    self, par_g: list | np.ndarray, no_parameters_lake: int
):  # ,kub,klb
    """Distribute parameters with a lumped K1 and lake parameters.

    Takes a list of parameters and distributes them horizontally on
    the number of cells given by a raster. All parameters are
    distributed except the lower zone coefficient (K1), which is
    lumped and appended at the end of the parameter list. Lake
    parameters are extracted from the end of the parameter list.

    Args:
        par_g: 1D list or numpy array of parameters. Each cell's
            distributed parameters are listed sequentially, followed
            by the lumped K1 value, followed by lake parameters at
            the end. For example, with 14 cells and 11 distributed
            parameters: ``14 * 11 = 154 + 1 (K1) = 155``.
        no_parameters_lake: Number of lake parameters to extract
            from the end of `par_g`.
    """
    # parameters in array
    # remove a place for the lumped parameter (k1) lower zone coefficient
    no_parameters = self.no_parameters - 1

    # create a 2d array [no_parameters, no_cells]
    self.Par2d = np.ones((no_parameters, self.no_elem))  # type: ignore[assignment]

    # take the parameters from the generated parameters or the 1D list and
    # assign them to each cell
    for i in range(self.no_elem):
        self.Par2d[:, i] = par_g[
            i * no_parameters : (i * no_parameters) + no_parameters
        ]

    # create a list with the value of the lumped parameter(k1)
    # (stored at the end of the list of the parameters)
    pk1 = (
        np.ones((1, self.no_elem))
        * par_g[(np.shape(self.Par2d)[0] * np.shape(self.Par2d)[1])]
    )

    # put the list of parameter k1 at the 6 row
    self.Par2d = np.vstack([self.Par2d[:6, :], pk1, self.Par2d[6:, :]])

    # assign the parameters from the array (no_parameters, no_cells) to
    # the spatially corrected location in par2d
    for i in range(self.no_elem):
        self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

    # calculate the value of k(travelling time in muskingum based on value of
    # x and the position and upper, lower bound of k value
    # for i in range(self.no_elem):
    #     self.Par3d[self.celli[i],self.cellj[i],-2] = Parameters.calculateK(
    #     self.Par3d[self.celli[i],self.cellj[i],-1],self.Par3d[self.celli[i],self.cellj[i],-2],kub,klb)

    # lake parameters
    self.lake_par = par_g[len(par_g) - no_parameters_lake :]

par3d(par_g: list | np.ndarray)

Distribute parameters horizontally across grid cells.

Takes a list of parameters (saved as one column or generated as a 1D list from an optimization algorithm) and distributes them horizontally on the number of cells given by a raster.

Parameters:

Name	Type	Description	Default
par_g	list | ndarray	1D list or numpy array of parameters. For totally distributed parameters, the length should be no_elem * no_parameters. For lumped parameters, the lumped parameter values should be appended at the end.	required

Raises:

Type	Description
AssertionError	If the length of par_g does not match the expected number of parameters based on the number of elements and lumped parameters.

Source code in src/hapi/rrm/parameters.py

def par3d(self, par_g: list | np.ndarray):  # , kub=1,klb=0.5, Maskingum=True
    """Distribute parameters horizontally across grid cells.

    Takes a list of parameters (saved as one column or generated as a
    1D list from an optimization algorithm) and distributes them
    horizontally on the number of cells given by a raster.

    Args:
        par_g: 1D list or numpy array of parameters. For totally
            distributed parameters, the length should be
            ``no_elem * no_parameters``. For lumped parameters, the
            lumped parameter values should be appended at the end.

    Raises:
        AssertionError: If the length of `par_g` does not match the
            expected number of parameters based on the number of
            elements and lumped parameters.
    """
    # input data validation
    # data type
    # assert type(par_g)==np.ndarray or type(par_g)==list, "par_g should be of type 1d array or list"
    # assert isinstance(kub,numbers.Number) , " kub should be a number"
    # assert isinstance(klb,numbers.Number) , " klb should be a number"

    # input values
    if self.no_lumped_par > 0:
        par_no = (self.no_elem * self.no_parameters) + self.no_lumped_par

        assert len(par_g) == par_no, (
            f"As there is {self.no_lumped_par} lumped parameters, length of input parameters should be "
            f"{self.no_elem}"
            + f"*({self.no_parameters + self.no_lumped_par} - {self.no_lumped_par}) + {self.no_lumped_par} = "
            + f"{self.no_elem * (self.no_parameters - self.no_lumped_par) + self.no_lumped_par} not {len(par_g)}"
            + " probably you have to add the value of the lumped parameter at the end of the list"
        )
    else:
        # if there are no lumped parameters
        par_no = self.no_elem * self.no_parameters
        assert len(par_g) == par_no, (
            f"As there is no lumped parameters length of input parameters should be {self.no_elem} * "
            + f"{self.no_parameters} = {self.no_elem * self.no_parameters}"
        )

    # parameters in array
    # create a 2d array [no_parameters, no_cells]
    self.Par2d = np.ones((self.no_parameters, self.no_elem))  # type: ignore[assignment]
    # take the parameters from the generated parameters or the 1D list and
    # assign them to each cell
    for i in range(self.no_elem):
        self.Par2d[:, i] = par_g[
            i * self.no_parameters : (i * self.no_parameters) + self.no_parameters
        ]

    # lumped parameters
    if self.no_lumped_par > 0:
        for i in range(self.no_lumped_par):
            # create a list with the value of the lumped parameter(k1)
            # (stored at the end of the list of the parameters)
            pk1 = (
                np.ones((1, self.no_elem))
                * par_g[(self.no_parameters * np.shape(self.Par2d)[1]) + i]
            )
            # put the list of parameter k1 at the 6th row.
            self.Par2d = np.vstack(
                [
                    self.Par2d[: self.lumped_par_pos[i], :],
                    pk1,
                    self.Par2d[self.lumped_par_pos[i] :, :],
                ]
            )

    # assign the parameters from the array (no_parameters, no_cells) to
    # the spatially corrected location in par2d
    for i in range(self.no_elem):
        self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

par3d_lumped(par_g: list | np.ndarray)

Distribute lumped parameters horizontally across grid cells.

Takes a list of parameters (saved as one column or generated as a 1D list from an optimization algorithm) and distributes them horizontally on the number of cells given by a raster, where all parameters are lumped (same value for every cell).

Parameters:

Name	Type	Description	Default
par_g	list | ndarray	1D list or numpy array of lumped parameters. The length should equal no_parameters.	required

Raises:

Type	Description
ValueError	If par_g is not a numpy ndarray or a list.

Source code in src/hapi/rrm/parameters.py

def par3d_lumped(
    self, par_g: list | np.ndarray
):  # , kub=1, klb=0.5, Maskingum = True
    r"""Distribute lumped parameters horizontally across grid cells.

    Takes a list of parameters (saved as one column or generated as a
    1D list from an optimization algorithm) and distributes them
    horizontally on the number of cells given by a raster, where all
    parameters are lumped (same value for every cell).

    Args:
        par_g: 1D list or numpy array of lumped parameters.
            The length should equal ``no_parameters``.

    Raises:
        ValueError: If `par_g` is not a numpy ndarray or a list.
    """
    # input data validation
    # data type
    if not (isinstance(par_g, np.ndarray) or isinstance(par_g, list)):
        raise ValueError("par_g should be of type 1d array or list")
    # assert isinstance(kub,numbers.Number) , " kub should be a number"
    # assert isinstance(klb,numbers.Number) , " klb should be a number"

    # take the parameters from the generated parameters or the 1D list and
    # assign them to each cell
    for i in range(self.no_elem):
        self.Par2d[:, i] = par_g

    # assign the parameters from the array (no_parameters, no_cells) to
    # the spatially corrected location in par2d
    for i in range(self.no_elem):
        self.Par3d[self.celli[i], self.cellj[i], :] = self.Par2d[:, i]

parameters_number()

Calculate the total number of parameters for the optimization.

Calculates the number of parameters that the optimization algorithm will search for. Use this only in case of totally distributed catchment parameters. In case of lumped parameters, the number of parameters is the same as the number of parameters of the conceptual model.

The result is stored in the ParametersNO attribute.

Note

The Parameters object should have the following attributes before calling this method: raster, no_parameters, no_lumped_par, and HRUs.

Source code in src/hapi/rrm/parameters.py

def parameters_number(self):
    """Calculate the total number of parameters for the optimization.

    Calculates the number of parameters that the optimization
    algorithm will search for. Use this only in case of totally
    distributed catchment parameters. In case of lumped parameters,
    the number of parameters is the same as the number of parameters
    of the conceptual model.

    The result is stored in the ``ParametersNO`` attribute.

    Note:
        The Parameters object should have the following attributes
        before calling this method: ``raster``, ``no_parameters``,
        ``no_lumped_par``, and ``HRUs``.
    """
    if not self.HRUs:
        if self.no_lumped_par > 0:
            # self.ParametersNO = (self.no_elem *( self.no_parameters - self.no_lumped_par)) + self.no_lumped_par
            self.ParametersNO = (
                self.no_elem * self.no_parameters
            ) + self.no_lumped_par
        else:
            # if there is no lumped parameters
            self.ParametersNO = self.no_elem * self.no_parameters
    else:
        if self.no_lumped_par > 0:
            # self.ParametersNO = (self.no_elem * (self.no_parameters - self.no_lumped_par)) + self.no_lumped_par
            self.ParametersNO = (
                self.no_elem * self.no_parameters
            ) + self.no_lumped_par
        else:
            # if there is no lumped parameters
            self.ParametersNO = self.no_elem * self.no_parameters

save_parameters(path: str | None)

Save distributed parameters as raster files.

Takes the generated 3D parameter array and saves each parameter layer as a separate GeoTIFF raster file.

Parameters:

Name	Type	Description	Default
path	str | None	Path to the folder where the parameter rasters will be saved.	required

Raises:

Type	Description
TypeError	path is not a str. Output names are built by concatenation (path + name), which a Path does not support.
FileNotFoundError	The output directory does not exist. Checked up front so the failure does not surface midway through writing.

Note

The Parameters object should have the following attributes set before calling this method: DistParFn, raster, Par, no_parameters, snow, kub, and klb.

Source code in src/hapi/rrm/parameters.py

def save_parameters(self, path: str | None):
    """Save distributed parameters as raster files.

    Takes the generated 3D parameter array and saves each parameter
    layer as a separate GeoTIFF raster file.

    Args:
        path: Path to the folder where the parameter rasters will
            be saved.

    Raises:
        TypeError: `path` is not a ``str``. Output names are built by concatenation
            (``path + name``), which a ``Path`` does not support.
        FileNotFoundError: The output directory does not exist. Checked up front so
            the failure does not surface midway through writing.

    Note:
        The Parameters object should have the following attributes
        set before calling this method: ``DistParFn``, ``raster``,
        ``Par``, ``no_parameters``, ``snow``, ``kub``, and ``klb``.
    """
    # Not delegated to pyramids: the output names are built by string concatenation
    # below (`path + name`), so this genuinely needs a str, and the directory must
    # exist before the first raster is written. Raised rather than asserted so the
    # checks survive `python -O`.
    if not isinstance(path, str):
        raise TypeError(f"path should be of type string, given: {type(path)}")
    if not os.path.exists(path):
        raise FileNotFoundError(f"{path} you have provided does not exist")

    # save
    if self.Snow == 0:  # now snow subroutine
        pnme = [
            "01_rfcf",
            "02_FC",
            "03_BETA",
            "04_ETF",
            "05_LP",
            "06_K0",
            "07_K1",
            "08_K2",
            "09_UZL",
            "10_PERC",
            "11_Kmuskingum",
            "12_Xmuskingum",
        ]
    else:  # there is snow subtoutine
        pnme = [
            "01_ltt",
            "02_rfcf",
            "03_sfcf",
            "04_cfmax",
            "05_cwh",
            "06_cfr",
            "07_fc",
            "08_beta",
            "09_etf",
            "10_lp",
            "11_k0",
            "12_k1",
            "13_k2",
            "14_uzl",
            "18_perc",
        ]

    if path is not None:
        pnme = [
            path + i + "_" + str(dt.datetime.now())[0:10] + ".tif" for i in pnme
        ]

    for i in range(np.shape(self.Par3d)[2]):
        Dataset.dataset_like(self.raster, self.Par3d[:, :, i], path=pnme[i])