Bank Erosion Module #

The Bank Erosion module is responsible for calculating bank erosion based on hydrodynamic data and detected bank lines. It is one of the core components of the D-FAST Bank Erosion software.

Overview #

The Bank Erosion module calculates the amount of bank material that will be eroded during the first year and until equilibrium, based on hydrodynamic simulation results and detected bank lines. It takes into account various factors such as river geometry, discharge levels, and shipping parameters.

Components #

The Bank Erosion module consists of the following components:

Main Classes #

`dfastbe.bank_erosion.bank_erosion` #

This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation version 2.1.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with this library; if not, see http://www.gnu.org/licenses/.

contact: delft3d.support@deltares.nl Stichting Deltares P.O. Box 177 2600 MH Delft, The Netherlands

INFORMATION This file is part of D-FAST Bank Erosion: https://github.com/Deltares/D-FAST_Bank_Erosion

`Erosion` #

Class to handle the bank erosion calculations.

Source code in src/dfastbe/bank_erosion/bank_erosion.py

class Erosion:
    """Class to handle the bank erosion calculations."""

    def __init__(self, config_file: ConfigFile, gui: bool = False):
        """Initialize the Erosion class."""
        self.root_dir = config_file.root_dir
        self._config_file = config_file
        self.gui = gui

        self.river_data = ErosionRiverData(config_file)
        self.river_center_line_arr = self.river_data.river_center_line.as_array()
        self.simulation_data = self.river_data.simulation_data()
        self.sim_files, self.p_discharge = self.river_data.get_erosion_sim_data(
            self.river_data.num_discharge_levels
        )
        self.bl_processor = BankLinesProcessor(self.river_data)
        self.debugger = Debugger(config_file.crs, self.river_data.output_dir)
        self.erosion_calculator = ErosionCalculator()

    @property
    def config_file(self) -> ConfigFile:
        """Configuration file object."""
        return self._config_file

    def get_ship_parameters(
        self, num_stations_per_bank: List[int]
    ) -> Dict[str, List[np.ndarray]]:
        """Get ship parameters from the configuration file."""
        ship_relative_velocity = self.config_file.get_parameter(
            "Erosion", "VShip", num_stations_per_bank, positive=True, onefile=True
        )
        num_ships_year = self.config_file.get_parameter(
            "Erosion", "NShip", num_stations_per_bank, positive=True, onefile=True
        )
        num_waves_p_ship = self.config_file.get_parameter(
            "Erosion",
            "NWave",
            num_stations_per_bank,
            default=5,
            positive=True,
            onefile=True,
        )
        ship_draught = self.config_file.get_parameter(
            "Erosion", "Draught", num_stations_per_bank, positive=True, onefile=True
        )
        ship_type = self.config_file.get_parameter(
            "Erosion", "ShipType", num_stations_per_bank, valid=[1, 2, 3], onefile=True
        )
        parslope0 = self.config_file.get_parameter(
            "Erosion",
            "Slope",
            num_stations_per_bank,
            default=20,
            positive=True,
            ext="slp",
        )
        reed_wave_damping_coeff = self.config_file.get_parameter(
            "Erosion",
            "Reed",
            num_stations_per_bank,
            default=0,
            positive=True,
            ext="rdd",
        )

        ship_data = {
            "vship0": ship_relative_velocity,
            "Nship0": num_ships_year,
            "nwave0": num_waves_p_ship,
            "Tship0": ship_draught,
            "ship0": ship_type,
            "parslope0": parslope0,
            "parreed0": reed_wave_damping_coeff,
        }
        return ship_data

    def _process_river_axis_by_center_line(self) -> LineGeometry:
        """Process the river axis by the center line.

        Intersect the river center line with the river axis to map the stations from the first to the latter
        then clip the river axis by the first and last station of the centerline.
        """
        river_axis = LineGeometry(self.river_data.river_axis, crs=self.config_file.crs)
        river_axis_numpy = river_axis.as_array()
        # optional sorting --> see 04_Waal_D3D example
        # check: sum all distances and determine maximum distance ...
        # if maximum > alpha * sum then perform sort
        # Waal OK: 0.0082 ratio max/sum, Waal NotOK: 0.13 - Waal: 2500 points,
        # so even when OK still some 21 times more than 1/2500 = 0.0004
        dist2 = (np.diff(river_axis_numpy, axis=0) ** 2).sum(axis=1)
        alpha = dist2.max() / dist2.sum()
        if alpha > 0.03:
            print("The river axis needs sorting!!")

        # map km to axis points, further using axis
        log_text("chainage_to_axis")
        river_axis_km = river_axis.intersect_with_line(self.river_center_line_arr)

        # clip river axis to reach of interest (get the closest point to the first and last station)
        i1 = np.argmin(
            ((self.river_center_line_arr[0, :2] - river_axis_numpy) ** 2).sum(axis=1)
        )
        i2 = np.argmin(
            ((self.river_center_line_arr[-1, :2] - river_axis_numpy) ** 2).sum(axis=1)
        )
        if i1 < i2:
            river_axis_km = river_axis_km[i1 : i2 + 1]
            river_axis_numpy = river_axis_numpy[i1 : i2 + 1]
        else:
            # reverse river axis
            river_axis_km = river_axis_km[i2 : i1 + 1][::-1]
            river_axis_numpy = river_axis_numpy[i2 : i1 + 1][::-1]

        # river_axis = LineString(river_axis_numpy)
        river_axis = LineGeometry(river_axis_numpy, crs=self.config_file.crs)
        river_axis.add_data(data={"stations": river_axis_km})
        return river_axis

    def _get_fairway_data(
        self,
        river_axis: LineGeometry,
        mesh_data: MeshData,
    ):
        # map km to fairway points, further using axis
        log_text("chainage_to_fairway")
        # intersect fairway and mesh
        fairway_intersection_coords, fairway_face_indices = intersect_line_mesh(
            river_axis.as_array(), mesh_data
        )
        if self.river_data.debug:
            arr = (
                fairway_intersection_coords[:-1] + fairway_intersection_coords[1:]
            ) / 2
            line_geom = LineGeometry(arr, crs=self.config_file.crs)
            line_geom.to_file(
                file_name=f"{str(self.river_data.output_dir)}{os.sep}fairway_face_indices.shp",
                data={"iface": fairway_face_indices},
            )

        return FairwayData(fairway_face_indices, fairway_intersection_coords)

    def calculate_fairway_bank_line_distance(
        self,
        bank_data: BankData,
        fairway_data: FairwayData,
        simulation_data: ErosionSimulationData,
    ):
        """Map bank data to fairway data.

        Args:
            bank_data (BankData):
            fairway_data (FairwayData):
            simulation_data (ErosionSimulationData):

        Returns:
            FairwayData:
                The method updates the following attributes in the `bank_data` instance
                    - fairway_face_indices
                    - fairway_distances
            BankData:
                the following attributes in the `fairway_data` instance
                    - fairway_initial_water_levels
        """
        # distance fairway-bankline (bank-fairway)
        log_text("bank_distance_fairway")

        num_fairway_face_ind = len(fairway_data.fairway_face_indices)

        for bank_i, single_bank in enumerate(bank_data):
            bank_coords = single_bank.bank_line_coords
            coords_mid = (bank_coords[:-1] + bank_coords[1:]) / 2
            bank_fairway_dist = np.zeros(len(coords_mid))
            bp_fw_face_idx = np.zeros(len(coords_mid), dtype=int)

            for ind, coord_i in enumerate(coords_mid):
                # find closest fairway support node
                closest_ind = np.argmin(
                    ((coord_i - fairway_data.intersection_coords) ** 2).sum(axis=1)
                )
                fairway_coord = fairway_data.intersection_coords[closest_ind]
                fairway_bank_distance = ((coord_i - fairway_coord) ** 2).sum() ** 0.5
                # If fairway support node is also the closest projected fairway point, then it likely
                # that that point is one of the original support points (a corner) of the fairway path
                # and located inside a grid cell. The segments before and after that point will then
                # both be located inside that same grid cell, so let's pick the segment before the point.
                # If the point happens to coincide with a grid edge and the two segments are located
                # in different grid cells, then we could either simply choose one or add complexity to
                # average the values of the two grid cells. Let's go for the simplest approach ...
                iseg = max(closest_ind - 1, 0)
                if closest_ind > 0:
                    alpha = calculate_alpha(
                        fairway_data.intersection_coords,
                        closest_ind,
                        closest_ind - 1,
                        coord_i,
                    )
                    if 0 < alpha < 1:
                        fwp1 = fairway_data.intersection_coords[
                            closest_ind - 1
                        ] + alpha * (
                            fairway_data.intersection_coords[closest_ind]
                            - fairway_data.intersection_coords[closest_ind - 1]
                        )
                        d1 = ((coord_i - fwp1) ** 2).sum() ** 0.5
                        if d1 < fairway_bank_distance:
                            fairway_bank_distance = d1
                            # projected point located on segment before, which corresponds to initial choice: iseg = ifw - 1
                if closest_ind < num_fairway_face_ind:
                    alpha = calculate_alpha(
                        fairway_data.intersection_coords,
                        closest_ind + 1,
                        closest_ind,
                        coord_i,
                    )
                    if 0 < alpha < 1:
                        fwp1 = fairway_data.intersection_coords[closest_ind] + alpha * (
                            fairway_data.intersection_coords[closest_ind + 1]
                            - fairway_data.intersection_coords[closest_ind]
                        )
                        d1 = ((coord_i - fwp1) ** 2).sum() ** 0.5
                        if d1 < fairway_bank_distance:
                            fairway_bank_distance = d1
                            iseg = closest_ind

                bp_fw_face_idx[ind] = fairway_data.fairway_face_indices[iseg]
                bank_fairway_dist[ind] = fairway_bank_distance

            if self.river_data.debug:
                line_geom = LineGeometry(coords_mid, crs=self.config_file.crs)
                line_geom.to_file(
                    file_name=f"{self.river_data.output_dir}/bank_{bank_i + 1}_chainage_and_fairway_face_idx.shp",
                    data={
                        "chainage": single_bank.bank_chainage_midpoints,
                        "iface_fw": bp_fw_face_idx[bank_i],
                    },
                )

            single_bank.fairway_face_indices = bp_fw_face_idx
            single_bank.fairway_distances = bank_fairway_dist

        # water level at fairway
        water_level_fairway_ref = []
        for single_bank in bank_data:
            ii = single_bank.fairway_face_indices
            water_level_fairway_ref.append(simulation_data.water_level_face[ii])
        fairway_data.fairway_initial_water_levels = water_level_fairway_ref

    def _prepare_initial_conditions(
        self,
        config_file: ConfigFile,
        num_stations_per_bank: List[int],
        fairway_data: FairwayData,
    ) -> ErosionInputs:
        # wave reduction s0, s1
        wave_fairway_distance_0 = config_file.get_parameter(
            "Erosion",
            "Wave0",
            num_stations_per_bank,
            default=200,
            positive=True,
            onefile=True,
        )
        wave_fairway_distance_1 = config_file.get_parameter(
            "Erosion",
            "Wave1",
            num_stations_per_bank,
            default=150,
            positive=True,
            onefile=True,
        )

        # save 1_banklines
        # read vship, nship, nwave, draught (tship), shiptype ... independent of level number
        shipping_data = self.get_ship_parameters(num_stations_per_bank)

        # read classes flag (yes: banktype = taucp, no: banktype = tauc) and banktype (taucp: 0-4 ... or ... tauc = critical shear value)
        classes = config_file.get_bool("Erosion", "Classes")
        if classes:
            bank_type = config_file.get_parameter(
                "Erosion",
                "BankType",
                num_stations_per_bank,
                default=0,
                ext=".btp",
            )
            tauc = []
            for bank in bank_type:
                tauc.append(ErosionInputs.taucls[bank])
        else:
            tauc = config_file.get_parameter(
                "Erosion",
                "BankType",
                num_stations_per_bank,
                default=0,
                ext=".btp",
            )
            thr = (ErosionInputs.taucls[:-1] + ErosionInputs.taucls[1:]) / 2
            bank_type = [None] * len(thr)
            for ib, shear_stress in enumerate(tauc):
                bt = np.zeros(shear_stress.size)
                for thr_i in thr:
                    bt[shear_stress < thr_i] += 1
                bank_type[ib] = bt

        # read bank protection level dike_height
        zss_miss = -1000
        dike_height = config_file.get_parameter(
            "Erosion",
            "ProtectionLevel",
            num_stations_per_bank,
            default=zss_miss,
            ext=".bpl",
        )
        # if dike_height undefined, set dike_height equal to water_level_fairway_ref - 1
        for ib, one_zss in enumerate(dike_height):
            mask = one_zss == zss_miss
            one_zss[mask] = fairway_data.fairway_initial_water_levels[ib][mask] - 1

        data = dict(
            wave_fairway_distance_0=wave_fairway_distance_0,
            wave_fairway_distance_1=wave_fairway_distance_1,
            bank_protection_level=dike_height,
            tauc=tauc,
        )
        return ErosionInputs.from_column_arrays(
            data, SingleErosion, shipping_data=shipping_data, bank_type=bank_type
        )

    def _calculate_bank_height(self, bank_data: BankData, simulation_data: ErosionSimulationData) -> BankData:
        # bank height = maximum bed elevation per cell
        for bank_i in bank_data:
            bank_i.height = simulation_data.calculate_bank_height(bank_i, self.river_data.zb_dx)

        return bank_data

    def _process_discharge_levels(
        self,
        km_mid,
        km_bin,
        config_file: ConfigFile,
        erosion_inputs: ErosionInputs,
        bank_data: BankData,
        fairway_data: FairwayData,
    ) -> Tuple[WaterLevelData, ErosionResults]:

        num_levels = self.river_data.num_discharge_levels
        num_km = len(km_mid)

        # initialize arrays for erosion loop over all discharges
        discharge_levels = []

        log_text("total_time", data={"t": self.river_data.erosion_time})

        for level_i in range(num_levels):
            log_text(
                "discharge_header",
                data={
                    "i": level_i + 1,
                    "p": self.p_discharge[level_i],
                    "t": self.p_discharge[level_i] * self.river_data.erosion_time,
                },
            )

            log_text("read_q_params", indent="  ")
            # 1) read level-specific parameters
            # read ship_velocity, num_ship, nwave, draught, ship_type, slope, reed, fairway_depth, ... (level specific values)
            level_parameters = self._read_discharge_parameters(
                level_i, erosion_inputs.shipping_data, bank_data.num_stations_per_bank
            )

            # 2) load FM result
            log_text("-", indent="  ")
            log_text(
                "read_simdata", data={"file": self.sim_files[level_i]}, indent="  "
            )
            simulation_data = ErosionSimulationData.read(
                self.sim_files[level_i], indent="  "
            )
            log_text("bank_erosion", indent="  ")

            if level_i == 0:
                bank_data = self._calculate_bank_height(bank_data, simulation_data)

            single_level, dvol_bank = self.compute_erosion_per_level(
                level_i,
                bank_data,
                simulation_data,
                fairway_data,
                level_parameters,
                erosion_inputs,
                km_bin,
                num_km,
            )

            discharge_levels.append(single_level)

            error_vol_file = config_file.get_str(
                "Erosion", f"EroVol{level_i + 1}", default=f"erovolQ{level_i + 1}.evo"
            )
            log_text("save_error_vol", data={"file": error_vol_file}, indent="  ")
            write_km_eroded_volumes(
                km_mid, dvol_bank, f"{self.river_data.output_dir}/{error_vol_file}"
            )

        # shape is (num_levels, 2, (num_stations_per_bank))
        # if num_levels = 13 and the num_stations_per_bank = [10, 15]
        # then shape = (13, 2, (10, 15)) list of 13 elements, each element is a list of 2 elements
        # first an array of 10 elements, and the second is array of 15 elements
        discharge_levels = DischargeLevels(discharge_levels)
        flow_erosion_dist = discharge_levels.accumulate("erosion_distance_flow")
        ship_erosion_dist = discharge_levels.accumulate("erosion_distance_shipping")
        total_erosion_dist = discharge_levels.accumulate("erosion_distance_tot")
        total_eroded_vol = discharge_levels.accumulate("erosion_volume_tot")

        erosion_results = ErosionResults(
            erosion_time=self.river_data.erosion_time,
            flow_erosion_dist=flow_erosion_dist,
            ship_erosion_dist=ship_erosion_dist,
            total_erosion_dist=total_erosion_dist,
            total_eroded_vol=total_eroded_vol,
            eq_erosion_dist=discharge_levels._get_attr_both_sides_level(
                "erosion_distance_eq", num_levels - 1
            ),
            eq_eroded_vol=discharge_levels._get_attr_both_sides_level(
                "erosion_volume_eq", num_levels - 1
            ),
        )
        water_level_data = discharge_levels.get_water_level_data()

        bank_data.left.bank_line_size, bank_data.right.bank_line_size = (
            bank_data.left.segment_length,
            bank_data.right.segment_length,
        )

        return water_level_data, erosion_results

    def _postprocess_erosion_results(
        self,
        km_bin: Tuple[float, float, float],
        km_mid,
        bank_data: BankData,
        erosion_results: ErosionResults,
    ) -> Tuple[List[LineString], List[LineString], List[LineString]]:
        """Postprocess the erosion results to get the new bank lines and volumes."""
        log_text("=")
        avg_erosion_rate = np.zeros(bank_data.n_bank_lines)
        dn_max = np.zeros(bank_data.n_bank_lines)
        d_nav_flow = np.zeros(bank_data.n_bank_lines)
        d_nav_ship = np.zeros(bank_data.n_bank_lines)
        d_nav_eq = np.zeros(bank_data.n_bank_lines)
        dn_max_eq = np.zeros(bank_data.n_bank_lines)
        eq_eroded_vol_per_km = np.zeros((len(km_mid), bank_data.n_bank_lines))
        total_eroded_vol_per_km = np.zeros((len(km_mid), bank_data.n_bank_lines))
        xy_line_new_list = []
        bankline_new_list = []
        xy_line_eq_list = []
        bankline_eq_list = []
        for ib, single_bank in enumerate(bank_data):
            bank_coords = single_bank.bank_line_coords
            avg_erosion_rate[ib] = (
                erosion_results.total_erosion_dist[ib] * single_bank.bank_line_size
            ).sum() / single_bank.bank_line_size.sum()
            dn_max[ib] = erosion_results.total_erosion_dist[ib].max()
            d_nav_flow[ib] = (
                erosion_results.flow_erosion_dist[ib] * single_bank.bank_line_size
            ).sum() / single_bank.bank_line_size.sum()
            d_nav_ship[ib] = (
                erosion_results.ship_erosion_dist[ib] * single_bank.bank_line_size
            ).sum() / single_bank.bank_line_size.sum()
            d_nav_eq[ib] = (
                erosion_results.eq_erosion_dist[ib] * single_bank.bank_line_size
            ).sum() / single_bank.bank_line_size.sum()
            dn_max_eq[ib] = erosion_results.eq_erosion_dist[ib].max()
            log_text("bank_dnav", data={"ib": ib + 1, "v": avg_erosion_rate[ib]})
            log_text("bank_dnavflow", data={"v": d_nav_flow[ib]})
            log_text("bank_dnavship", data={"v": d_nav_ship[ib]})
            log_text("bank_dnmax", data={"v": dn_max[ib]})
            log_text("bank_dnaveq", data={"v": d_nav_eq[ib]})
            log_text("bank_dnmaxeq", data={"v": dn_max_eq[ib]})

            xy_line_new = move_line(
                bank_coords,
                erosion_results.total_erosion_dist[ib],
                single_bank.is_right_bank,
            )
            xy_line_new_list.append(xy_line_new)
            bankline_new_list.append(LineString(xy_line_new))

            xy_line_eq = move_line(
                bank_coords,
                erosion_results.eq_erosion_dist[ib],
                single_bank.is_right_bank,
            )
            xy_line_eq_list.append(xy_line_eq)
            bankline_eq_list.append(LineString(xy_line_eq))

            dvol_eq = get_km_eroded_volume(
                single_bank.bank_chainage_midpoints,
                erosion_results.eq_eroded_vol[ib],
                km_bin,
            )
            eq_eroded_vol_per_km[:, ib] = dvol_eq
            dvol_tot = get_km_eroded_volume(
                single_bank.bank_chainage_midpoints,
                erosion_results.total_eroded_vol[ib],
                km_bin,
            )
            total_eroded_vol_per_km[:, ib] = dvol_tot
            if ib < bank_data.n_bank_lines - 1:
                log_text("-")

        erosion_results.avg_erosion_rate = avg_erosion_rate
        erosion_results.eq_eroded_vol_per_km = eq_eroded_vol_per_km
        erosion_results.total_eroded_vol_per_km = total_eroded_vol_per_km

        return bankline_new_list, bankline_eq_list, xy_line_eq_list

    def _get_param(
        self, name: str, default_val, iq_str, num_stations_per_bank, **kwargs
    ):
        return self.config_file.get_parameter(
            "Erosion",
            f"{name}{iq_str}",
            num_stations_per_bank,
            default=default_val,
            **kwargs,
        )

    def _read_discharge_parameters(
        self,
        level_i: int,
        shipping_data: Dict[str, List[np.ndarray]],
        num_stations_per_bank: List[int],
    ) -> SingleLevelParameters:
        """Read Discharge level parameters.

        Read all discharge-specific input arrays for level *iq*.
        Returns a dict with keys: vship, num_ship, n_wave, t_ship, ship_type,
        mu_slope, mu_reed, par_slope, par_reed.
        """
        iq_str = f"{level_i + 1}"

        ship_velocity = self._get_param(
            "VShip",
            shipping_data["vship0"],
            iq_str,
            num_stations_per_bank,
        )
        num_ship = self._get_param(
            "NShip",
            shipping_data["Nship0"],
            iq_str,
            num_stations_per_bank,
        )
        num_waves_per_ship = self._get_param(
            "NWave",
            shipping_data["nwave0"],
            iq_str,
            num_stations_per_bank,
        )
        ship_draught = self._get_param(
            "Draught",
            shipping_data["Tship0"],
            iq_str,
            num_stations_per_bank,
        )
        ship_type = self._get_param(
            "ShipType",
            shipping_data["ship0"],
            iq_str,
            num_stations_per_bank,
            valid=[1, 2, 3],
            onefile=True,
        )
        par_slope = self._get_param(
            "Slope",
            shipping_data["parslope0"],
            iq_str,
            num_stations_per_bank,
            positive=True,
            ext="slp",
        )
        par_reed = self._get_param(
            "Reed",
            shipping_data["parreed0"],
            iq_str,
            num_stations_per_bank,
            positive=True,
            ext="rdd",
        )

        mu_slope, mu_reed = [], []
        for ps, pr in zip(par_slope, par_reed):
            mus = ps.copy()
            mus[mus > 0] = 1.0 / mus[mus > 0]  # 1/slope for non-zero values
            mu_slope.append(mus)
            mu_reed.append(8.5e-4 * pr**0.8)  # empirical damping coefficient

        return SingleLevelParameters.from_column_arrays(
            {
                "id": level_i,
                "ship_velocity": ship_velocity,
                "num_ship": num_ship,
                "num_waves_per_ship": num_waves_per_ship,
                "ship_draught": ship_draught,
                "ship_type": ship_type,
                "par_slope": par_slope,
                "par_reed": par_reed,
                "mu_slope": mu_slope,
                "mu_reed": mu_reed,
            },
            SingleParameters,
        )

    def compute_erosion_per_level(
        self,
        level_i: int,
        bank_data: BankData,
        simulation_data: ErosionSimulationData,
        fairway_data: FairwayData,
        single_parameters: SingleLevelParameters,
        erosion_inputs: ErosionInputs,
        km_bin: Tuple[float, float, float],
        num_km: int,
    ) -> Tuple[SingleDischargeLevel, np.ndarray]:
        """Compute the bank erosion for a given level."""
        num_levels = self.river_data.num_discharge_levels
        dvol_bank = np.zeros((num_km, 2))
        hfw_max_level = 0
        par_list = []
        for ind, bank_i in enumerate(bank_data):

            single_calculation = SingleCalculation()
            # bank_i = 0: left bank, bank_i = 1: right bank
            # calculate velocity along banks ...
            single_calculation.bank_velocity = simulation_data.calculate_bank_velocity(
                bank_i, self.river_data.vel_dx
            )

            # get fairway face indices
            fairway_face_indices = bank_i.fairway_face_indices
            data = simulation_data.get_fairway_data(fairway_face_indices)
            single_calculation.water_level = data["water_level"]
            single_calculation.chezy = data["chezy"]
            single_calculation.water_depth = data["water_depth"]

            # get water depth along the fair-way
            hfw_max_level = max(hfw_max_level, single_calculation.water_depth.max())

            # last discharge level
            if level_i == num_levels - 1:
                erosion_distance_eq, erosion_volume_eq = self.erosion_calculator.comp_erosion_eq(
                    bank_i.height,
                    bank_i.segment_length,
                    fairway_data.fairway_initial_water_levels[ind],
                    single_parameters.get_bank(ind),
                    bank_i.fairway_distances,
                    single_calculation.water_depth,
                    erosion_inputs.get_bank(ind),
                )
                single_calculation.erosion_distance_eq = erosion_distance_eq
                single_calculation.erosion_volume_eq = erosion_volume_eq

            single_calculation = self.erosion_calculator.compute_bank_erosion_dynamics(
                single_calculation,
                bank_i.height,
                bank_i.segment_length,
                bank_i.fairway_distances,
                fairway_data.fairway_initial_water_levels[ind],
                single_parameters.get_bank(ind),
                self.river_data.erosion_time * self.p_discharge[level_i],
                erosion_inputs.get_bank(ind),
            )

            # accumulate eroded volumes per km
            volume_per_discharge = get_km_eroded_volume(
                bank_i.bank_chainage_midpoints, single_calculation.erosion_volume_tot, km_bin
            )
            single_calculation.volume_per_discharge = volume_per_discharge
            par_list.append(single_calculation)

            dvol_bank[:, ind] += volume_per_discharge

            if self.river_data.debug:
                self._debug_output(
                    level_i,
                    ind,
                    bank_data,
                    fairway_data,
                    erosion_inputs,
                    single_parameters,
                    num_levels,
                    single_calculation,
                )

        level_calculation = SingleDischargeLevel(
            left=par_list[0], right=par_list[1], hfw_max=hfw_max_level
        )

        return level_calculation, dvol_bank

    def _debug_output(
        self,
        level_i,
        ind,
        bank_data: BankData,
        fairway_data: FairwayData,
        erosion_inputs: ErosionInputs,
        single_level_parameters: SingleLevelParameters,
        num_levels: int,
        single_calculation: SingleCalculation,
    ):
        if level_i == num_levels - 1:
            # EQ debug
            self.debugger.last_discharge_level(
                ind,
                bank_data.get_bank(ind),
                fairway_data,
                erosion_inputs.get_bank(ind),
                single_level_parameters.get_bank(ind),
                single_calculation,
            )
        # Q-specific debug
        self.debugger.middle_levels(
            ind,
            level_i,
            bank_data.get_bank(ind),
            fairway_data,
            erosion_inputs.get_bank(ind),
            single_level_parameters.get_bank(ind),
            single_calculation,
        )

    def run(self) -> None:
        """Run the bank erosion analysis for a specified configuration."""
        timed_logger("-- start analysis --")
        log_text(
            "header_bankerosion",
            data={
                "version": __version__,
                "location": "https://github.com/Deltares/D-FAST_Bank_Erosion",
            },
        )
        log_text("-")
        config_file = self.config_file

        log_text("derive_topology")

        mesh_data = self.simulation_data.compute_mesh_topology()
        river_axis = self._process_river_axis_by_center_line()

        # map to the output interval
        km_bin = (
            river_axis.data["stations"].min(),
            river_axis.data["stations"].max(),
            self.river_data.output_intervals,
        )
        km_mid = get_km_bins(km_bin, station_type="mid")  # get mid-points

        fairway_data = self._get_fairway_data(river_axis, mesh_data)

        # map bank lines to mesh cells
        log_text("intersect_bank_mesh")
        bank_data = self.bl_processor.intersect_with_mesh(mesh_data)
        # map the bank data to the fairway data (the bank_data and fairway_data will be updated inside the `_map_bank_to_fairway` function)
        self.calculate_fairway_bank_line_distance(
            bank_data, fairway_data, self.simulation_data
        )

        erosion_inputs = self._prepare_initial_conditions(
            config_file, bank_data.num_stations_per_bank, fairway_data
        )

        # initialize arrays for erosion loop over all discharges
        water_level_data, erosion_results = self._process_discharge_levels(
            km_mid,
            km_bin,
            config_file,
            erosion_inputs,
            bank_data,
            fairway_data,
        )

        bankline_new_list, bankline_eq_list, xy_line_eq_list = (
            self._postprocess_erosion_results(
                km_bin,
                km_mid,
                bank_data,
                erosion_results,
            )
        )

        self._write_bankline_shapefiles(
            bankline_new_list, bankline_eq_list, config_file
        )
        self._write_volume_outputs(erosion_results, km_mid)

        # create various plots
        self._generate_plots(
            river_axis.data["stations"],
            self.simulation_data,
            xy_line_eq_list,
            km_mid,
            self.river_data.output_intervals,
            erosion_inputs,
            water_level_data,
            mesh_data,
            bank_data,
            erosion_results,
        )
        log_text("end_bankerosion")
        timed_logger("-- end analysis --")

    def _write_bankline_shapefiles(
        self, bankline_new_list, bankline_eq_list, config_file: ConfigFile
    ):
        bankline_new_series = GeoSeries(bankline_new_list, crs=config_file.crs)
        bank_lines_new = GeoDataFrame(geometry=bankline_new_series)
        bank_name = self.config_file.get_str("General", "BankFile", "bankfile")

        bank_file = self.river_data.output_dir / f"{bank_name}_new.shp"
        log_text("save_banklines", data={"file": str(bank_file)})
        bank_lines_new.to_file(bank_file)

        bankline_eq_series = GeoSeries(bankline_eq_list, crs=config_file.crs)
        banklines_eq = GeoDataFrame(geometry=bankline_eq_series)

        bank_file = self.river_data.output_dir / f"{bank_name}_eq.shp"
        log_text("save_banklines", data={"file": str(bank_file)})
        banklines_eq.to_file(bank_file)

    def _write_volume_outputs(self, erosion_results: ErosionResults, km_mid):
        erosion_vol_file = self.config_file.get_str(
            "Erosion", "EroVol", default="erovol.evo"
        )
        log_text("save_tot_erovol", data={"file": erosion_vol_file})
        write_km_eroded_volumes(
            km_mid,
            erosion_results.total_eroded_vol_per_km,
            str(self.river_data.output_dir / erosion_vol_file),
        )

        # write eroded volumes per km (equilibrium)
        erosion_vol_file = self.config_file.get_str(
            "Erosion", "EroVolEqui", default="erovol_eq.evo"
        )
        log_text("save_eq_erovol", data={"file": erosion_vol_file})
        write_km_eroded_volumes(
            km_mid,
            erosion_results.eq_eroded_vol_per_km,
            str(self.river_data.output_dir / erosion_vol_file),
        )

    def _generate_plots(
        self,
        river_axis_km,
        simulation_data: ErosionSimulationData,
        xy_line_eq_list,
        km_mid,
        km_step,
        erosion_inputs: ErosionInputs,
        water_level_data: WaterLevelData,
        mesh_data: MeshData,
        bank_data: BankData,
        erosion_results: ErosionResults,
    ):
        # create various plots
        if self.river_data.plot_flags["plot_data"]:
            log_text("=")
            log_text("create_figures")
            fig_i = 0
            bbox = self.river_data.get_bbox(self.river_center_line_arr)

            if self.river_data.plot_flags["save_plot_zoomed"]:
                bank_coords_mid = []
                for ib in range(bank_data.n_bank_lines):
                    bank_coords_mid.append(
                        (
                            bank_data.bank_line_coords[ib][:-1, :]
                            + bank_data.bank_line_coords[ib][1:, :]
                        )
                        / 2
                    )
                km_zoom, xy_zoom = get_zoom_extends(
                    river_axis_km.min(),
                    river_axis_km.max(),
                    self.river_data.plot_flags["zoom_km_step"],
                    bank_coords_mid,
                    bank_data.bank_chainage_midpoints,
                )

            fig, ax = df_plt.plot1_waterdepth_and_banklines(
                bbox,
                self.river_center_line_arr,
                bank_data.bank_lines,
                simulation_data.face_node,
                simulation_data.n_nodes,
                simulation_data.x_node,
                simulation_data.y_node,
                simulation_data.water_depth_face,
                1.1 * water_level_data.hfw_max,
                X_AXIS_TITLE,
                Y_AXIS_TITLE,
                "water depth and initial bank lines",
                "water depth [m]",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    f"{self.river_data.plot_flags['fig_dir']}{os.sep}{fig_i}_banklines"
                )

                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_xy_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        xy_zoom,
                    )

                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            fig, ax = df_plt.plot2_eroded_distance_and_equilibrium(
                bbox,
                self.river_center_line_arr,
                bank_data.bank_line_coords,
                erosion_results.total_erosion_dist,
                bank_data.is_right_bank,
                erosion_results.avg_erosion_rate,
                xy_line_eq_list,
                mesh_data.x_edge_coords,
                mesh_data.y_edge_coords,
                X_AXIS_TITLE,
                Y_AXIS_TITLE,
                "eroded distance and equilibrium bank location",
                f"eroded during {erosion_results.erosion_time} year",
                "eroded distance [m]",
                "equilibrium location",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = f"{self.river_data.plot_flags['fig_dir']}{os.sep}{fig_i}_erosion_sensitivity"

                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_xy_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        xy_zoom,
                    )

                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            fig, ax = df_plt.plot3_eroded_volume(
                km_mid,
                km_step,
                "river chainage [km]",
                water_level_data.vol_per_discharge,
                "eroded volume [m^3]",
                f"eroded volume per {km_step} chainage km ({erosion_results.erosion_time} years)",
                "Q{iq}",
                "Bank {ib}",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = f"{self.river_data.plot_flags['fig_dir']}{os.sep}{fig_i}_eroded_volume"

                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_x_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        km_zoom,
                    )

                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            fig, ax = df_plt.plot3_eroded_volume_subdivided_1(
                km_mid,
                km_step,
                "river chainage [km]",
                water_level_data.vol_per_discharge,
                "eroded volume [m^3]",
                f"eroded volume per {km_step} chainage km ({erosion_results.erosion_time} years)",
                "Q{iq}",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    self.river_data.plot_flags["fig_dir"]
                    + os.sep
                    + str(fig_i)
                    + "_eroded_volume_per_discharge"
                )
                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_x_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        km_zoom,
                    )
                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            fig, ax = df_plt.plot3_eroded_volume_subdivided_2(
                km_mid,
                km_step,
                "river chainage [km]",
                water_level_data.vol_per_discharge,
                "eroded volume [m^3]",
                f"eroded volume per {km_step} chainage km ({erosion_results.erosion_time} years)",
                "Bank {ib}",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    self.river_data.plot_flags["fig_dir"]
                    + os.sep
                    + str(fig_i)
                    + "_eroded_volume_per_bank"
                )
                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_x_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        km_zoom,
                    )
                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            fig, ax = df_plt.plot4_eroded_volume_eq(
                km_mid,
                km_step,
                "river chainage [km]",
                erosion_results.eq_eroded_vol_per_km,
                "eroded volume [m^3]",
                f"eroded volume per {km_step} chainage km (equilibrium)",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    self.river_data.plot_flags["fig_dir"]
                    + os.sep
                    + str(fig_i)
                    + "_eroded_volume_eq"
                )
                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_x_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        km_zoom,
                    )
                fig_path = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_path)

            figlist, axlist = df_plt.plot5series_waterlevels_per_bank(
                bank_data.bank_chainage_midpoints,
                "river chainage [km]",
                water_level_data.water_level,
                water_level_data.ship_wave_max,
                water_level_data.ship_wave_min,
                "water level at Q{iq}",
                "average water level",
                "wave influenced range",
                bank_data.height,
                "level of bank",
                erosion_inputs.bank_protection_level,
                "bank protection level",
                "elevation",
                "(water)levels along bank line {ib}",
                "[m NAP]",
            )
            if self.river_data.plot_flags["save_plot"]:
                for ib, fig in enumerate(figlist):
                    fig_i = fig_i + 1
                    fig_base = f"{self.river_data.plot_flags['fig_dir']}/{fig_i}_levels_bank_{ib + 1}"

                    if self.river_data.plot_flags["save_plot_zoomed"]:
                        df_plt.zoom_x_and_save(
                            fig,
                            axlist[ib],
                            fig_base,
                            self.river_data.plot_flags["plot_ext"],
                            km_zoom,
                        )
                    fig_file = f"{fig_base}{self.river_data.plot_flags['plot_ext']}"
                    df_plt.savefig(fig, fig_file)

            figlist, axlist = df_plt.plot6series_velocity_per_bank(
                bank_data.bank_chainage_midpoints,
                "river chainage [km]",
                water_level_data.velocity,
                "velocity at Q{iq}",
                erosion_inputs.tauc,
                water_level_data.chezy[0],
                "critical velocity",
                "velocity",
                "velocity along bank line {ib}",
                "[m/s]",
            )
            if self.river_data.plot_flags["save_plot"]:
                for ib, fig in enumerate(figlist):
                    fig_i = fig_i + 1
                    fig_base = f"{self.river_data.plot_flags['fig_dir']}{os.sep}{fig_i}_velocity_bank_{ib + 1}"

                    if self.river_data.plot_flags["save_plot_zoomed"]:
                        df_plt.zoom_x_and_save(
                            fig,
                            axlist[ib],
                            fig_base,
                            self.river_data.plot_flags["plot_ext"],
                            km_zoom,
                        )

                    fig_file = fig_base + self.river_data.plot_flags["plot_ext"]
                    df_plt.savefig(fig, fig_file)

            fig, ax = df_plt.plot7_banktype(
                bbox,
                self.river_center_line_arr,
                bank_data.bank_line_coords,
                erosion_inputs.bank_type,
                erosion_inputs.taucls_str,
                X_AXIS_TITLE,
                Y_AXIS_TITLE,
                "bank type",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    self.river_data.plot_flags["fig_dir"]
                    + os.sep
                    + str(fig_i)
                    + "_banktype"
                )
                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_xy_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        xy_zoom,
                    )
                fig_file = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_file)

            fig, ax = df_plt.plot8_eroded_distance(
                bank_data.bank_chainage_midpoints,
                "river chainage [km]",
                erosion_results.total_erosion_dist,
                "Bank {ib}",
                erosion_results.eq_erosion_dist,
                "Bank {ib} (eq)",
                "eroded distance",
                "[m]",
            )
            if self.river_data.plot_flags["save_plot"]:
                fig_i = fig_i + 1
                fig_base = (
                    self.river_data.plot_flags["fig_dir"]
                    + os.sep
                    + str(fig_i)
                    + "_erodis"
                )
                if self.river_data.plot_flags["save_plot_zoomed"]:
                    df_plt.zoom_x_and_save(
                        fig,
                        ax,
                        fig_base,
                        self.river_data.plot_flags["plot_ext"],
                        km_zoom,
                    )
                fig_file = fig_base + self.river_data.plot_flags["plot_ext"]
                df_plt.savefig(fig, fig_file)

            if self.river_data.plot_flags["close_plot"]:
                plt.close("all")
            else:
                plt.show(block=not self.gui)

`config_file: ConfigFile` `property` #

Configuration file object.

`init(config_file: ConfigFile, gui: bool = False)` #

Initialize the Erosion class.

Source code in src/dfastbe/bank_erosion/bank_erosion.py

def __init__(self, config_file: ConfigFile, gui: bool = False):
    """Initialize the Erosion class."""
    self.root_dir = config_file.root_dir
    self._config_file = config_file
    self.gui = gui

    self.river_data = ErosionRiverData(config_file)
    self.river_center_line_arr = self.river_data.river_center_line.as_array()
    self.simulation_data = self.river_data.simulation_data()
    self.sim_files, self.p_discharge = self.river_data.get_erosion_sim_data(
        self.river_data.num_discharge_levels
    )
    self.bl_processor = BankLinesProcessor(self.river_data)
    self.debugger = Debugger(config_file.crs, self.river_data.output_dir)
    self.erosion_calculator = ErosionCalculator()

`calculate_fairway_bank_line_distance(bank_data: BankData, fairway_data: FairwayData, simulation_data: ErosionSimulationData)` #

Map bank data to fairway data.

Parameters:

Name	Type	Default
`bank_data`	`BankData`	required
`fairway_data`	`FairwayData`	required
`simulation_data`	`ErosionSimulationData`	required

Returns:

Name	Type	Description
`FairwayData`		The method updates the following attributes in the `bank_data` instance - fairway_face_indices - fairway_distances
`BankData`		the following attributes in the `fairway_data` instance - fairway_initial_water_levels

Source code in src/dfastbe/bank_erosion/bank_erosion.py

def calculate_fairway_bank_line_distance(
    self,
    bank_data: BankData,
    fairway_data: FairwayData,
    simulation_data: ErosionSimulationData,
):
    """Map bank data to fairway data.

    Args:
        bank_data (BankData):
        fairway_data (FairwayData):
        simulation_data (ErosionSimulationData):

    Returns:
        FairwayData:
            The method updates the following attributes in the `bank_data` instance
                - fairway_face_indices
                - fairway_distances
        BankData:
            the following attributes in the `fairway_data` instance
                - fairway_initial_water_levels
    """
    # distance fairway-bankline (bank-fairway)
    log_text("bank_distance_fairway")

    num_fairway_face_ind = len(fairway_data.fairway_face_indices)

    for bank_i, single_bank in enumerate(bank_data):
        bank_coords = single_bank.bank_line_coords
        coords_mid = (bank_coords[:-1] + bank_coords[1:]) / 2
        bank_fairway_dist = np.zeros(len(coords_mid))
        bp_fw_face_idx = np.zeros(len(coords_mid), dtype=int)

        for ind, coord_i in enumerate(coords_mid):
            # find closest fairway support node
            closest_ind = np.argmin(
                ((coord_i - fairway_data.intersection_coords) ** 2).sum(axis=1)
            )
            fairway_coord = fairway_data.intersection_coords[closest_ind]
            fairway_bank_distance = ((coord_i - fairway_coord) ** 2).sum() ** 0.5
            # If fairway support node is also the closest projected fairway point, then it likely
            # that that point is one of the original support points (a corner) of the fairway path
            # and located inside a grid cell. The segments before and after that point will then
            # both be located inside that same grid cell, so let's pick the segment before the point.
            # If the point happens to coincide with a grid edge and the two segments are located
            # in different grid cells, then we could either simply choose one or add complexity to
            # average the values of the two grid cells. Let's go for the simplest approach ...
            iseg = max(closest_ind - 1, 0)
            if closest_ind > 0:
                alpha = calculate_alpha(
                    fairway_data.intersection_coords,
                    closest_ind,
                    closest_ind - 1,
                    coord_i,
                )
                if 0 < alpha < 1:
                    fwp1 = fairway_data.intersection_coords[
                        closest_ind - 1
                    ] + alpha * (
                        fairway_data.intersection_coords[closest_ind]
                        - fairway_data.intersection_coords[closest_ind - 1]
                    )
                    d1 = ((coord_i - fwp1) ** 2).sum() ** 0.5
                    if d1 < fairway_bank_distance:
                        fairway_bank_distance = d1
                        # projected point located on segment before, which corresponds to initial choice: iseg = ifw - 1
            if closest_ind < num_fairway_face_ind:
                alpha = calculate_alpha(
                    fairway_data.intersection_coords,
                    closest_ind + 1,
                    closest_ind,
                    coord_i,
                )
                if 0 < alpha < 1:
                    fwp1 = fairway_data.intersection_coords[closest_ind] + alpha * (
                        fairway_data.intersection_coords[closest_ind + 1]
                        - fairway_data.intersection_coords[closest_ind]
                    )
                    d1 = ((coord_i - fwp1) ** 2).sum() ** 0.5
                    if d1 < fairway_bank_distance:
                        fairway_bank_distance = d1
                        iseg = closest_ind

            bp_fw_face_idx[ind] = fairway_data.fairway_face_indices[iseg]
            bank_fairway_dist[ind] = fairway_bank_distance

        if self.river_data.debug:
            line_geom = LineGeometry(coords_mid, crs=self.config_file.crs)
            line_geom.to_file(
                file_name=f"{self.river_data.output_dir}/bank_{bank_i + 1}_chainage_and_fairway_face_idx.shp",
                data={
                    "chainage": single_bank.bank_chainage_midpoints,
                    "iface_fw": bp_fw_face_idx[bank_i],
                },
            )

        single_bank.fairway_face_indices = bp_fw_face_idx
        single_bank.fairway_distances = bank_fairway_dist

    # water level at fairway
    water_level_fairway_ref = []
    for single_bank in bank_data:
        ii = single_bank.fairway_face_indices
        water_level_fairway_ref.append(simulation_data.water_level_face[ii])
    fairway_data.fairway_initial_water_levels = water_level_fairway_ref

`compute_erosion_per_level(level_i: int, bank_data: BankData, simulation_data: ErosionSimulationData, fairway_data: FairwayData, single_parameters: SingleLevelParameters, erosion_inputs: ErosionInputs, km_bin: Tuple[float, float, float], num_km: int) -> Tuple[SingleDischargeLevel, np.ndarray]` #

Compute the bank erosion for a given level.

Source code in src/dfastbe/bank_erosion/bank_erosion.py

def compute_erosion_per_level(
    self,
    level_i: int,
    bank_data: BankData,
    simulation_data: ErosionSimulationData,
    fairway_data: FairwayData,
    single_parameters: SingleLevelParameters,
    erosion_inputs: ErosionInputs,
    km_bin: Tuple[float, float, float],
    num_km: int,
) -> Tuple[SingleDischargeLevel, np.ndarray]:
    """Compute the bank erosion for a given level."""
    num_levels = self.river_data.num_discharge_levels
    dvol_bank = np.zeros((num_km, 2))
    hfw_max_level = 0
    par_list = []
    for ind, bank_i in enumerate(bank_data):

        single_calculation = SingleCalculation()
        # bank_i = 0: left bank, bank_i = 1: right bank
        # calculate velocity along banks ...
        single_calculation.bank_velocity = simulation_data.calculate_bank_velocity(
            bank_i, self.river_data.vel_dx
        )

        # get fairway face indices
        fairway_face_indices = bank_i.fairway_face_indices
        data = simulation_data.get_fairway_data(fairway_face_indices)
        single_calculation.water_level = data["water_level"]
        single_calculation.chezy = data["chezy"]
        single_calculation.water_depth = data["water_depth"]

        # get water depth along the fair-way
        hfw_max_level = max(hfw_max_level, single_calculation.water_depth.max())

        # last discharge level
        if level_i == num_levels - 1:
            erosion_distance_eq, erosion_volume_eq = self.erosion_calculator.comp_erosion_eq(
                bank_i.height,
                bank_i.segment_length,
                fairway_data.fairway_initial_water_levels[ind],
                single_parameters.get_bank(ind),
                bank_i.fairway_distances,
                single_calculation.water_depth,
                erosion_inputs.get_bank(ind),
            )
            single_calculation.erosion_distance_eq = erosion_distance_eq
            single_calculation.erosion_volume_eq = erosion_volume_eq

        single_calculation = self.erosion_calculator.compute_bank_erosion_dynamics(
            single_calculation,
            bank_i.height,
            bank_i.segment_length,
            bank_i.fairway_distances,
            fairway_data.fairway_initial_water_levels[ind],
            single_parameters.get_bank(ind),
            self.river_data.erosion_time * self.p_discharge[level_i],
            erosion_inputs.get_bank(ind),
        )

        # accumulate eroded volumes per km
        volume_per_discharge = get_km_eroded_volume(
            bank_i.bank_chainage_midpoints, single_calculation.erosion_volume_tot, km_bin
        )
        single_calculation.volume_per_discharge = volume_per_discharge
        par_list.append(single_calculation)

        dvol_bank[:, ind] += volume_per_discharge

        if self.river_data.debug:
            self._debug_output(
                level_i,
                ind,
                bank_data,
                fairway_data,
                erosion_inputs,
                single_parameters,
                num_levels,
                single_calculation,
            )

    level_calculation = SingleDischargeLevel(
        left=par_list[0], right=par_list[1], hfw_max=hfw_max_level
    )

    return level_calculation, dvol_bank

`get_ship_parameters(num_stations_per_bank: List[int]) -> Dict[str, List[np.ndarray]]` #

Get ship parameters from the configuration file.

Source code in src/dfastbe/bank_erosion/bank_erosion.py

def get_ship_parameters(
    self, num_stations_per_bank: List[int]
) -> Dict[str, List[np.ndarray]]:
    """Get ship parameters from the configuration file."""
    ship_relative_velocity = self.config_file.get_parameter(
        "Erosion", "VShip", num_stations_per_bank, positive=True, onefile=True
    )
    num_ships_year = self.config_file.get_parameter(
        "Erosion", "NShip", num_stations_per_bank, positive=True, onefile=True
    )
    num_waves_p_ship = self.config_file.get_parameter(
        "Erosion",
        "NWave",
        num_stations_per_bank,
        default=5,
        positive=True,
        onefile=True,
    )
    ship_draught = self.config_file.get_parameter(
        "Erosion", "Draught", num_stations_per_bank, positive=True, onefile=True
    )
    ship_type = self.config_file.get_parameter(
        "Erosion", "ShipType", num_stations_per_bank, valid=[1, 2, 3], onefile=True
    )
    parslope0 = self.config_file.get_parameter(
        "Erosion",
        "Slope",
        num_stations_per_bank,
        default=20,
        positive=True,
        ext="slp",
    )
    reed_wave_damping_coeff = self.config_file.get_parameter(
        "Erosion",
        "Reed",
        num_stations_per_bank,
        default=0,
        positive=True,
        ext="rdd",
    )

    ship_data = {
        "vship0": ship_relative_velocity,
        "Nship0": num_ships_year,
        "nwave0": num_waves_p_ship,
        "Tship0": ship_draught,
        "ship0": ship_type,
        "parslope0": parslope0,
        "parreed0": reed_wave_damping_coeff,
    }
    return ship_data

`run() -> None` #

Run the bank erosion analysis for a specified configuration.

Source code in src/dfastbe/bank_erosion/bank_erosion.py

def run(self) -> None:
    """Run the bank erosion analysis for a specified configuration."""
    timed_logger("-- start analysis --")
    log_text(
        "header_bankerosion",
        data={
            "version": __version__,
            "location": "https://github.com/Deltares/D-FAST_Bank_Erosion",
        },
    )
    log_text("-")
    config_file = self.config_file

    log_text("derive_topology")

    mesh_data = self.simulation_data.compute_mesh_topology()
    river_axis = self._process_river_axis_by_center_line()

    # map to the output interval
    km_bin = (
        river_axis.data["stations"].min(),
        river_axis.data["stations"].max(),
        self.river_data.output_intervals,
    )
    km_mid = get_km_bins(km_bin, station_type="mid")  # get mid-points

    fairway_data = self._get_fairway_data(river_axis, mesh_data)

    # map bank lines to mesh cells
    log_text("intersect_bank_mesh")
    bank_data = self.bl_processor.intersect_with_mesh(mesh_data)
    # map the bank data to the fairway data (the bank_data and fairway_data will be updated inside the `_map_bank_to_fairway` function)
    self.calculate_fairway_bank_line_distance(
        bank_data, fairway_data, self.simulation_data
    )

    erosion_inputs = self._prepare_initial_conditions(
        config_file, bank_data.num_stations_per_bank, fairway_data
    )

    # initialize arrays for erosion loop over all discharges
    water_level_data, erosion_results = self._process_discharge_levels(
        km_mid,
        km_bin,
        config_file,
        erosion_inputs,
        bank_data,
        fairway_data,
    )

    bankline_new_list, bankline_eq_list, xy_line_eq_list = (
        self._postprocess_erosion_results(
            km_bin,
            km_mid,
            bank_data,
            erosion_results,
        )
    )

    self._write_bankline_shapefiles(
        bankline_new_list, bankline_eq_list, config_file
    )
    self._write_volume_outputs(erosion_results, km_mid)

    # create various plots
    self._generate_plots(
        river_axis.data["stations"],
        self.simulation_data,
        xy_line_eq_list,
        km_mid,
        self.river_data.output_intervals,
        erosion_inputs,
        water_level_data,
        mesh_data,
        bank_data,
        erosion_results,
    )
    log_text("end_bankerosion")
    timed_logger("-- end analysis --")

Mesh Processing #

`dfastbe.bank_erosion.mesh_processor` #

module for processing mesh-related operations.

`enlarge(old_array: np.ndarray, new_shape: Tuple)` #

Copy the values of the old array to a new, larger array of specified shape.

Arguments #

old_array : numpy.ndarray Array containing the values. new_shape : Tuple New shape of the array.

Returns #

new_array : numpy.ndarray Array of shape "new_shape" with the 'first entries filled by the same values as contained in "old_array". The data type of the new array is equal to that of the old array.

Source code in src/dfastbe/bank_erosion/mesh_processor.py

def enlarge(
    old_array: np.ndarray,
    new_shape: Tuple
):
    """
    Copy the values of the old array to a new, larger array of specified shape.

    Arguments
    ---------
    old_array : numpy.ndarray
        Array containing the values.
    new_shape : Tuple
        New shape of the array.

    Returns
    -------
    new_array : numpy.ndarray
        Array of shape "new_shape" with the 'first entries filled by the same
        values as contained in "old_array". The data type of the new array is
        equal to that of the old array.
    """
    old_shape = old_array.shape
    print("old: ", old_shape)
    print("new: ", new_shape)
    new_array = np.zeros(new_shape, dtype=old_array.dtype)
    if len(new_shape)==1:
        new_array[:old_shape[0]] = old_array
    elif len(new_shape)==2:
        new_array[:old_shape[0], :old_shape[1]] = old_array
    return new_array

`get_slices_ab(X0: np.ndarray, Y0: np.ndarray, X1: np.ndarray, Y1: np.ndarray, xi0: float, yi0: float, xi1: float, yi1: float, bmin: float, bmax1: bool = True) -> Tuple[np.ndarray, np.ndarray, np.ndarray]` #

Get the relative locations a and b at which a segment intersects/slices a number of edges.

Arguments #

X0 : np.ndarray Array containing the x-coordinates of the start point of each edge. X1 : np.ndarray Array containing the x-coordinates of the end point of each edge. Y0 : np.ndarray Array containing the y-coordinates of the start point of each edge. Y1 : np.ndarray Array containing the y-coordinates of the end point of each edge. xi0 : float x-coordinate of start point of the segment. xi1 : float x-coordinate of end point of the segment. yi0 : float y-coordinate of start point of the segment. yi1 : float y-coordinate of end point of the segment. bmin : float Minimum relative distance from bpj1 at which slice should occur. bmax1 : bool Flag indicating whether the the relative distance along the segment bpj1-bpj should be limited to 1.

Returns #

a : np.ndarray Array containing relative distance along each edge. b : np.ndarray Array containing relative distance along the segment for each edge. slices : np.ndarray Array containing a flag indicating whether the edge is sliced at a valid location.

Source code in src/dfastbe/bank_erosion/mesh_processor.py

def get_slices_ab(
    X0: np.ndarray,
    Y0: np.ndarray,
    X1: np.ndarray,
    Y1: np.ndarray,
    xi0: float,
    yi0: float,
    xi1: float,
    yi1: float,
    bmin: float,
    bmax1: bool = True,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Get the relative locations a and b at which a segment intersects/slices a number of edges.

    Arguments
    ---------
    X0 : np.ndarray
        Array containing the x-coordinates of the start point of each edge.
    X1 : np.ndarray
        Array containing the x-coordinates of the end point of each edge.
    Y0 : np.ndarray
        Array containing the y-coordinates of the start point of each edge.
    Y1 : np.ndarray
        Array containing the y-coordinates of the end point of each edge.
    xi0 : float
        x-coordinate of start point of the segment.
    xi1 : float
        x-coordinate of end point of the segment.
    yi0 : float
        y-coordinate of start point of the segment.
    yi1 : float
        y-coordinate of end point of the segment.
    bmin : float
        Minimum relative distance from bpj1 at which slice should occur.
    bmax1 : bool
        Flag indicating whether the the relative distance along the segment bpj1-bpj should be limited to 1.

    Returns
    -------
    a : np.ndarray
        Array containing relative distance along each edge.
    b : np.ndarray
        Array containing relative distance along the segment for each edge.
    slices : np.ndarray
        Array containing a flag indicating whether the edge is sliced at a valid location.
    """
    dX = X1 - X0
    dY = Y1 - Y0
    dxi = xi1 - xi0
    dyi = yi1 - yi0
    det = dX * dyi - dY * dxi
    det[det == 0] = 1e-10
    a = (dyi * (xi0 - X0) - dxi * (yi0 - Y0)) / det  # along mesh edge
    b = (dY * (xi0 - X0) - dX * (yi0 - Y0)) / det  # along bank line
    if bmax1:
        slices = np.nonzero((b > bmin) & (b <= 1) & (a >= 0) & (a <= 1))[0]
    else:
        slices = np.nonzero((b > bmin) & (a >= 0) & (a <= 1))[0]
    a = a[slices]
    b = b[slices]
    return a, b, slices

`intersect_line_mesh(bp: np.ndarray, mesh_data: MeshData, d_thresh: float = 0.001) -> Tuple[np.ndarray, np.ndarray]` #

Intersects a line with an unstructured mesh and returns the intersection coordinates and mesh face indices.

This function determines where a given line (e.g., a bank line) intersects the faces of an unstructured mesh. It calculates the intersection points and identifies the mesh faces corresponding to each segment of the line.

Parameters:

Name	Type	Description	Default
`bp`	`ndarray`	A 2D array of shape (N, 2) containing the x, y coordinates of the line to be intersected with the mesh.	required
`mesh_data`	`MeshData`	An instance of the `MeshData` class containing mesh-related data, such as face coordinates, edge coordinates, and connectivity information.	required
`d_thresh`	`float`	A distance threshold for filtering out very small segments. Segments shorter than this threshold will be removed. Defaults to 0.001.	`0.001`

Returns:

Type	Description
`Tuple[ndarray, ndarray]`	Tuple[np.ndarray, np.ndarray]: A tuple containing: - `crds` (np.ndarray): A 2D array of shape (M, 2) containing the x, y coordinates of the intersection points. - `idx` (np.ndarray): A 1D array of shape (M-1,) containing the indices of the mesh faces corresponding to each segment of the intersected line.

Raises:

Type	Description
`Exception`	If the line starts outside the mesh and cannot be associated with any mesh face, or if the line crosses ambiguous regions (e.g., edges shared by multiple faces).

Notes

The function uses Shapely geometry operations to determine whether points are inside polygons or on edges.
The function handles cases where the line starts outside the mesh, crosses multiple edges, or ends on a node.
Tiny segments shorter than d_thresh are removed from the output.

Source code in src/dfastbe/bank_erosion/mesh_processor.py

def intersect_line_mesh(
    bp: np.ndarray,
    mesh_data: MeshData,
    d_thresh: float = 0.001,
) -> Tuple[np.ndarray, np.ndarray]:
    """Intersects a line with an unstructured mesh and returns the intersection coordinates and mesh face indices.

    This function determines where a given line (e.g., a bank line) intersects the faces of an unstructured mesh.
    It calculates the intersection points and identifies the mesh faces corresponding to each segment of the line.

    Args:
        bp (np.ndarray):
            A 2D array of shape (N, 2) containing the x, y coordinates of the line to be intersected with the mesh.
        mesh_data (MeshData):
            An instance of the `MeshData` class containing mesh-related data, such as face coordinates, edge coordinates,
            and connectivity information.
        d_thresh (float, optional):
            A distance threshold for filtering out very small segments. Segments shorter than this threshold will be removed.
            Defaults to 0.001.

    Returns:
        Tuple[np.ndarray, np.ndarray]:
            A tuple containing:
            - `crds` (np.ndarray): A 2D array of shape (M, 2) containing the x, y coordinates of the intersection points.
            - `idx` (np.ndarray): A 1D array of shape (M-1,) containing the indices of the mesh faces corresponding to
              each segment of the intersected line.

    Raises:
        Exception:
            If the line starts outside the mesh and cannot be associated with any mesh face, or if the line crosses
            ambiguous regions (e.g., edges shared by multiple faces).

    Notes:
        - The function uses Shapely geometry operations to determine whether points are inside polygons or on edges.
        - The function handles cases where the line starts outside the mesh, crosses multiple edges, or ends on a node.
        - Tiny segments shorter than `d_thresh` are removed from the output.
    """
    crds = np.zeros((len(bp), 2))
    idx = np.zeros(len(bp), dtype=np.int64)
    verbose = False
    ind = 0
    index: int
    vindex: np.ndarray
    for j, bpj in enumerate(bp):
        if verbose:
            print("Current location: {}, {}".format(bpj[0], bpj[1]))
        if j == 0:
            # first bp inside or outside?
            dx = mesh_data.x_face_coords - bpj[0]
            dy = mesh_data.y_face_coords - bpj[1]
            possible_cells = np.nonzero(
                ~(
                        (dx < 0).all(axis=1)
                        | (dx > 0).all(axis=1)
                        | (dy < 0).all(axis=1)
                        | (dy > 0).all(axis=1)
                )
            )[0]
            if len(possible_cells) == 0:
                # no cells found ... it must be outside
                index = -1
                if verbose:
                    print("starting outside mesh")
            else:
                # one or more possible cells, check whether it's really inside one of them
                # using np math might be faster, but since it's should only be for a few points let's use Shapely
                # a point on the edge of a polygon is not contained in the polygon.
                # a point on the edge of two polygons will thus be considered outside the mesh whereas it actually isn't.
                pnt = Point(bp[0])
                for k in possible_cells:
                    polygon_k = Polygon(
                        np.concatenate(
                            (
                                mesh_data.x_face_coords[
                                k : k + 1, : mesh_data.n_nodes[k]
                                ],
                                mesh_data.y_face_coords[
                                k : k + 1, : mesh_data.n_nodes[k]
                                ],
                            ),
                            axis=0,
                        ).T
                    )
                    if polygon_k.contains(pnt):
                        index = k
                        if verbose:
                            print("starting in {}".format(index))
                        break
                else:
                    on_edge: List[int] = []
                    for k in possible_cells:
                        nd = np.concatenate(
                            (
                                mesh_data.x_face_coords[
                                k : k + 1, : mesh_data.n_nodes[k]
                                ],
                                mesh_data.y_face_coords[
                                k : k + 1, : mesh_data.n_nodes[k]
                                ],
                            ),
                            axis=0,
                        ).T
                        line_k = LineString(np.concatenate(nd, nd[0:1], axis=0))
                        if line_k.contains(pnt):
                            on_edge.append(k)
                    if not on_edge:
                        index = -1
                        if verbose:
                            print("starting outside mesh")
                    else:
                        if len(on_edge) == 1:
                            index = on_edge[0]
                        else:
                            index = -2
                            vindex = on_edge
                        if verbose:
                            print("starting on edge of {}".format(on_edge))
                        raise Exception("determine direction!")
            crds[ind] = bpj
            if index == -2:
                idx[ind] = vindex[0]
            else:
                idx[ind] = index
            ind += 1
        else:
            # second or later point
            bpj1 = bp[j - 1]
            prev_b = 0
            while True:
                if index == -2:
                    b = np.zeros(0)
                    edges = np.zeros(0, dtype=np.int64)
                    nodes = np.zeros(0, dtype=np.int64)
                    index_src = np.zeros(0, dtype=np.int64)
                    for i in vindex:
                        b1, edges1, nodes1 = _get_slices(
                            i,
                            prev_b,
                            bpj,
                            bpj1,
                            mesh_data,
                        )
                        b = np.concatenate((b, b1), axis=0)
                        edges = np.concatenate((edges, edges1), axis=0)
                        nodes = np.concatenate((nodes, nodes1), axis=0)
                        index_src = np.concatenate((index_src, i + 0 * edges1), axis=0)
                    edges, id_edges = np.unique(edges, return_index=True)
                    b = b[id_edges]
                    nodes = nodes[id_edges]
                    index_src = index_src[id_edges]
                    if len(index_src) == 1:
                        index = index_src[0]
                        vindex = index_src[0:1]
                elif (bpj == bpj1).all():
                    # this is a segment of length 0, skip it since it takes us nowhere
                    break
                else:
                    b, edges, nodes = _get_slices(
                        index,
                        prev_b,
                        bpj,
                        bpj1,
                        mesh_data,
                    )

                if len(edges) == 0:
                    # rest of segment associated with same face
                    if verbose:
                        if prev_b > 0:
                            print(
                                "{}: -- no further slices along this segment --".format(
                                    j
                                )
                            )
                        else:
                            print("{}: -- no slices along this segment --".format(j))
                        if index >= 0:
                            pnt = Point(bpj)
                            polygon_k = Polygon(
                                np.concatenate(
                                    (
                                        mesh_data.x_face_coords[
                                        index : index + 1,
                                        : mesh_data.n_nodes[index],
                                        ],
                                        mesh_data.y_face_coords[
                                        index : index + 1,
                                        : mesh_data.n_nodes[index],
                                        ],
                                    ),
                                    axis=0,
                                ).T
                            )
                            if not polygon_k.contains(pnt):
                                raise Exception(
                                    "{}: ERROR: point actually not contained within {}!".format(
                                        j, index
                                    )
                                )
                    if ind == crds.shape[0]:
                        crds = enlarge(crds, (2 * ind, 2))
                        idx = enlarge(idx, (2 * ind,))
                    crds[ind] = bpj
                    idx[ind] = index
                    ind += 1
                    break
                else:
                    index0 = None
                    if len(edges) > 1:
                        # line segment crosses the edge list multiple times
                        # - moving out of a cell at a corner node
                        # - moving into and out of the mesh from outside
                        # select first crossing ...
                        bmin = b == np.amin(b)
                        b = b[bmin]
                        edges = edges[bmin]
                        nodes = nodes[bmin]

                    # slice location identified ...
                    node = nodes[0]
                    edge = edges[0]
                    faces = mesh_data.edge_face_connectivity[edge]
                    prev_b = b[0]

                    if node >= 0:
                        # if we slice at a node ...
                        if verbose:
                            print(
                                "{}: moving via node {} on edges {} at {}".format(
                                    j, node, edges, b[0]
                                )
                            )
                        # figure out where we will be heading afterwards ...
                        all_node_edges = np.nonzero(
                            (mesh_data.edge_node == node).any(axis=1)
                        )[0]
                        if b[0] < 1.0:
                            # segment passes through node and enter non-neighbouring cell ...
                            # direction of current segment from bpj1 to bpj
                            theta = math.atan2(bpj[1] - bpj1[1], bpj[0] - bpj1[0])
                        else:
                            if b[0] == 1.0 and j == len(bp) - 1:
                                # catch case of last segment
                                if verbose:
                                    print("{}: last point ends in a node".format(j))
                                if ind == crds.shape[0]:
                                    crds = enlarge(crds, (ind + 1, 2))
                                    idx = enlarge(idx, (ind + 1,))
                                crds[ind] = bpj
                                if index == -2:
                                    idx[ind] = vindex[0]
                                else:
                                    idx[ind] = index
                                ind += 1
                                break
                            else:
                                # this segment ends in the node, so check next segment ...
                                # direction of next segment from bpj to bp[j+1]
                                theta = math.atan2(
                                    bp[j + 1][1] - bpj[1], bp[j + 1][0] - bp[j][0]
                                )
                        if verbose:
                            print("{}: moving in direction theta = {}".format(j, theta))
                        twopi = 2 * math.pi
                        left_edge = -1
                        left_dtheta = twopi
                        right_edge = -1
                        right_dtheta = twopi
                        if verbose:
                            print(
                                "{}: the edges connected to node {} are {}".format(
                                    j, node, all_node_edges
                                )
                            )
                        for ie in all_node_edges:
                            if mesh_data.edge_node[ie, 0] == node:
                                theta_edge = math.atan2(
                                    mesh_data.y_edge_coords[ie, 1]
                                    - mesh_data.y_edge_coords[ie, 0],
                                    mesh_data.x_edge_coords[ie, 1]
                                    - mesh_data.x_edge_coords[ie, 0],
                                    )
                            else:
                                theta_edge = math.atan2(
                                    mesh_data.y_edge_coords[ie, 0]
                                    - mesh_data.y_edge_coords[ie, 1],
                                    mesh_data.x_edge_coords[ie, 0]
                                    - mesh_data.x_edge_coords[ie, 1],
                                    )
                            if verbose:
                                print(
                                    "{}: edge {} connects {}".format(
                                        j, ie, mesh_data.edge_node[ie, :]
                                    )
                                )
                                print(
                                    "{}: edge {} theta is {}".format(j, ie, theta_edge)
                                )
                            dtheta = theta_edge - theta
                            if dtheta > 0:
                                if dtheta < left_dtheta:
                                    left_edge = ie
                                    left_dtheta = dtheta
                                if twopi - dtheta < right_dtheta:
                                    right_edge = ie
                                    right_dtheta = twopi - dtheta
                            elif dtheta < 0:
                                dtheta = -dtheta
                                if twopi - dtheta < left_dtheta:
                                    left_edge = ie
                                    left_dtheta = twopi - dtheta
                                if dtheta < right_dtheta:
                                    right_edge = ie
                                    right_dtheta = dtheta
                            else:
                                # aligned with edge
                                if verbose:
                                    print(
                                        "{}: line is aligned with edge {}".format(j, ie)
                                    )
                                left_edge = ie
                                right_edge = ie
                                break
                        if verbose:
                            print(
                                "{}: the edge to the left is edge {}".format(
                                    j, left_edge
                                )
                            )
                            print(
                                "{}: the edge to the right is edge {}".format(
                                    j, left_edge
                                )
                            )
                        if left_edge == right_edge:
                            if verbose:
                                print("{}: continue along edge {}".format(j, left_edge))
                            index0 = mesh_data.edge_face_connectivity[left_edge, :]
                        else:
                            if verbose:
                                print(
                                    "{}: continue between edges {} on the left and {} on the right".format(
                                        j, left_edge, right_edge
                                    )
                                )
                            left_faces = mesh_data.edge_face_connectivity[left_edge, :]
                            right_faces = mesh_data.edge_face_connectivity[
                                          right_edge, :
                                          ]
                            if (
                                    left_faces[0] in right_faces
                                    and left_faces[1] in right_faces
                            ):
                                # the two edges are shared by two faces ... check first face
                                fn1 = mesh_data.face_node[left_faces[0]]
                                fe1 = mesh_data.face_edge_connectivity[left_faces[0]]
                                if verbose:
                                    print(
                                        "{}: those edges are shared by two faces: {}".format(
                                            j, left_faces
                                        )
                                    )
                                    print(
                                        "{}: face {} has nodes: {}".format(
                                            j, left_faces[0], fn1
                                        )
                                    )
                                    print(
                                        "{}: face {} has edges: {}".format(
                                            j, left_faces[0], fe1
                                        )
                                    )
                                # here we need that the nodes of the face are listed in clockwise order
                                # and that edges[i] is the edge connecting node[i-1] with node[i]
                                # the latter is guaranteed by batch.derive_topology_arrays
                                if fe1[fn1 == node] == right_edge:
                                    index0 = left_faces[0]
                                else:
                                    index0 = left_faces[1]
                            elif left_faces[0] in right_faces:
                                index0 = left_faces[0]
                            elif left_faces[1] in right_faces:
                                index0 = left_faces[1]
                            else:
                                raise Exception(
                                    "Shouldn't come here .... left edge {} and right edge {} don't share any face".format(
                                        left_edge, right_edge
                                    )
                                )

                    elif b[0] == 1:
                        # ending at slice point, so ending on an edge ...
                        if verbose:
                            print("{}: ending on edge {} at {}".format(j, edge, b[0]))
                        # figure out where we will be heading afterwards ...
                        if j == len(bp) - 1:
                            # catch case of last segment
                            if verbose:
                                print("{}: last point ends on an edge".format(j))
                            if ind == crds.shape[0]:
                                crds = enlarge(crds, (ind + 1, 2))
                                idx = enlarge(idx, (ind + 1,))
                            crds[ind] = bpj
                            if index == -2:
                                idx[ind] = vindex[0]
                            else:
                                idx[ind] = index
                            ind += 1
                            break
                        else:
                            # this segment ends on the edge, so check next segment ...
                            # direction of next segment from bpj to bp[j+1]
                            theta = math.atan2(
                                bp[j + 1][1] - bpj[1], bp[j + 1][0] - bp[j][0]
                            )
                        if verbose:
                            print("{}: moving in direction theta = {}".format(j, theta))
                        theta_edge = math.atan2(
                            mesh_data.y_edge_coords[edge, 1]
                            - mesh_data.y_edge_coords[edge, 0],
                            mesh_data.x_edge_coords[edge, 1]
                            - mesh_data.x_edge_coords[edge, 0],
                            )
                        if theta == theta_edge or theta == -theta_edge:
                            # aligned with edge
                            if verbose:
                                print("{}: continue along edge {}".format(j, edge))
                            index0 = faces
                        else:
                            # check whether the (extended) segment slices any edge of faces[0]
                            fe1 = mesh_data.face_edge_connectivity[faces[0]]
                            a, b, edges = _get_slices_core(
                                fe1, mesh_data, bpj, bp[j + 1], 0.0, False
                            )
                            if len(edges) > 0:
                                # yes, a slice (typically 1, but could be 2 if it slices at a node
                                # but that doesn't matter) ... so, we continue towards faces[0]
                                index0 = faces[0]
                            else:
                                # no slice for faces[0], so we must be going in the other direction
                                index0 = faces[1]

                    if index0 is not None:
                        if verbose:
                            if index == -1:
                                print(
                                    "{}: moving from outside via node {} to {} at b = {}".format(
                                        j, node, index0, prev_b
                                    )
                                )
                            elif index == -2:
                                print(
                                    "{}: moving from edge between {} via node {} to {} at b = {}".format(
                                        j, vindex, node, index0, prev_b
                                    )
                                )
                            else:
                                print(
                                    "{}: moving from {} via node {} to {} at b = {}".format(
                                        j, index, node, index0, prev_b
                                    )
                                )

                        if (
                                isinstance(index0, int)
                                or isinstance(index0, np.int32)
                                or isinstance(index0, np.int64)
                        ):
                            index = index0
                        elif len(index0) == 1:
                            index = index0[0]
                        else:
                            index = -2
                            vindex = index0

                    elif faces[0] == index:
                        if verbose:
                            print(
                                "{}: moving from {} via edge {} to {} at b = {}".format(
                                    j, index, edge, faces[1], prev_b
                                )
                            )
                        index = faces[1]
                    elif faces[1] == index:
                        if verbose:
                            print(
                                "{}: moving from {} via edge {} to {} at b = {}".format(
                                    j, index, edge, faces[0], prev_b
                                )
                            )
                        index = faces[0]
                    else:
                        raise Exception(
                            "Shouldn't come here .... index {} differs from both faces {} and {} associated with slicing edge {}".format(
                                index, faces[0], faces[1], edge
                            )
                        )
                    if ind == crds.shape[0]:
                        crds = enlarge(crds, (2 * ind, 2))
                        idx = enlarge(idx, (2 * ind,))
                    crds[ind] = bpj1 + prev_b * (bpj - bpj1)
                    if index == -2:
                        idx[ind] = vindex[0]
                    else:
                        idx[ind] = index
                    ind += 1
                    if prev_b == 1:
                        break

    # clip to actual length (idx refers to segments, so we can ignore the last value)
    crds = crds[:ind]
    idx = idx[: ind - 1]

    # remove tiny segments
    d = np.sqrt((np.diff(crds, axis=0) ** 2).sum(axis=1))
    mask = np.concatenate((np.ones((1), dtype="bool"), d > d_thresh))
    crds = crds[mask, :]
    idx = idx[mask[1:]]

    # since index refers to segments, don't return the first one
    return crds, idx

For more details, see Bank Erosion Mesh Processor.

Debugging Utilities #

`dfastbe.bank_erosion.debugger` #

Bank Erosion Debugger.

`Debugger` #

Class to handle debugging and output of bank erosion calculations.

Source code in src/dfastbe/bank_erosion/debugger.py

class Debugger:
    """Class to handle debugging and output of bank erosion calculations."""

    def __init__(self, crs: str, output_dir: str):
        """Debugger constructor."""
        self.crs = crs
        self.output_dir = output_dir

    def last_discharge_level(
        self,
        bank_index: int,
        single_bank: SingleBank,
        fairway_data: FairwayData,
        erosion_inputs: SingleErosion,
        single_parameters: SingleParameters,
        single_calculation: SingleCalculation,
    ):
        """Write the last discharge level to a shapefile and CSV file."""
        bank_coords_mind = single_bank.get_mid_points()
        params = {
            "chainage": single_bank.bank_chainage_midpoints,
            "x": bank_coords_mind[:, 0],
            "y": bank_coords_mind[:, 1],
            "iface_fw": single_bank.fairway_face_indices,
            "iface_bank": single_bank.bank_face_indices,
            "bank_height": single_bank.height,
            "segment_length": single_bank.segment_length,
            "zw0": fairway_data.fairway_initial_water_levels[bank_index],
            "ship_velocity": single_parameters.ship_velocity,
            "ship_type": single_parameters.ship_type,
            "draught": single_parameters.ship_draught,
            "mu_slp": single_parameters.mu_slope,
            "bank_fairway_dist": single_bank.fairway_distances,
            "fairway_wave_reduction_distance": erosion_inputs.wave_fairway_distance_0,
            "fairway_wave_disappear_distance": erosion_inputs.wave_fairway_distance_1,
            "water_depth_fairway": single_calculation.water_depth,
            "dike_height": erosion_inputs.bank_protection_level,
            "erosion_distance": single_calculation.erosion_distance_eq,
            "erosion_volume": single_calculation.erosion_volume_eq,
        }

        path = f"{str(self.output_dir)}/debug.EQ.B{bank_index + 1}"
        bank_coords_geo = single_bank.get_mid_points(as_geo_series=True, crs=self.crs)
        self._write_data(bank_coords_geo, params, path)

    def middle_levels(
        self,
        bank_ind: int,
        q_level: int,
        single_bank: SingleBank,
        fairway_data: FairwayData,
        erosion_inputs: SingleErosion,
        single_parameters: SingleParameters,
        single_calculation: SingleCalculation,
    ):
        """Write the middle levels to a shapefile and CSV file."""
        bank_coords_mind = single_bank.get_mid_points()
        params = {
            "chainage": single_bank.bank_chainage_midpoints,
            "x": bank_coords_mind[:, 0],
            "y": bank_coords_mind[:, 1],
            "iface_fw": single_bank.fairway_face_indices,
            "iface_bank": single_bank.bank_face_indices,
            "velocity": single_calculation.bank_velocity,
            "bank_height": single_bank.height,
            "segment_length": single_bank.segment_length,
            "zw": single_calculation.water_level,
            "zw0": fairway_data.fairway_initial_water_levels[bank_ind],
            "tauc": erosion_inputs.tauc,
            "num_ship": single_parameters.num_ship,
            "ship_velocity": single_parameters.ship_velocity,
            "num_waves_per_ship": single_parameters.num_waves_per_ship,
            "ship_type": single_parameters.ship_type,
            "draught": single_parameters.ship_draught,
            "mu_slp": single_parameters.mu_slope,
            "mu_reed": single_parameters.mu_reed,
            "dist_fw": single_bank.fairway_distances,
            "fairway_wave_reduction_distance": erosion_inputs.wave_fairway_distance_0,
            "fairway_wave_disappear_distance": erosion_inputs.wave_fairway_distance_1,
            "water_depth_fairway": single_calculation.water_depth,
            "chez": single_calculation.chezy,
            "dike_height": erosion_inputs.bank_protection_level,
            "erosion_distance": single_calculation.erosion_distance_tot,
            "erosion_volume": single_calculation.erosion_volume_tot,
            "erosion_distance_shipping": single_calculation.erosion_distance_shipping,
            "erosion_distance_flow": single_calculation.erosion_distance_flow,
        }
        path = f"{str(self.output_dir)}/debug.Q{q_level + 1}.B{bank_ind + 1}"
        bank_coords_geo = single_bank.get_mid_points(as_geo_series=True, crs=self.crs)
        self._write_data(bank_coords_geo, params, path)

    @staticmethod
    def _write_data(coords: GeoSeries, data: Dict[str, np.ndarray], path: str):
        """Write the data to a shapefile and CSV file."""
        csv_path = f"{path}.csv"
        shp_path = f"{path}.shp"
        _write_shp(coords, data, shp_path)
        _write_csv(data, csv_path)

`init(crs: str, output_dir: str)` #

Debugger constructor.

Source code in src/dfastbe/bank_erosion/debugger.py

def __init__(self, crs: str, output_dir: str):
    """Debugger constructor."""
    self.crs = crs
    self.output_dir = output_dir

`last_discharge_level(bank_index: int, single_bank: SingleBank, fairway_data: FairwayData, erosion_inputs: SingleErosion, single_parameters: SingleParameters, single_calculation: SingleCalculation)` #

Write the last discharge level to a shapefile and CSV file.

Source code in src/dfastbe/bank_erosion/debugger.py

def last_discharge_level(
    self,
    bank_index: int,
    single_bank: SingleBank,
    fairway_data: FairwayData,
    erosion_inputs: SingleErosion,
    single_parameters: SingleParameters,
    single_calculation: SingleCalculation,
):
    """Write the last discharge level to a shapefile and CSV file."""
    bank_coords_mind = single_bank.get_mid_points()
    params = {
        "chainage": single_bank.bank_chainage_midpoints,
        "x": bank_coords_mind[:, 0],
        "y": bank_coords_mind[:, 1],
        "iface_fw": single_bank.fairway_face_indices,
        "iface_bank": single_bank.bank_face_indices,
        "bank_height": single_bank.height,
        "segment_length": single_bank.segment_length,
        "zw0": fairway_data.fairway_initial_water_levels[bank_index],
        "ship_velocity": single_parameters.ship_velocity,
        "ship_type": single_parameters.ship_type,
        "draught": single_parameters.ship_draught,
        "mu_slp": single_parameters.mu_slope,
        "bank_fairway_dist": single_bank.fairway_distances,
        "fairway_wave_reduction_distance": erosion_inputs.wave_fairway_distance_0,
        "fairway_wave_disappear_distance": erosion_inputs.wave_fairway_distance_1,
        "water_depth_fairway": single_calculation.water_depth,
        "dike_height": erosion_inputs.bank_protection_level,
        "erosion_distance": single_calculation.erosion_distance_eq,
        "erosion_volume": single_calculation.erosion_volume_eq,
    }

    path = f"{str(self.output_dir)}/debug.EQ.B{bank_index + 1}"
    bank_coords_geo = single_bank.get_mid_points(as_geo_series=True, crs=self.crs)
    self._write_data(bank_coords_geo, params, path)

`middle_levels(bank_ind: int, q_level: int, single_bank: SingleBank, fairway_data: FairwayData, erosion_inputs: SingleErosion, single_parameters: SingleParameters, single_calculation: SingleCalculation)` #

Write the middle levels to a shapefile and CSV file.

Source code in src/dfastbe/bank_erosion/debugger.py

def middle_levels(
    self,
    bank_ind: int,
    q_level: int,
    single_bank: SingleBank,
    fairway_data: FairwayData,
    erosion_inputs: SingleErosion,
    single_parameters: SingleParameters,
    single_calculation: SingleCalculation,
):
    """Write the middle levels to a shapefile and CSV file."""
    bank_coords_mind = single_bank.get_mid_points()
    params = {
        "chainage": single_bank.bank_chainage_midpoints,
        "x": bank_coords_mind[:, 0],
        "y": bank_coords_mind[:, 1],
        "iface_fw": single_bank.fairway_face_indices,
        "iface_bank": single_bank.bank_face_indices,
        "velocity": single_calculation.bank_velocity,
        "bank_height": single_bank.height,
        "segment_length": single_bank.segment_length,
        "zw": single_calculation.water_level,
        "zw0": fairway_data.fairway_initial_water_levels[bank_ind],
        "tauc": erosion_inputs.tauc,
        "num_ship": single_parameters.num_ship,
        "ship_velocity": single_parameters.ship_velocity,
        "num_waves_per_ship": single_parameters.num_waves_per_ship,
        "ship_type": single_parameters.ship_type,
        "draught": single_parameters.ship_draught,
        "mu_slp": single_parameters.mu_slope,
        "mu_reed": single_parameters.mu_reed,
        "dist_fw": single_bank.fairway_distances,
        "fairway_wave_reduction_distance": erosion_inputs.wave_fairway_distance_0,
        "fairway_wave_disappear_distance": erosion_inputs.wave_fairway_distance_1,
        "water_depth_fairway": single_calculation.water_depth,
        "chez": single_calculation.chezy,
        "dike_height": erosion_inputs.bank_protection_level,
        "erosion_distance": single_calculation.erosion_distance_tot,
        "erosion_volume": single_calculation.erosion_volume_tot,
        "erosion_distance_shipping": single_calculation.erosion_distance_shipping,
        "erosion_distance_flow": single_calculation.erosion_distance_flow,
    }
    path = f"{str(self.output_dir)}/debug.Q{q_level + 1}.B{bank_ind + 1}"
    bank_coords_geo = single_bank.get_mid_points(as_geo_series=True, crs=self.crs)
    self._write_data(bank_coords_geo, params, path)

For more details, see Bank Erosion Debugger.

Parameters:

Name	Type	Description	Default
`is_right_bank`	`List[bool]`	List indicating if the bank is right or not.	required
`bank_chainage_midpoints`	`List[ndarray]`	River chainage for the midpoints of each segment of the bank line	required
`bank_line_coords`	`List[ndarray]`	Coordinates of the bank lines.	required
`bank_face_indices`	`List[ndarray]`	Indices of the faces associated with the banks.	required
`bank_lines`	`GeoDataFrame`	GeoDataFrame containing the bank lines.	`GeoDataFrame()`
`n_bank_lines`	`int`	Number of bank lines.	`0`
`bank_line_size`	`List[ndarray]`	Size of each individual bank line.	required
`fairway_distances`	`List[ndarray]`	The distance of each bank line point to the closest fairway point.	required
`fairway_face_indices`	`List[ndarray]`	The face index of the closest fairway point for each bank line point.	required

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class BankData(BaseBank[SingleBank]):
    """Class to hold bank-related data.

    args:
        is_right_bank (List[bool]):
            List indicating if the bank is right or not.
        bank_chainage_midpoints (List[np.ndarray]):
            River chainage for the midpoints of each segment of the bank line
        bank_line_coords (List[np.ndarray]):
            Coordinates of the bank lines.
        bank_face_indices (List[np.ndarray]):
            Indices of the faces associated with the banks.
        bank_lines (GeoDataFrame):
            GeoDataFrame containing the bank lines.
        n_bank_lines (int):
            Number of bank lines.
        bank_line_size (List[np.ndarray]):
            Size of each individual bank line.
        fairway_distances (List[np.ndarray]):
            The distance of each bank line point to the closest fairway point.
        fairway_face_indices (List[np.ndarray]):
            The face index of the closest fairway point for each bank line point.
    """
    bank_lines: GeoDataFrame = field(default_factory=GeoDataFrame)
    n_bank_lines: int = 0

    @classmethod
    def from_column_arrays(
        cls,
        data: dict,
        bank_cls: Type["SingleBank"],
        bank_lines: GeoDataFrame,
        n_bank_lines: int,
        bank_order: Tuple[str, str] = ("left", "right")
    ) -> "BankData":
        # Only include fields that belong to the bank-specific data
        base_fields = {k: v for k, v in data.items() if k != "id"}
        base = BaseBank.from_column_arrays(
            {"id": data.get("id"), **base_fields}, bank_cls, bank_order=bank_order
        )

        return cls(
            id=base.id,
            left=base.left,
            right=base.right,
            bank_lines=bank_lines,
            n_bank_lines=n_bank_lines,
        )

    @property
    def bank_line_coords(self) -> List[np.ndarray]:
        """Get the coordinates of the bank lines."""
        return [self.left.bank_line_coords, self.right.bank_line_coords]

    @property
    def is_right_bank(self) -> List[bool]:
        """Get the bank direction."""
        return [self.left.is_right_bank, self.right.is_right_bank]

    @property
    def bank_chainage_midpoints(self) -> List[np.ndarray]:
        """Get the chainage midpoints of the bank lines."""
        return [self.left.bank_chainage_midpoints, self.right.bank_chainage_midpoints]

    @property
    def num_stations_per_bank(self) -> List[int]:
        """Get the number of stations per bank."""
        return [self.left.length, self.right.length]

    @property
    def height(self) -> List[np.ndarray]:
        """Get the bank height."""
        return [self.left.height, self.right.height]

@dataclass
class BaseBank(Generic[GenericType]):
    left: GenericType
    right: GenericType
    id: Optional[int] = field(default=None)

    def get_bank(self, bank_index: int) -> GenericType:
        if bank_index == 0:
            return self.left
        elif bank_index == 1:
            return self.right
        else:
            raise ValueError("bank_index must be 0 (left) or 1 (right)")

    @classmethod
    def from_column_arrays(
        cls: Type["BaseBank[GenericType]"],
        data: Dict[str, Any],
        bank_cls: Type[GenericType],
        bank_order: Tuple[str, str] = ("left", "right")
    ) -> "BaseBank[GenericType]":
        if set(bank_order) != {"left", "right"}:
            raise ValueError("bank_order must be a permutation of ('left', 'right')")

        id_val = data.get("id")

        # Extract the first and second array for each parameter (excluding id)
        first_args = {}
        second_args = {}
        for key, value in data.items():
            if key == "id":
                continue
            if not isinstance(value, list) or len(value) != 2:
                raise ValueError(f"Expected 2-column array for key '{key}', got shape {value.shape}")

            split = dict(zip(bank_order, value))
            first_args[key] = split["left"]
            second_args[key] = split["right"]

        left = bank_cls(**first_args)
        right = bank_cls(**second_args)

        return cls(id=id_val, left=left, right=right)

    def __iter__(self) -> Iterator[GenericType]:
        """Iterate over the banks."""
        return iter([self.left, self.right])

`iter() -> Iterator[GenericType]` #

Iterate over the banks.

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

def __iter__(self) -> Iterator[GenericType]:
    """Iterate over the banks."""
    return iter([self.left, self.right])

`DischargeLevels` #

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

class DischargeLevels:

    def __init__(self, levels: List[SingleDischargeLevel]):
        self.levels = levels

    def __getitem__(self, index: int) -> SingleDischargeLevel:
        return self.levels[index]

    def __len__(self) -> int:
        return len(self.levels)

    def append(self, level_calc: SingleDischargeLevel):
        self.levels.append(level_calc)

    def get_max_hfw_level(self) -> float:
        return max(level.hfw_max for level in self.levels)

    def total_erosion_volume(self) -> float:
        return sum(
            np.sum(level.left.erosion_volume_tot) + np.sum(level.right.erosion_volume_tot)
            for level in self.levels
        )

    def __iter__(self):
        return iter(self.levels)

    def accumulate(self, attribute_name: str, bank_side: Union[str, List[str]] = None) -> List[np.ndarray]:
        if bank_side is None:
            bank_side = ["left", "right"]
        elif isinstance(bank_side, str):
            bank_side = [bank_side]

        if not all(side in ["left", "right"] for side in bank_side):
            raise ValueError("bank_side must be 'left', 'right', or a list of these.")

        total = [
            self._accumulate_attribute_side(attribute_name, side) for side in bank_side
        ]
        return total

    def _accumulate_attribute_side(self, attribute_name: str, bank_side: str) -> np.ndarray:
        for i, level in enumerate(self.levels):
            bank = getattr(level, bank_side)
            attr = getattr(bank, attribute_name, None)
            if attr is None:
                raise AttributeError(f"{attribute_name} not found in {bank_side} bank of level with id={level.id}")
            if i == 0:
                total = attr
            else:
                total += attr
        return total

    def _get_attr_both_sides_level(self, attribute_name: str, level) -> List[np.ndarray]:
        """Get the attributes of the levels for both left and right bank."""
        sides = [getattr(self.levels[level], side) for side in ["left", "right"]]
        attr = [getattr(side, attribute_name, None) for side in sides]
        return attr

    def get_attr_level(self, attribute_name: str) -> List[List[np.ndarray]]:
        """Get the attributes of the levels for both left and right bank."""
        return [self._get_attr_both_sides_level(attribute_name, level) for level in range(len(self.levels))]

    def get_water_level_data(self) -> WaterLevelData:
        return WaterLevelData(
            hfw_max=self.levels[-1].hfw_max,
            water_level=self.get_attr_level("water_level"),
            ship_wave_max=self.get_attr_level("ship_wave_max"),
            ship_wave_min=self.get_attr_level("ship_wave_min"),
            velocity=self.get_attr_level("bank_velocity"),
            chezy=self.get_attr_level("chezy"),
            vol_per_discharge=self.get_attr_level("volume_per_discharge"),
        )

`get_attr_level(attribute_name: str) -> List[List[np.ndarray]]` #

Get the attributes of the levels for both left and right bank.

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

def get_attr_level(self, attribute_name: str) -> List[List[np.ndarray]]:
    """Get the attributes of the levels for both left and right bank."""
    return [self._get_attr_both_sides_level(attribute_name, level) for level in range(len(self.levels))]

`ErosionInputs` `dataclass` #

Bases: BaseBank[SingleErosion]

Class to hold erosion inputs.

Parameters:

Name	Type	Description	Default
`shipping_data`	`Dict[str, ndarray]`	Data on all the vessels that travel through the river.	`dict()`
`wave_fairway_distance_0`	`List[ndarray]`	Threshold fairway distance 0 for wave attenuation.	required
`wave_fairway_distance_1`	`List[ndarray]`	Threshold fairway distance 1 for wave attenuation.	required
`bank_protection_level`	`List[ndarray]`	Bank protection level.	required
`tauc`	`List[ndarray]`	Critical bank shear stress values.	required
`bank_type`	`List[ndarray]`	Integer representation of the bank type. Represents an index into the taucls_str array.	`lambda: array([])()`
`taucls`	`ndarray`	Critical bank shear stress values for different bank types.	required
`taucls_str`	`Tuple[str]`	String representation for different bank types.	required

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class ErosionInputs(BaseBank[SingleErosion]):
    """Class to hold erosion inputs.

    args:
        shipping_data (Dict[str, np.ndarray]):
            Data on all the vessels that travel through the river.
        wave_fairway_distance_0 (List[np.ndarray]):
            Threshold fairway distance 0 for wave attenuation.
        wave_fairway_distance_1 (List[np.ndarray]):
            Threshold fairway distance 1 for wave attenuation.
        bank_protection_level (List[np.ndarray]):
            Bank protection level.
        tauc (List[np.ndarray]):
            Critical bank shear stress values.
        bank_type (List[np.ndarray]):
            Integer representation of the bank type. Represents an index into the taucls_str array.
        taucls (np.ndarray):
            Critical bank shear stress values for different bank types.
        taucls_str (Tuple[str]):
            String representation for different bank types.
    """
    shipping_data: Dict[str, List[np.ndarray]] = field(default_factory=dict)
    bank_type: np.ndarray = field(default_factory=lambda: np.array([]))
    taucls: ClassVar[np.ndarray] = np.array([1e20, 95, 3.0, 0.95, 0.15])
    taucls_str: ClassVar[Tuple[str]] = (
        "protected",
        "vegetation",
        "good clay",
        "moderate/bad clay",
        "sand",
    )

    @classmethod
    def from_column_arrays(
        cls, data: dict, bank_cls: Type["SingleErosion"], shipping_data: Dict[str, List[np.ndarray]],
        bank_type: np.ndarray, bank_order: Tuple[str, str] = ("left", "right")
    ) -> "ErosionInputs":
        # Only include fields that belong to the bank-specific data
        base_fields = {k: v for k, v in data.items() if k != "id"}
        base = BaseBank.from_column_arrays(
            {"id": data.get("id"), **base_fields}, bank_cls, bank_order=bank_order
        )

        return cls(
            id=base.id,
            left=base.left,
            right=base.right,
            shipping_data=shipping_data,
            bank_type=bank_type,
        )

    @property
    def bank_protection_level(self) -> List[np.ndarray]:
        """Get the bank protection level."""
        return [self.left.bank_protection_level, self.right.bank_protection_level]

    @property
    def tauc(self) -> List[np.ndarray]:
        """Get the critical bank shear stress values."""
        return [self.left.tauc, self.right.tauc]

Name	Type	Description	Default
`eq_erosion_dist`	`List[ndarray]`	Erosion distance at equilibrium for each bank line.	required
`total_erosion_dist`	`List[ndarray]`	Total erosion distance for each bank line.	required
`flow_erosion_dist`	`List[ndarray]`	Total erosion distance caused by flow for each bank line.	required
`ship_erosion_dist`	`List[ndarray]`	Total erosion distance caused by ship waves for each bank line.	required
`eq_eroded_vol`	`List[ndarray]`	Eroded volume at equilibrium for each bank line.	required
`total_eroded_vol`	`List[ndarray]`	Total eroded volume for each bank line.	required
`erosion_time`	`int`	Time over which erosion is calculated.	required
`avg_erosion_rate`	`ndarray`	Average erosion rate data.	`lambda: empty(0)()`
`eq_eroded_vol_per_km`	`ndarray`	Equilibrium eroded volume calculated per kilometer bin.	`lambda: empty(0)()`
`total_eroded_vol_per_km`	`ndarray`	Total eroded volume calculated per kilometer bin.	`lambda: empty(0)()`

Examples:

You can create an instance of the ErosionResults class as follows:

>>> from dfastbe.bank_erosion.data_models.calculation import ErosionResults
>>> import numpy as np
>>> erosion_results = ErosionResults(
...     eq_erosion_dist=[np.array([0.1, 0.2])],
...     total_erosion_dist=[np.array([0.3, 0.4])],
...     flow_erosion_dist=[np.array([0.5, 0.6])],
...     ship_erosion_dist=[np.array([0.7, 0.8])],
...     eq_eroded_vol=[np.array([1.1, 1.2])],
...     total_eroded_vol=[np.array([1.3, 1.4])],
...     erosion_time=10,
...     avg_erosion_rate=np.array([0.1, 0.2]),
...     eq_eroded_vol_per_km=np.array([0.3, 0.4]),
...     total_eroded_vol_per_km=np.array([0.5, 0.6]),
... )
>>> print(erosion_results)
ErosionResults(eq_erosion_dist=[array([0.1, 0.2])], total_erosion_dist=[array([0.3, 0.4])], flow_erosion_dist=[array([0.5, 0.6])], ship_erosion_dist=[array([0.7, 0.8])], eq_eroded_vol=[array([1.1, 1.2])], total_eroded_vol=[array([1.3, 1.4])], erosion_time=10, avg_erosion_rate=array([0.1, 0.2]), eq_eroded_vol_per_km=array([0.3, 0.4]), total_eroded_vol_per_km=array([0.5, 0.6]))

The avg_erosion_rate, eq_eroded_vol_per_km, and total_eroded_vol_per_km attributes are optional and can be set to empty arrays if not needed.

>>> from dfastbe.bank_erosion.data_models.calculation import ErosionResults
>>> import numpy as np
>>> erosion_results = ErosionResults(
...     eq_erosion_dist=[np.array([0.1, 0.2])],
...     total_erosion_dist=[np.array([0.3, 0.4])],
...     flow_erosion_dist=[np.array([0.5, 0.6])],
...     ship_erosion_dist=[np.array([0.7, 0.8])],
...     eq_eroded_vol=[np.array([1.1, 1.2])],
...     total_eroded_vol=[np.array([1.3, 1.4])],
...     erosion_time=10,
... )
>>> print(erosion_results)
ErosionResults(eq_erosion_dist=[array([0.1, 0.2])], total_erosion_dist=[array([0.3, 0.4])], flow_erosion_dist=[array([0.5, 0.6])], ship_erosion_dist=[array([0.7, 0.8])], eq_eroded_vol=[array([1.1, 1.2])], total_eroded_vol=[array([1.3, 1.4])], erosion_time=10, avg_erosion_rate=array([], dtype=float64), eq_eroded_vol_per_km=array([], dtype=float64), total_eroded_vol_per_km=array([], dtype=float64))

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class ErosionResults:
    """Class to hold erosion results.

    args:
        eq_erosion_dist (List[np.ndarray]):
            Erosion distance at equilibrium for each bank line.
        total_erosion_dist (List[np.ndarray]):
            Total erosion distance for each bank line.
        flow_erosion_dist (List[np.ndarray]):
            Total erosion distance caused by flow for each bank line.
        ship_erosion_dist (List[np.ndarray]):
            Total erosion distance caused by ship waves for each bank line.
        eq_eroded_vol (List[np.ndarray]):
            Eroded volume at equilibrium for each bank line.
        total_eroded_vol (List[np.ndarray]):
            Total eroded volume for each bank line.
        erosion_time (int):
            Time over which erosion is calculated.
        avg_erosion_rate (np.ndarray):
            Average erosion rate data.
        eq_eroded_vol_per_km (np.ndarray):
            Equilibrium eroded volume calculated per kilometer bin.
        total_eroded_vol_per_km (np.ndarray):
            Total eroded volume calculated per kilometer bin.

    Examples:
        - You can create an instance of the ErosionResults class as follows:
        ```python
        >>> from dfastbe.bank_erosion.data_models.calculation import ErosionResults
        >>> import numpy as np
        >>> erosion_results = ErosionResults(
        ...     eq_erosion_dist=[np.array([0.1, 0.2])],
        ...     total_erosion_dist=[np.array([0.3, 0.4])],
        ...     flow_erosion_dist=[np.array([0.5, 0.6])],
        ...     ship_erosion_dist=[np.array([0.7, 0.8])],
        ...     eq_eroded_vol=[np.array([1.1, 1.2])],
        ...     total_eroded_vol=[np.array([1.3, 1.4])],
        ...     erosion_time=10,
        ...     avg_erosion_rate=np.array([0.1, 0.2]),
        ...     eq_eroded_vol_per_km=np.array([0.3, 0.4]),
        ...     total_eroded_vol_per_km=np.array([0.5, 0.6]),
        ... )
        >>> print(erosion_results)
        ErosionResults(eq_erosion_dist=[array([0.1, 0.2])], total_erosion_dist=[array([0.3, 0.4])], flow_erosion_dist=[array([0.5, 0.6])], ship_erosion_dist=[array([0.7, 0.8])], eq_eroded_vol=[array([1.1, 1.2])], total_eroded_vol=[array([1.3, 1.4])], erosion_time=10, avg_erosion_rate=array([0.1, 0.2]), eq_eroded_vol_per_km=array([0.3, 0.4]), total_eroded_vol_per_km=array([0.5, 0.6]))

        ```

        - The `avg_erosion_rate`, `eq_eroded_vol_per_km`, and `total_eroded_vol_per_km` attributes are optional and
        can be set to empty arrays if not needed.

        ```python
        >>> from dfastbe.bank_erosion.data_models.calculation import ErosionResults
        >>> import numpy as np
        >>> erosion_results = ErosionResults(
        ...     eq_erosion_dist=[np.array([0.1, 0.2])],
        ...     total_erosion_dist=[np.array([0.3, 0.4])],
        ...     flow_erosion_dist=[np.array([0.5, 0.6])],
        ...     ship_erosion_dist=[np.array([0.7, 0.8])],
        ...     eq_eroded_vol=[np.array([1.1, 1.2])],
        ...     total_eroded_vol=[np.array([1.3, 1.4])],
        ...     erosion_time=10,
        ... )
        >>> print(erosion_results)
        ErosionResults(eq_erosion_dist=[array([0.1, 0.2])], total_erosion_dist=[array([0.3, 0.4])], flow_erosion_dist=[array([0.5, 0.6])], ship_erosion_dist=[array([0.7, 0.8])], eq_eroded_vol=[array([1.1, 1.2])], total_eroded_vol=[array([1.3, 1.4])], erosion_time=10, avg_erosion_rate=array([], dtype=float64), eq_eroded_vol_per_km=array([], dtype=float64), total_eroded_vol_per_km=array([], dtype=float64))

        ```
    """

    eq_erosion_dist: List[np.ndarray]
    total_erosion_dist: List[np.ndarray]
    flow_erosion_dist: List[np.ndarray]
    ship_erosion_dist: List[np.ndarray]
    eq_eroded_vol: List[np.ndarray]
    total_eroded_vol: List[np.ndarray]
    erosion_time: int
    avg_erosion_rate: np.ndarray = field(default_factory=lambda : np.empty(0))
    eq_eroded_vol_per_km: np.ndarray = field(default_factory=lambda : np.empty(0))
    total_eroded_vol_per_km: np.ndarray = field(default_factory=lambda : np.empty(0))

`FairwayData` `dataclass` #

Class to hold fairway-related data.

Parameters:

Name	Type	Description	Default
`fairway_face_indices`	`ndarray`	Mesh face indices matching to the fairway points.	required
`intersection_coords`	`ndarray`	The x, y coordinates of the intersection points of the fairway with the simulation mesh.	required
`fairway_initial_water_levels`	`List[ndarray]`	Reference water level at the fairway	`list()`

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class FairwayData:
    """Class to hold fairway-related data.

    args:
        fairway_face_indices (np.ndarray):
            Mesh face indices matching to the fairway points.
        intersection_coords (np.ndarray):
            The x, y coordinates of the intersection points of the fairway with the simulation mesh.
        fairway_initial_water_levels (List[np.ndarray]):
            Reference water level at the fairway
    """

    fairway_face_indices: np.ndarray
    intersection_coords: np.ndarray
    fairway_initial_water_levels: List[np.ndarray] = field(default_factory=list)

`MeshData` `dataclass` #

Class to hold mesh-related data.

Parameters:

Name	Type	Description	Default
`x_face_coords`	`ndarray`	X-coordinates of the mesh faces.	required
`y_face_coords`	`ndarray`	Y-coordinates of the mesh faces.	required
`x_edge_coords`	`ndarray`	X-coordinates of the mesh edges.	required
`y_edge_coords`	`ndarray`	Y-coordinates of the mesh edges.	required
`face_node`	`ndarray`	Node connectivity for each face.	required
`n_nodes`	`ndarray`	Number of nodes in the mesh.	required
`edge_node`	`ndarray`	Node connectivity for each edge.	required
`edge_face_connectivity`	`ndarray`	Per edge a list of the indices of the faces on the left and right side of that edge.	required
`face_edge_connectivity`	`ndarray`	Per face a list of indices of the edges that together form the boundary of that face.	required
`boundary_edge_nrs`	`ndarray`	List of edge indices that together form the boundary of the whole mesh.	required

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class MeshData:
    """Class to hold mesh-related data.

    args:
        x_face_coords (np.ndarray):
            X-coordinates of the mesh faces.
        y_face_coords (np.ndarray):
            Y-coordinates of the mesh faces.
        x_edge_coords (np.ndarray):
            X-coordinates of the mesh edges.
        y_edge_coords (np.ndarray):
            Y-coordinates of the mesh edges.
        face_node (np.ndarray):
            Node connectivity for each face.
        n_nodes (np.ndarray):
            Number of nodes in the mesh.
        edge_node (np.ndarray):
            Node connectivity for each edge.
        edge_face_connectivity (np.ndarray):
            Per edge a list of the indices of the faces on the left and right side of that edge.
        face_edge_connectivity (np.ndarray):
            Per face a list of indices of the edges that together form the boundary of that face.
        boundary_edge_nrs (np.ndarray):
            List of edge indices that together form the boundary of the whole mesh.
    """

    x_face_coords: np.ndarray
    y_face_coords: np.ndarray
    x_edge_coords: np.ndarray
    y_edge_coords: np.ndarray
    face_node: np.ndarray
    n_nodes: np.ndarray
    edge_node: np.ndarray
    edge_face_connectivity: np.ndarray
    face_edge_connectivity: np.ndarray
    boundary_edge_nrs: np.ndarray

`SingleBank` `dataclass` #

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class SingleBank:
    is_right_bank: bool
    bank_line_coords: np.ndarray
    bank_face_indices: np.ndarray
    bank_line_size: np.ndarray = field(default_factory=lambda: np.array([]))
    fairway_distances: np.ndarray = field(default_factory=lambda: np.array([]))
    fairway_face_indices: np.ndarray = field(default_factory=lambda: np.array([]))
    bank_chainage_midpoints: np.ndarray = field(default_factory=lambda: np.array([]))

    segment_length: np.ndarray = field(init=False)
    dx: np.ndarray = field(init=False)
    dy: np.ndarray = field(init=False)
    length: int = field(init=False)

    # bank height is calculated at the first discharge level only.
    height: Optional[np.ndarray] = field(default=lambda : np.array([]))

    def __post_init__(self):
        """Post-initialization to ensure bank_line_coords is a list of numpy arrays."""
        self.segment_length = self._segment_length()
        self.dx = self._dx()
        self.dy = self._dy()
        self.length = len(self.bank_chainage_midpoints)

    def _segment_length(self) -> np.ndarray:
        """Calculate the length of each segment in the bank line.

        Returns:
            List[np.ndarray]: Length of each segment in the bank line.
        """
        return np.linalg.norm(np.diff(self.bank_line_coords, axis=0), axis=1)

    def _dx(self) -> np.ndarray:
        """Calculate the distance between each bank line point.

        Returns:
            List[np.ndarray]: Distance to the closest fairway point for each bank line point.
        """
        return np.diff(self.bank_line_coords[:, 0])

    def _dy(self) -> np.ndarray:
        """Calculate the distance between each bank line point.

        Returns:
            List[np.ndarray]: Distance to the closest fairway point for each bank line point.
        """
        return np.diff(self.bank_line_coords[:, 1])

    def get_mid_points(self, as_geo_series: bool = False, crs: str = None) -> Union[GeoSeries, np.ndarray]:
        """Band line midpoints.

        Args:
            as_geo_series (bool):
                bool indicating if the output should be a GeoSeries or not.
            crs (str):
                coordinate reference system.
        Returns:
            the midpoints of the bank line coordinates as a GeoSeries or numpy array.
        """
        bank_coords = self.bank_line_coords
        bank_coords_mind = (bank_coords[:-1] + bank_coords[1:]) / 2

        if as_geo_series:
            bank_coords_mind = [Point(xy) for xy in bank_coords_mind]
            bank_coords_mind = GeoSeries(bank_coords_mind, crs=crs)
        return bank_coords_mind

`__post_init__()` #

Post-initialization to ensure bank_line_coords is a list of numpy arrays.

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

def __post_init__(self):
    """Post-initialization to ensure bank_line_coords is a list of numpy arrays."""
    self.segment_length = self._segment_length()
    self.dx = self._dx()
    self.dy = self._dy()
    self.length = len(self.bank_chainage_midpoints)

`get_mid_points(as_geo_series: bool = False, crs: str = None) -> Union[GeoSeries, np.ndarray]` #

Band line midpoints.

Parameters:

Name	Type	Description	Default
`as_geo_series`	`bool`	bool indicating if the output should be a GeoSeries or not.	`False`
`crs`	`str`	coordinate reference system.	`None`

Returns: the midpoints of the bank line coordinates as a GeoSeries or numpy array.

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

def get_mid_points(self, as_geo_series: bool = False, crs: str = None) -> Union[GeoSeries, np.ndarray]:
    """Band line midpoints.

    Args:
        as_geo_series (bool):
            bool indicating if the output should be a GeoSeries or not.
        crs (str):
            coordinate reference system.
    Returns:
        the midpoints of the bank line coordinates as a GeoSeries or numpy array.
    """
    bank_coords = self.bank_line_coords
    bank_coords_mind = (bank_coords[:-1] + bank_coords[1:]) / 2

    if as_geo_series:
        bank_coords_mind = [Point(xy) for xy in bank_coords_mind]
        bank_coords_mind = GeoSeries(bank_coords_mind, crs=crs)
    return bank_coords_mind

`WaterLevelData` `dataclass` #

Class to hold water level data.

Parameters:

Name	Type	Description	Default
`hfw_max`	`float`	Maximum water depth along the fairway.	required
`water_level`	`List[List[ndarray]]`	Water level data.	required
`ship_wave_max`	`List[List[ndarray]]`	Maximum bank height subject to ship waves [m]	required
`ship_wave_min`	`List[List[ndarray]]`	Minimum bank height subject to ship waves [m]	required
`velocity`	`List[List[ndarray]]`	Flow velocity magnitude along the bank [m/s]	required
`bank_height`	`List[ndarray]`	Bank height data.	required
`chezy`	`List[List[ndarray]]`	Chezy coefficient data.	required
`vol_per_discharge`	`List[List[ndarray]]`	Eroded volume per discharge level for each bank line.	required

Source code in src/dfastbe/bank_erosion/data_models/calculation.py

@dataclass
class WaterLevelData:
    """Class to hold water level data.

    args:
        hfw_max (float): Maximum water depth along the fairway.
        water_level (List[List[np.ndarray]]): Water level data.
        ship_wave_max (List[List[np.ndarray]]): Maximum bank height subject to ship waves [m]
        ship_wave_min (List[List[np.ndarray]]): Minimum bank height subject to ship waves [m]
        velocity (List[List[np.ndarray]]): Flow velocity magnitude along the bank [m/s]
        bank_height (List[np.ndarray]): Bank height data.
        chezy (List[List[np.ndarray]]): Chezy coefficient data.
        vol_per_discharge (List[List[np.ndarray]]):
            Eroded volume per discharge level for each bank line.
    """

    hfw_max: float
    water_level: List[List[np.ndarray]]
    ship_wave_max: List[List[np.ndarray]]
    ship_wave_min: List[List[np.ndarray]]
    velocity: List[List[np.ndarray]]
    chezy: List[List[np.ndarray]]
    vol_per_discharge: List[List[np.ndarray]]

For more details, see Bank Erosion Calculation Data Models.

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

class BankLinesResultsError(Exception):
    """Custom exception for BankLine results errors."""

    pass

`ErosionRiverData` #

Bases: BaseRiverData

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

class ErosionRiverData(BaseRiverData):

    def __init__(self, config_file: ConfigFile):
        super().__init__(config_file)
        self.bank_dir = self._get_bank_line_dir()
        self.output_dir = config_file.get_output_dir("erosion")
        self.debug = config_file.debug
        # set plotting flags
        self.plot_flags = config_file.get_plotting_flags(config_file.root_dir)
        # get filter settings for bank levels and flow velocities along banks
        self.zb_dx = config_file.get_float("Erosion", "BedFilterDist", 0.0, positive=True)
        self.vel_dx = config_file.get_float("Erosion", "VelFilterDist", 0.0, positive=True)
        log_text("get_levels")
        self.num_discharge_levels = config_file.get_int("Erosion", "NLevel")
        self.output_intervals = config_file.get_float("Erosion", "OutputInterval", 1.0)
        self.bank_lines = config_file.read_bank_lines(str(self.bank_dir))
        self.river_axis = self._read_river_axis()
        self.erosion_time = self.config_file.get_int("Erosion", "TErosion", positive=True)

    def simulation_data(self) -> ErosionSimulationData:
        """Simulation Data."""
        ref_level = self.config_file.get_int("Erosion", "RefLevel") - 1
        # read simulation data (get_sim_data)
        sim_file = self.config_file.get_sim_file("Erosion", str(ref_level + 1))
        log_text("-")
        log_text("read_simdata", data={"file": sim_file})
        log_text("-")
        simulation_data = ErosionSimulationData.read(sim_file)

        return simulation_data

    def _get_bank_output_dir(self) -> Path:
        bank_output_dir = self.config_file.get_str("General", "BankDir")
        log_text("bank_dir_out", data={"dir": bank_output_dir})
        if os.path.exists(bank_output_dir):
            log_text("overwrite_dir", data={"dir": bank_output_dir})
        else:
            os.makedirs(bank_output_dir)

        return Path(bank_output_dir)

    def _get_bank_line_dir(self) -> Path:
        bank_dir = self.config_file.get_str("General", "BankDir")
        log_text("bank_dir_in", data={"dir": bank_dir})
        bank_dir = Path(bank_dir)
        if not bank_dir.exists():
            log_text("missing_dir", data={"dir": bank_dir})
            raise BankLinesResultsError(
                f"Required bank line directory:{bank_dir} does not exist. please use the banklines command to run the "
                "bankline detection tool first it."
            )
        else:
            return bank_dir

    def _read_river_axis(self) -> LineString:
        """Get the river axis from the analysis settings."""
        river_axis_file = self.config_file.get_str("Erosion", "RiverAxis")
        log_text("read_river_axis", data={"file": river_axis_file})
        river_axis = XYCModel.read(river_axis_file)
        return river_axis

`simulation_data() -> ErosionSimulationData` #

Simulation Data.

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

def simulation_data(self) -> ErosionSimulationData:
    """Simulation Data."""
    ref_level = self.config_file.get_int("Erosion", "RefLevel") - 1
    # read simulation data (get_sim_data)
    sim_file = self.config_file.get_sim_file("Erosion", str(ref_level + 1))
    log_text("-")
    log_text("read_simdata", data={"file": sim_file})
    log_text("-")
    simulation_data = ErosionSimulationData.read(sim_file)

    return simulation_data

`ErosionSimulationData` #

Bases: BaseSimulationData

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

class ErosionSimulationData(BaseSimulationData):

    def compute_mesh_topology(self) -> MeshData:
        """Derive secondary topology arrays from the face-node connectivity of the mesh.

        This function computes the edge-node, edge-face, and face-edge connectivity arrays,
        as well as the boundary edges of the mesh, based on the face-node connectivity provided
        in the simulation data.

        Returns:
            MeshData: a dataclass containing the following attributes:
                - `x_face_coords`: x-coordinates of face nodes
                - `y_face_coords`: y-coordinates of face nodes
                - `x_edge_coords`: x-coordinates of edge nodes
                - `y_edge_coords`: y-coordinates of edge nodes
                - `face_node`: the node indices for each of the mesh faces.
                - `n_nodes`: number of nodes per face
                - `edge_node`: the node indices for each of the mesh edges.
                - `edge_face_connectivity`: the face indices for each of the mesh edge
                - `face_edge_connectivity`: the edge indices for each of the mesh face
                - `boundary_edge_nrs`: indices of boundary edges

        Raises:
            KeyError:
                If required keys (e.g., `face_node`, `nnodes`, `x_node`, `y_node`) are missing from the `sim` object.

        Notes:
            - The function identifies unique edges by sorting and comparing node indices.
            - Boundary edges are identified as edges that belong to only one face.
            - The function assumes that the mesh is well-formed, with consistent face-node connectivity.
        """

        # get a sorted list of edge node connections (shared edges occur twice)
        # face_nr contains the face index to which the edge belongs
        n_faces = self.face_node.shape[0]
        n_edges = sum(self.n_nodes)
        edge_node = np.zeros((n_edges, 2), dtype=int)
        face_nr = np.zeros((n_edges,), dtype=int)
        i = 0
        for face_i in range(n_faces):
            num_edges = self.n_nodes[face_i]  # note: nEdges = nNodes
            for edge_i in range(num_edges):
                if edge_i == 0:
                    edge_node[i, 1] = self.face_node[face_i, num_edges - 1]
                else:
                    edge_node[i, 1] = self.face_node[face_i, edge_i - 1]
                edge_node[i, 0] = self.face_node[face_i, edge_i]
                face_nr[i] = face_i
                i = i + 1
        edge_node.sort(axis=1)
        i2 = np.argsort(edge_node[:, 1], kind="stable")
        i1 = np.argsort(edge_node[i2, 0], kind="stable")
        i12 = i2[i1]
        edge_node = edge_node[i12, :]
        face_nr = face_nr[i12]

        # detect which edges are equal to the previous edge, and get a list of all unique edges
        numpy_true = np.array([True])
        equal_to_previous = np.concatenate(
            (~numpy_true, (np.diff(edge_node, axis=0) == 0).all(axis=1))
        )
        unique_edge = ~equal_to_previous
        n_unique_edges = np.sum(unique_edge)
        # reduce the edge node connections to only the unique edges
        edge_node = edge_node[unique_edge, :]

        # number the edges
        edge_nr = np.zeros(n_edges, dtype=int)
        edge_nr[unique_edge] = np.arange(n_unique_edges, dtype=int)
        edge_nr[equal_to_previous] = edge_nr[
            np.concatenate((equal_to_previous[1:], equal_to_previous[:1]))
        ]

        # if two consecutive edges are unique, the first one occurs only once and represents a boundary edge
        is_boundary_edge = unique_edge & np.concatenate((unique_edge[1:], numpy_true))
        boundary_edge_nrs = edge_nr[is_boundary_edge]

        # go back to the original face order
        edge_nr_in_face_order = np.zeros(n_edges, dtype=int)
        edge_nr_in_face_order[i12] = edge_nr
        # create the face edge connectivity array
        face_edge_connectivity = np.zeros(self.face_node.shape, dtype=int)

        i = 0
        for face_i in range(n_faces):
            num_edges = self.n_nodes[face_i]  # note: num_edges = n_nodes
            for edge_i in range(num_edges):
                face_edge_connectivity[face_i, edge_i] = edge_nr_in_face_order[i]
                i = i + 1

        # determine the edge face connectivity
        edge_face = -np.ones((n_unique_edges, 2), dtype=int)
        edge_face[edge_nr[unique_edge], 0] = face_nr[unique_edge]
        edge_face[edge_nr[equal_to_previous], 1] = face_nr[equal_to_previous]

        x_face_coords = self.apply_masked_indexing(
            self.x_node, self.face_node
        )
        y_face_coords = self.apply_masked_indexing(
            self.y_node, self.face_node
        )
        x_edge_coords = self.x_node[edge_node]
        y_edge_coords = self.y_node[edge_node]

        return MeshData(
            x_face_coords=x_face_coords,
            y_face_coords=y_face_coords,
            x_edge_coords=x_edge_coords,
            y_edge_coords=y_edge_coords,
            face_node=self.face_node,
            n_nodes=self.n_nodes,
            edge_node=edge_node,
            edge_face_connectivity=edge_face,
            face_edge_connectivity=face_edge_connectivity,
            boundary_edge_nrs=boundary_edge_nrs,
        )

    @staticmethod
    def apply_masked_indexing(
        x0: np.array, idx: np.ma.masked_array
    ) -> np.ma.masked_array:
        """
        Index one array by another transferring the mask.

        Args:
            x0 : np.ndarray
                A linear array.
            idx : np.ma.masked_array
                An index array with possibly masked indices.

        returns:
            x1: np.ma.masked_array
                An array with same shape as idx, with mask.
        """
        idx_safe = idx.copy()
        idx_safe.data[np.ma.getmask(idx)] = 0
        x1 = np.ma.masked_where(np.ma.getmask(idx), x0[idx_safe])
        return x1

    def calculate_bank_velocity(self, single_bank: "SingleBank", vel_dx) -> np.ndarray:
        from dfastbe.bank_erosion.utils import moving_avg
        bank_face_indices = single_bank.bank_face_indices
        vel_bank = (
                np.abs(
                    self.velocity_x_face[bank_face_indices] * single_bank.dx
                    + self.velocity_y_face[bank_face_indices] * single_bank.dy
                )
                / single_bank.segment_length
        )

        if vel_dx > 0.0:
            vel_bank = moving_avg(
                single_bank.bank_chainage_midpoints, vel_bank, vel_dx
            )

        return vel_bank

    def calculate_bank_height(self, single_bank: SingleBank, zb_dx):
        # bank height = maximum bed elevation per cell
        from dfastbe.bank_erosion.utils import moving_avg
        bank_index = single_bank.bank_face_indices
        if self.bed_elevation_location == "node":
            zb_nodes = self.bed_elevation_values
            zb_all = self.apply_masked_indexing(
                zb_nodes, self.face_node[bank_index, :]
            )
            zb_bank = zb_all.max(axis=1)
            if zb_dx > 0.0:
                zb_bank = moving_avg(
                    single_bank.bank_chainage_midpoints, zb_bank, zb_dx,
                )
        else:
            # don't know ... need to check neighbouring cells ...
            zb_bank = None

        return zb_bank

    def get_fairway_data(self, fairway_face_indices):
        # get fairway face indices
        fairway_face_indices = fairway_face_indices

        # get water depth along the fair-way
        water_depth_fairway = self.water_depth_face[fairway_face_indices]
        water_level = self.water_level_face[fairway_face_indices]
        chez_face = self.chezy_face[fairway_face_indices]
        chezy = 0 * chez_face + chez_face.mean()

        data = {
            "water_depth": water_depth_fairway,
            "water_level": water_level,
            "chezy": chezy,
        }
        return data

`apply_masked_indexing(x0: np.array, idx: np.ma.masked_array) -> np.ma.masked_array` `staticmethod` #

Index one array by another transferring the mask.

Parameters:

Name	Type	Description	Default
`x0`		np.ndarray A linear array.	required
`idx`		np.ma.masked_array An index array with possibly masked indices.	required

Returns:

Name	Type	Description
`x1`	`masked_array`	np.ma.masked_array An array with same shape as idx, with mask.

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

@staticmethod
def apply_masked_indexing(
    x0: np.array, idx: np.ma.masked_array
) -> np.ma.masked_array:
    """
    Index one array by another transferring the mask.

    Args:
        x0 : np.ndarray
            A linear array.
        idx : np.ma.masked_array
            An index array with possibly masked indices.

    returns:
        x1: np.ma.masked_array
            An array with same shape as idx, with mask.
    """
    idx_safe = idx.copy()
    idx_safe.data[np.ma.getmask(idx)] = 0
    x1 = np.ma.masked_where(np.ma.getmask(idx), x0[idx_safe])
    return x1

`compute_mesh_topology() -> MeshData` #

Derive secondary topology arrays from the face-node connectivity of the mesh.

This function computes the edge-node, edge-face, and face-edge connectivity arrays, as well as the boundary edges of the mesh, based on the face-node connectivity provided in the simulation data.

Returns:

Name	Type	Description
`MeshData`	`MeshData`	a dataclass containing the following attributes: - `x_face_coords`: x-coordinates of face nodes - `y_face_coords`: y-coordinates of face nodes - `x_edge_coords`: x-coordinates of edge nodes - `y_edge_coords`: y-coordinates of edge nodes - `face_node`: the node indices for each of the mesh faces. - `n_nodes`: number of nodes per face - `edge_node`: the node indices for each of the mesh edges. - `edge_face_connectivity`: the face indices for each of the mesh edge - `face_edge_connectivity`: the edge indices for each of the mesh face - `boundary_edge_nrs`: indices of boundary edges

Raises:

Type	Description
`KeyError`	If required keys (e.g., `face_node`, `nnodes`, `x_node`, `y_node`) are missing from the `sim` object.

Notes

The function identifies unique edges by sorting and comparing node indices.
Boundary edges are identified as edges that belong to only one face.
The function assumes that the mesh is well-formed, with consistent face-node connectivity.

Source code in src/dfastbe/bank_erosion/data_models/inputs.py

def compute_mesh_topology(self) -> MeshData:
    """Derive secondary topology arrays from the face-node connectivity of the mesh.

    This function computes the edge-node, edge-face, and face-edge connectivity arrays,
    as well as the boundary edges of the mesh, based on the face-node connectivity provided
    in the simulation data.

    Returns:
        MeshData: a dataclass containing the following attributes:
            - `x_face_coords`: x-coordinates of face nodes
            - `y_face_coords`: y-coordinates of face nodes
            - `x_edge_coords`: x-coordinates of edge nodes
            - `y_edge_coords`: y-coordinates of edge nodes
            - `face_node`: the node indices for each of the mesh faces.
            - `n_nodes`: number of nodes per face
            - `edge_node`: the node indices for each of the mesh edges.
            - `edge_face_connectivity`: the face indices for each of the mesh edge
            - `face_edge_connectivity`: the edge indices for each of the mesh face
            - `boundary_edge_nrs`: indices of boundary edges

    Raises:
        KeyError:
            If required keys (e.g., `face_node`, `nnodes`, `x_node`, `y_node`) are missing from the `sim` object.

    Notes:
        - The function identifies unique edges by sorting and comparing node indices.
        - Boundary edges are identified as edges that belong to only one face.
        - The function assumes that the mesh is well-formed, with consistent face-node connectivity.
    """

    # get a sorted list of edge node connections (shared edges occur twice)
    # face_nr contains the face index to which the edge belongs
    n_faces = self.face_node.shape[0]
    n_edges = sum(self.n_nodes)
    edge_node = np.zeros((n_edges, 2), dtype=int)
    face_nr = np.zeros((n_edges,), dtype=int)
    i = 0
    for face_i in range(n_faces):
        num_edges = self.n_nodes[face_i]  # note: nEdges = nNodes
        for edge_i in range(num_edges):
            if edge_i == 0:
                edge_node[i, 1] = self.face_node[face_i, num_edges - 1]
            else:
                edge_node[i, 1] = self.face_node[face_i, edge_i - 1]
            edge_node[i, 0] = self.face_node[face_i, edge_i]
            face_nr[i] = face_i
            i = i + 1
    edge_node.sort(axis=1)
    i2 = np.argsort(edge_node[:, 1], kind="stable")
    i1 = np.argsort(edge_node[i2, 0], kind="stable")
    i12 = i2[i1]
    edge_node = edge_node[i12, :]
    face_nr = face_nr[i12]

    # detect which edges are equal to the previous edge, and get a list of all unique edges
    numpy_true = np.array([True])
    equal_to_previous = np.concatenate(
        (~numpy_true, (np.diff(edge_node, axis=0) == 0).all(axis=1))
    )
    unique_edge = ~equal_to_previous
    n_unique_edges = np.sum(unique_edge)
    # reduce the edge node connections to only the unique edges
    edge_node = edge_node[unique_edge, :]

    # number the edges
    edge_nr = np.zeros(n_edges, dtype=int)
    edge_nr[unique_edge] = np.arange(n_unique_edges, dtype=int)
    edge_nr[equal_to_previous] = edge_nr[
        np.concatenate((equal_to_previous[1:], equal_to_previous[:1]))
    ]

    # if two consecutive edges are unique, the first one occurs only once and represents a boundary edge
    is_boundary_edge = unique_edge & np.concatenate((unique_edge[1:], numpy_true))
    boundary_edge_nrs = edge_nr[is_boundary_edge]

    # go back to the original face order
    edge_nr_in_face_order = np.zeros(n_edges, dtype=int)
    edge_nr_in_face_order[i12] = edge_nr
    # create the face edge connectivity array
    face_edge_connectivity = np.zeros(self.face_node.shape, dtype=int)

    i = 0
    for face_i in range(n_faces):
        num_edges = self.n_nodes[face_i]  # note: num_edges = n_nodes
        for edge_i in range(num_edges):
            face_edge_connectivity[face_i, edge_i] = edge_nr_in_face_order[i]
            i = i + 1

    # determine the edge face connectivity
    edge_face = -np.ones((n_unique_edges, 2), dtype=int)
    edge_face[edge_nr[unique_edge], 0] = face_nr[unique_edge]
    edge_face[edge_nr[equal_to_previous], 1] = face_nr[equal_to_previous]

    x_face_coords = self.apply_masked_indexing(
        self.x_node, self.face_node
    )
    y_face_coords = self.apply_masked_indexing(
        self.y_node, self.face_node
    )
    x_edge_coords = self.x_node[edge_node]
    y_edge_coords = self.y_node[edge_node]

    return MeshData(
        x_face_coords=x_face_coords,
        y_face_coords=y_face_coords,
        x_edge_coords=x_edge_coords,
        y_edge_coords=y_edge_coords,
        face_node=self.face_node,
        n_nodes=self.n_nodes,
        edge_node=edge_node,
        edge_face_connectivity=edge_face,
        face_edge_connectivity=face_edge_connectivity,
        boundary_edge_nrs=boundary_edge_nrs,
    )

For more details, see Bank Erosion Input Data Models.

Workflow #

The typical workflow for bank erosion calculation is:

Initialize the Erosion object with a configuration file
Call the run method to start the erosion calculation process
The run method orchestrates the entire process:
Processes the river axis
Gets fairway data
Calculates bank-fairway distance
Prepares initial conditions
Processes discharge levels
Computes erosion per level
Post-processes results
Writes output files
Generates plots

Usage Example #

from dfastbe.io.config import ConfigFile
from dfastbe.bank_erosion.bank_erosion import Erosion

# Load configuration file
config_file = ConfigFile.read("config.cfg")

# Initialize Erosion object
erosion = Erosion(config_file)

# Run erosion calculation
erosion.run()

For more details on the specific methods and classes, refer to the API reference below.