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

Copyright (C) 2025 Stichting Deltares.

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

All indications and logos of, and references to, "Delft3D" and "Deltares" are registered trademarks of Stichting Deltares, and remain the property of Stichting Deltares. All rights reserved.

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

Erosion

Bases: BaseCalculator

Class to handle the bank erosion calculations.

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

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

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

        self.river_data = ErosionRiverData(config_file)
        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.debugger = Debugger(config_file.crs, self.river_data.output_dir)
        self.erosion_calculator = ErosionCalculator()

    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,
        num_stations_per_bank: List[int],
        fairway_data: FairwayData,
    ) -> ErosionInputs:
        # wave reduction s0, s1
        parameters = self._get_parameters(num_stations_per_bank)
        wave_fairway_distance_0 = parameters["Wave0"]
        wave_fairway_distance_1 = parameters["Wave1"]

        # save 1_banklines
        # read vship, nship, nwave, draught (tship), shiptype ... independent of level number
        ships_parameters = ShipsParameters.get_ship_data(
            num_stations_per_bank, self.config_file
        )

        # read classes flag (yes: banktype = taucp, no: banktype = tauc) and banktype (taucp: 0-4 ... or ... tauc = critical shear value)
        classes = self.config_file.get_bool("Erosion", "Classes")
        if classes:
            bank_type = parameters["BankType"]
            tauc = []
            for bank in bank_type:
                tauc.append(ErosionInputs.taucls[bank])
        else:
            tauc = parameters["BankType"]
            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

        dike_height_default = -1000
        dike_height = parameters["ProtectionLevel"]
        # 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 == dike_height_default
            one_zss[mask] = fairway_data.fairway_initial_water_levels[ib][mask] - 1

        data = {
            '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=ships_parameters, bank_type=bank_type
        )

    def _get_erosion_input_parameters(self) -> List[Parameters]:
        return [
            Parameters(name="Wave0", default=200, valid=True, onefile=True, positive=None, ext=None),
            Parameters(name="Wave1", default=150, valid=True, onefile=True, positive=None, ext=None),
            Parameters(name="BankType", default=0, valid=None, onefile=None, positive=None, ext=".btp"),
            Parameters(name="ProtectionLevel", default=-1000, valid=None, onefile=None, positive=None, ext=".bpl"),
        ]

    def _get_parameters(self, num_stations_per_bank) -> Dict[str, Any]:
        """Get a parameter from the configuration file."""
        data = {}
        for parameter in self._get_erosion_input_parameters():
            data[parameter.name] = self.config_file.get_parameter(
                "Erosion",
                parameter.name,
                num_stations_per_bank,
                default=parameter.default,
                positive=parameter.positive,
                onefile=parameter.onefile,
                ext=parameter.ext,
            )
        return data

    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,
        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 = erosion_inputs.shipping_data.read_discharge_parameters(
                level_i, 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 = self.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 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 get_mesh_processor(self):
        log_text("derive_topology")
        mesh_data = self.simulation_data.compute_mesh_topology(verbose=False)

        return MeshProcessor(self.river_data, mesh_data, verbose=self.config_file.debug)

    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("-")

        mesh_processor = self.get_mesh_processor()

        river_axis = self.river_data.process_river_axis_by_center_line()
        fairway_data = mesh_processor.get_fairway_data(river_axis)

        # 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

        # map bank lines to mesh cells
        log_text("intersect_bank_mesh")
        bank_data = mesh_processor.get_bank_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(
            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,
            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.results = {
            "river_axis": river_axis,
            "bank_data": bank_data,
            "water_level_data": water_level_data,
            "erosion_results": erosion_results,
            "erosion_inputs": erosion_inputs,
            "bankline_new_list": bankline_new_list,
            "bankline_eq_list": bankline_eq_list,
            "xy_line_eq_list": xy_line_eq_list,
            "km_mid": km_mid,
        }

    def plot(self):
        # create various plots
        if self.river_data.plot_flags.plot_data:
            plotter = ErosionPlotter(
                self.gui,
                self.river_data.plot_flags,
                self.results["erosion_results"],
                self.results["bank_data"],
                self.results["water_level_data"],
                self.results["erosion_inputs"],
            )
            plotter.plot_all(
                self.results["river_axis"].data["stations"],
                self.results["xy_line_eq_list"],
                self.results["km_mid"],
                self.river_data.output_intervals,
                self.river_data.river_center_line.as_array(),
                self.simulation_data,
            )

    def save(self):
        self._write_bankline_shapefiles(
            self.results["bankline_new_list"],
            self.results["bankline_eq_list"],
            self.config_file,
        )
        self._write_volume_outputs(
            self.results["erosion_results"], self.results["km_mid"]
        )

    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),
        )

__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."""
    super().__init__(config_file, gui)

    self.river_data = ErosionRiverData(config_file)
    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.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	Description	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

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("-")

    mesh_processor = self.get_mesh_processor()

    river_axis = self.river_data.process_river_axis_by_center_line()
    fairway_data = mesh_processor.get_fairway_data(river_axis)

    # 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

    # map bank lines to mesh cells
    log_text("intersect_bank_mesh")
    bank_data = mesh_processor.get_bank_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(
        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,
        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.results = {
        "river_axis": river_axis,
        "bank_data": bank_data,
        "water_level_data": water_level_data,
        "erosion_results": erosion_results,
        "erosion_inputs": erosion_inputs,
        "bankline_new_list": bankline_new_list,
        "bankline_eq_list": bankline_eq_list,
        "xy_line_eq_list": xy_line_eq_list,
        "km_mid": km_mid,
    }

calculate_alpha(coords: np.ndarray, ind_1: int, ind_2: int, bp: Tuple[int, Any])

Calculate the alpha value for the bank erosion model.

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

def calculate_alpha(coords: np.ndarray, ind_1: int, ind_2: int, bp: Tuple[int, Any]):
    """Calculate the alpha value for the bank erosion model."""
    alpha = (
        (coords[ind_1, 0] - coords[ind_2, 0]) * (bp[0] - coords[ind_2, 0])
        + (coords[ind_1, 1] - coords[ind_2, 1]) * (bp[1] - coords[ind_2, 1])
    ) / (
        (coords[ind_1, 0] - coords[ind_2, 0]) ** 2
        + (coords[ind_1, 1] - coords[ind_2, 1]) ** 2
    )

    return alpha

Mesh data models

dfastbe.bank_erosion.mesh.data_models

This module defines data structures and methods for handling mesh data and river segments.

Edges dataclass

Dataclass to hold edge candidates for left and right edges.

Parameters:

Name	Type	Description	Default
left	int	Index of the left edge.	required
left_theta	float	Angle of the left edge in radians.	required
right	int	Index of the right edge.	required
right_theta	float	Angle of the right edge in radians.	required
found	bool	Flag indicating whether a valid edge pair was found.	False

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

@dataclass
class Edges:
    """Dataclass to hold edge candidates for left and right edges.

    Args:
        left (int):
            Index of the left edge.
        left_theta (float):
            Angle of the left edge in radians.
        right (int):
            Index of the right edge.
        right_theta (float):
            Angle of the right edge in radians.
        found (bool):
            Flag indicating whether a valid edge pair was found.
    """

    left: int
    left_theta: float
    right: int
    right_theta: float
    found: bool = False

    def update_edges_by_angle(
        self, edge_index: int, dtheta: float, j, verbose: bool = False
    ):
        """Update the left and right edges based on the angle difference."""

        if dtheta > 0:
            if dtheta < self.left_theta:
                self.left = edge_index
                self.left_theta = dtheta
            if TWO_PI - dtheta < self.right_theta:
                self.right = edge_index
                self.right_theta = TWO_PI - dtheta
        elif dtheta < 0:
            dtheta = -dtheta
            if TWO_PI - dtheta < self.left_theta:
                self.left = edge_index
                self.left_theta = TWO_PI - dtheta
            if dtheta < self.right_theta:
                self.right = edge_index
                self.right_theta = dtheta
        else:
            # aligned with edge
            if verbose and j is not None:
                print(f"{j}: line is aligned with edge {edge_index}")

            self.left = edge_index
            self.right = edge_index
            self.found = True

update_edges_by_angle(edge_index: int, dtheta: float, j, verbose: bool = False)

Update the left and right edges based on the angle difference.

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

def update_edges_by_angle(
    self, edge_index: int, dtheta: float, j, verbose: bool = False
):
    """Update the left and right edges based on the angle difference."""

    if dtheta > 0:
        if dtheta < self.left_theta:
            self.left = edge_index
            self.left_theta = dtheta
        if TWO_PI - dtheta < self.right_theta:
            self.right = edge_index
            self.right_theta = TWO_PI - dtheta
    elif dtheta < 0:
        dtheta = -dtheta
        if TWO_PI - dtheta < self.left_theta:
            self.left = edge_index
            self.left_theta = TWO_PI - dtheta
        if dtheta < self.right_theta:
            self.right = edge_index
            self.right_theta = dtheta
    else:
        # aligned with edge
        if verbose and j is not None:
            print(f"{j}: line is aligned with edge {edge_index}")

        self.left = edge_index
        self.right = edge_index
        self.found = True

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/mesh/data_models.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
    verbose: bool = False

    def get_face_by_index(self, index: int, as_polygon: bool = False) -> np.ndarray | Polygon:
        """Returns the coordinates of the index-th mesh face as an (N, 2) array.

        Args:
            index:
                The face index.
            as_polygon:
                whither to return the face as a shapely polygon or not. Default is False

        Returns:
            np.ndarray:
                Array of shape (n_nodes, 2) with x, y coordinates.
        """
        x = self.x_face_coords[index : index + 1, : self.n_nodes[index]]
        y = self.y_face_coords[index : index + 1, : self.n_nodes[index]]
        face = np.concatenate((x, y), axis=0).T

        if as_polygon:
            face = Polygon(face)
        return face

    def locate_point(
        self, point: Point | np.ndarray | list | Tuple, face_index: int | np.ndarray
    ) -> int | List[int]:
        """Locate a point in the mesh faces.

        Args:
            point:
                The point to check.
            face_index:
                The index of the mesh face.

        Returns:
            indexes (int|list[int]):
                index if the face that the point is located in, or a list of indexes if the point is on the edge of
                multiple faces.
        """
        if not isinstance(point, Point) and isinstance(
            point, (list, tuple, np.ndarray)
        ):
            point = Point(point)
        else:
            raise TypeError(
                "point must be a Point object, a list, a tuple, or an np.ndarray of coordinates"
            )

        index_list = []
        for ind in face_index:
            face = self.get_face_by_index(ind)
            face_polygon = Polygon(face)
            if face_polygon.contains(point):
                return ind
            else:
                # create a closed line string from the face coordinates
                closed_line = LineString(np.vstack([face, face[0]]))
                if closed_line.contains(point):
                    index_list.append(ind)

        return index_list

    def find_segment_intersections(
        self,
        index: int,
        segment: RiverSegment,
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Calculate the intersection of a line segment with the edges of a mesh face.

        This function determines where a line segment (defined by two points) intersects the edges of a mesh face.
        It returns the relative distances along the segment and the edges where the intersections occur, as well as
        flags indicating whether the intersections occur at nodes.

        Args:
            index (int):
                Index of the current mesh face. If `index` is negative, the function assumes the segment intersects
                the boundary edges of the mesh.
            segment (RiverSegment):
                A `RiverSegment` object containing the previous and current points of the segment, as well as the
                minimum relative distance along the segment where the last intersection occurred.

        Returns:
            Tuple[np.ndarray, np.ndarray, np.ndarray]:
                - b (np.ndarray):
                    Relative distances along the segment `bpj1-bpj` where the intersections occur.
                - edges (np.ndarray):
                    Indices of the edges that are intersected by the segment.
                - nodes (np.ndarray):
                    Flags indicating whether the intersections occur at nodes. A value of `-1` indicates no
                    intersection at a node, while other values correspond to node indices.

        Raises:
            ValueError:
                If the input data is invalid or inconsistent.

        Notes:
            - If `index` is negative, the function assumes the segment intersects the boundary edges of the mesh.
            - The function uses the `get_slices_core` helper function to calculate the intersections.
            - Intersections at nodes are flagged in the `nodes` array, with the corresponding node indices.

        """
        if index < 0:
            edges = self.boundary_edge_nrs
        else:
            edges = self.face_edge_connectivity[index, : self.n_nodes[index]]
        edge_relative_dist, segment_relative_dist, edges = (
            self.calculate_edge_intersections(edges, segment, True)
        )
        is_intersected_at_node = -np.ones(edge_relative_dist.shape, dtype=np.int64)
        is_intersected_at_node[edge_relative_dist == 0] = self.edge_node[
            edges[edge_relative_dist == 0], 0
        ]
        is_intersected_at_node[edge_relative_dist == 1] = self.edge_node[
            edges[edge_relative_dist == 1], 1
        ]

        return segment_relative_dist, edges, is_intersected_at_node

    def calculate_edge_intersections(
        self,
        edges: np.ndarray,
        segment: RiverSegment,
        limit_relative_distance: bool = True,
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Calculate the intersection of a line segment with multiple mesh edges.

        This function determines where a line segment intersects a set of mesh edges.
        It calculates the relative distances along the segment and the edges where
        the intersections occur, and returns the indices of the intersected edges.

        Args:
            edges (np.ndarray):
                Array containing the indices of the edges to check for intersections.
            segment (RiverSegment):
                A `RiverSegment` object containing the previous and current points of the segment, as well as the
                minimum relative distance along the segment where the last intersection occurred.
            limit_relative_distance (bool, optional):
                If True, limits the relative distance along the segment `bpj1-bpj`
                to a maximum of 1. Defaults to True.

        Returns:
            Tuple[np.ndarray, np.ndarray, np.ndarray]:
                - a (np.ndarray): Relative distances along the edges where the
                  intersections occur.
                - b (np.ndarray): Relative distances along the segment `bpj1-bpj`
                  where the intersections occur.
                - edges (np.ndarray): Indices of the edges that are intersected
                  by the segment.

        Raises:
            ValueError:
                If the input data is invalid or inconsistent.

        Notes:
            - The function uses the `get_slices_ab` helper function to calculate the
              relative distances `a` and `b` for each edge.
            - The `bmin` parameter is used to filter out intersections that occur
              too close to the starting point of the segment.
            - If `limit_relative_distance` is True, intersections beyond the endpoint of the segment
              are ignored.
        """
        from dfastbe.bank_erosion.utils import calculate_segment_edge_intersections

        edge_relative_dist, segment_relative_dist, valid_intersections = (
            calculate_segment_edge_intersections(
                self.x_edge_coords[edges, 0],
                self.y_edge_coords[edges, 0],
                self.x_edge_coords[edges, 1],
                self.y_edge_coords[edges, 1],
                segment.previous_point[0],
                segment.previous_point[1],
                segment.current_point[0],
                segment.current_point[1],
                segment.min_relative_distance,
                limit_relative_distance,
            )
        )
        edges = edges[valid_intersections]
        return edge_relative_dist, segment_relative_dist, edges

    def calculate_edge_angle(self, edge: int, reverse: bool = False) -> float:
        """Calculate the angle of a mesh edge in radians.

        Args:
            edge (int):
                The edge index.
            reverse (bool):
                If True, computes the angle from end to start.

        Returns:
            float: The angle of the edge in radians.
        """
        start, end = (1, 0) if reverse else (0, 1)
        dx = self.x_edge_coords[edge, end] - self.x_edge_coords[edge, start]
        dy = self.y_edge_coords[edge, end] - self.y_edge_coords[edge, start]

        return math.atan2(dy, dx)

    def find_edges(self, theta, node, verbose_index: int = None) -> Edges:
        """
        Helper to find the left and right edges at a node based on the direction theta.

        Args:
            theta (float):
                Direction angle of the segment.
            node (int):
                The node index.
            verbose_index (int, optional):
                Step index for verbose output.

        Returns:
            Edges:
                A dataclass containing the left and right edge indices, their angle differences, and a found flag.
        """
        all_node_edges = np.nonzero((self.edge_node == node).any(axis=1))[0]

        if self.verbose and verbose_index is not None:
            print(
                f"{verbose_index}: the edges connected to node {node} are {all_node_edges}"
            )

        edges = Edges(left=-1, left_theta=TWO_PI, right=-1, right_theta=TWO_PI)

        for ie in all_node_edges:
            reverse = self.edge_node[ie, 0] != node
            theta_edge = self.calculate_edge_angle(ie, reverse=reverse)

            if self.verbose and verbose_index is not None:
                print(f"{verbose_index}: edge {ie} connects {self.edge_node[ie, :]}")
                print(f"{verbose_index}: edge {ie} theta is {theta_edge}")

            dtheta = theta_edge - theta

            edges.update_edges_by_angle(ie, dtheta, verbose_index)
            if edges.found:
                break

        if self.verbose and verbose_index is not None:
            print(f"{verbose_index}: the edge to the left is edge {edges.left}")
            print(f"{verbose_index}: the edge to the right is edge {edges.right}")

        return edges

    def resolve_next_face_from_edges(
        self, node, edges: Edges, verbose_index: int = None
    ) -> int:
        """
        Helper to resolve the next face index when traversing between two edges at a node.

        Args:
            node (int): The node index.
            edges (Edges):
                The edges connecting the node, containing left and right edge indices.
            verbose_index (int, optional):
                Step index for verbose output.

        Returns:
            next_face_index (int):
                The next face index.
        """
        left_faces = self.edge_face_connectivity[edges.left, :]
        right_faces = self.edge_face_connectivity[edges.right, :]

        if left_faces[0] in right_faces and left_faces[1] in right_faces:
            fn1 = self.face_node[left_faces[0]]
            fe1 = self.face_edge_connectivity[left_faces[0]]

            if self.verbose and verbose_index is not None:
                print(
                    f"{verbose_index}: those edges are shared by two faces: {left_faces}"
                )
                print(f"{verbose_index}: face {left_faces[0]} has nodes: {fn1}")
                print(f"{verbose_index}: face {left_faces[0]} has edges: {fe1}")

            # nodes of the face should be listed in clockwise order
            # 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] == edges.right:
                next_face_index = left_faces[0]
            else:
                next_face_index = left_faces[1]

        elif left_faces[0] in right_faces:
            next_face_index = left_faces[0]
        elif left_faces[1] in right_faces:
            next_face_index = left_faces[1]
        else:
            raise ValueError(
                f"Shouldn't come here .... left edge {edges.left}"
                f" and right edge {edges.right} don't share any face"
            )
        return next_face_index

calculate_edge_angle(edge: int, reverse: bool = False) -> float

Calculate the angle of a mesh edge in radians.

Parameters:

Name	Type	Description	Default
edge	int	The edge index.	required
reverse	bool	If True, computes the angle from end to start.	False

Returns:

Name	Type	Description
float	float	The angle of the edge in radians.

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

def calculate_edge_angle(self, edge: int, reverse: bool = False) -> float:
    """Calculate the angle of a mesh edge in radians.

    Args:
        edge (int):
            The edge index.
        reverse (bool):
            If True, computes the angle from end to start.

    Returns:
        float: The angle of the edge in radians.
    """
    start, end = (1, 0) if reverse else (0, 1)
    dx = self.x_edge_coords[edge, end] - self.x_edge_coords[edge, start]
    dy = self.y_edge_coords[edge, end] - self.y_edge_coords[edge, start]

    return math.atan2(dy, dx)

calculate_edge_intersections(edges: np.ndarray, segment: RiverSegment, limit_relative_distance: bool = True) -> Tuple[np.ndarray, np.ndarray, np.ndarray]

Calculate the intersection of a line segment with multiple mesh edges.

This function determines where a line segment intersects a set of mesh edges. It calculates the relative distances along the segment and the edges where the intersections occur, and returns the indices of the intersected edges.

Parameters:

Name	Type	Description	Default
edges	ndarray	Array containing the indices of the edges to check for intersections.	required
segment	RiverSegment	A RiverSegment object containing the previous and current points of the segment, as well as the minimum relative distance along the segment where the last intersection occurred.	required
limit_relative_distance	bool	If True, limits the relative distance along the segment bpj1-bpj to a maximum of 1. Defaults to True.	True

Returns:

Type	Description
Tuple[ndarray, ndarray, ndarray]	Tuple[np.ndarray, np.ndarray, np.ndarray]: - a (np.ndarray): Relative distances along the edges where the intersections occur. - b (np.ndarray): Relative distances along the segment bpj1-bpj where the intersections occur. - edges (np.ndarray): Indices of the edges that are intersected by the segment.

Raises:

Type	Description
ValueError	If the input data is invalid or inconsistent.

Notes

• The function uses the get_slices_ab helper function to calculate the relative distances a and b for each edge.
• The bmin parameter is used to filter out intersections that occur too close to the starting point of the segment.
• If limit_relative_distance is True, intersections beyond the endpoint of the segment are ignored.

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

def calculate_edge_intersections(
    self,
    edges: np.ndarray,
    segment: RiverSegment,
    limit_relative_distance: bool = True,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Calculate the intersection of a line segment with multiple mesh edges.

    This function determines where a line segment intersects a set of mesh edges.
    It calculates the relative distances along the segment and the edges where
    the intersections occur, and returns the indices of the intersected edges.

    Args:
        edges (np.ndarray):
            Array containing the indices of the edges to check for intersections.
        segment (RiverSegment):
            A `RiverSegment` object containing the previous and current points of the segment, as well as the
            minimum relative distance along the segment where the last intersection occurred.
        limit_relative_distance (bool, optional):
            If True, limits the relative distance along the segment `bpj1-bpj`
            to a maximum of 1. Defaults to True.

    Returns:
        Tuple[np.ndarray, np.ndarray, np.ndarray]:
            - a (np.ndarray): Relative distances along the edges where the
              intersections occur.
            - b (np.ndarray): Relative distances along the segment `bpj1-bpj`
              where the intersections occur.
            - edges (np.ndarray): Indices of the edges that are intersected
              by the segment.

    Raises:
        ValueError:
            If the input data is invalid or inconsistent.

    Notes:
        - The function uses the `get_slices_ab` helper function to calculate the
          relative distances `a` and `b` for each edge.
        - The `bmin` parameter is used to filter out intersections that occur
          too close to the starting point of the segment.
        - If `limit_relative_distance` is True, intersections beyond the endpoint of the segment
          are ignored.
    """
    from dfastbe.bank_erosion.utils import calculate_segment_edge_intersections

    edge_relative_dist, segment_relative_dist, valid_intersections = (
        calculate_segment_edge_intersections(
            self.x_edge_coords[edges, 0],
            self.y_edge_coords[edges, 0],
            self.x_edge_coords[edges, 1],
            self.y_edge_coords[edges, 1],
            segment.previous_point[0],
            segment.previous_point[1],
            segment.current_point[0],
            segment.current_point[1],
            segment.min_relative_distance,
            limit_relative_distance,
        )
    )
    edges = edges[valid_intersections]
    return edge_relative_dist, segment_relative_dist, edges

find_edges(theta, node, verbose_index: int = None) -> Edges

Helper to find the left and right edges at a node based on the direction theta.

Parameters:

Name	Type	Description	Default
theta	float	Direction angle of the segment.	required
node	int	The node index.	required
verbose_index	int	Step index for verbose output.	None

Returns:

Name	Type	Description
Edges	Edges	A dataclass containing the left and right edge indices, their angle differences, and a found flag.

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

def find_edges(self, theta, node, verbose_index: int = None) -> Edges:
    """
    Helper to find the left and right edges at a node based on the direction theta.

    Args:
        theta (float):
            Direction angle of the segment.
        node (int):
            The node index.
        verbose_index (int, optional):
            Step index for verbose output.

    Returns:
        Edges:
            A dataclass containing the left and right edge indices, their angle differences, and a found flag.
    """
    all_node_edges = np.nonzero((self.edge_node == node).any(axis=1))[0]

    if self.verbose and verbose_index is not None:
        print(
            f"{verbose_index}: the edges connected to node {node} are {all_node_edges}"
        )

    edges = Edges(left=-1, left_theta=TWO_PI, right=-1, right_theta=TWO_PI)

    for ie in all_node_edges:
        reverse = self.edge_node[ie, 0] != node
        theta_edge = self.calculate_edge_angle(ie, reverse=reverse)

        if self.verbose and verbose_index is not None:
            print(f"{verbose_index}: edge {ie} connects {self.edge_node[ie, :]}")
            print(f"{verbose_index}: edge {ie} theta is {theta_edge}")

        dtheta = theta_edge - theta

        edges.update_edges_by_angle(ie, dtheta, verbose_index)
        if edges.found:
            break

    if self.verbose and verbose_index is not None:
        print(f"{verbose_index}: the edge to the left is edge {edges.left}")
        print(f"{verbose_index}: the edge to the right is edge {edges.right}")

    return edges

find_segment_intersections(index: int, segment: RiverSegment) -> Tuple[np.ndarray, np.ndarray, np.ndarray]

Calculate the intersection of a line segment with the edges of a mesh face.

This function determines where a line segment (defined by two points) intersects the edges of a mesh face. It returns the relative distances along the segment and the edges where the intersections occur, as well as flags indicating whether the intersections occur at nodes.

Parameters:

Name	Type	Description	Default
index	int	Index of the current mesh face. If index is negative, the function assumes the segment intersects the boundary edges of the mesh.	required
segment	RiverSegment	A RiverSegment object containing the previous and current points of the segment, as well as the minimum relative distance along the segment where the last intersection occurred.	required

Returns:

Type	Description
Tuple[ndarray, ndarray, ndarray]	Tuple[np.ndarray, np.ndarray, np.ndarray]: - b (np.ndarray): Relative distances along the segment bpj1-bpj where the intersections occur. - edges (np.ndarray): Indices of the edges that are intersected by the segment. - nodes (np.ndarray): Flags indicating whether the intersections occur at nodes. A value of -1 indicates no intersection at a node, while other values correspond to node indices.

Raises:

Type	Description
ValueError	If the input data is invalid or inconsistent.

Notes

• If index is negative, the function assumes the segment intersects the boundary edges of the mesh.
• The function uses the get_slices_core helper function to calculate the intersections.
• Intersections at nodes are flagged in the nodes array, with the corresponding node indices.

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

def find_segment_intersections(
    self,
    index: int,
    segment: RiverSegment,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Calculate the intersection of a line segment with the edges of a mesh face.

    This function determines where a line segment (defined by two points) intersects the edges of a mesh face.
    It returns the relative distances along the segment and the edges where the intersections occur, as well as
    flags indicating whether the intersections occur at nodes.

    Args:
        index (int):
            Index of the current mesh face. If `index` is negative, the function assumes the segment intersects
            the boundary edges of the mesh.
        segment (RiverSegment):
            A `RiverSegment` object containing the previous and current points of the segment, as well as the
            minimum relative distance along the segment where the last intersection occurred.

    Returns:
        Tuple[np.ndarray, np.ndarray, np.ndarray]:
            - b (np.ndarray):
                Relative distances along the segment `bpj1-bpj` where the intersections occur.
            - edges (np.ndarray):
                Indices of the edges that are intersected by the segment.
            - nodes (np.ndarray):
                Flags indicating whether the intersections occur at nodes. A value of `-1` indicates no
                intersection at a node, while other values correspond to node indices.

    Raises:
        ValueError:
            If the input data is invalid or inconsistent.

    Notes:
        - If `index` is negative, the function assumes the segment intersects the boundary edges of the mesh.
        - The function uses the `get_slices_core` helper function to calculate the intersections.
        - Intersections at nodes are flagged in the `nodes` array, with the corresponding node indices.

    """
    if index < 0:
        edges = self.boundary_edge_nrs
    else:
        edges = self.face_edge_connectivity[index, : self.n_nodes[index]]
    edge_relative_dist, segment_relative_dist, edges = (
        self.calculate_edge_intersections(edges, segment, True)
    )
    is_intersected_at_node = -np.ones(edge_relative_dist.shape, dtype=np.int64)
    is_intersected_at_node[edge_relative_dist == 0] = self.edge_node[
        edges[edge_relative_dist == 0], 0
    ]
    is_intersected_at_node[edge_relative_dist == 1] = self.edge_node[
        edges[edge_relative_dist == 1], 1
    ]

    return segment_relative_dist, edges, is_intersected_at_node

get_face_by_index(index: int, as_polygon: bool = False) -> np.ndarray | Polygon

Returns the coordinates of the index-th mesh face as an (N, 2) array.

Parameters:

Name	Type	Description	Default
index	int	The face index.	required
as_polygon	bool	whither to return the face as a shapely polygon or not. Default is False	False

Returns:

Type	Description
ndarray | Polygon	np.ndarray: Array of shape (n_nodes, 2) with x, y coordinates.

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

def get_face_by_index(self, index: int, as_polygon: bool = False) -> np.ndarray | Polygon:
    """Returns the coordinates of the index-th mesh face as an (N, 2) array.

    Args:
        index:
            The face index.
        as_polygon:
            whither to return the face as a shapely polygon or not. Default is False

    Returns:
        np.ndarray:
            Array of shape (n_nodes, 2) with x, y coordinates.
    """
    x = self.x_face_coords[index : index + 1, : self.n_nodes[index]]
    y = self.y_face_coords[index : index + 1, : self.n_nodes[index]]
    face = np.concatenate((x, y), axis=0).T

    if as_polygon:
        face = Polygon(face)
    return face

locate_point(point: Point | np.ndarray | list | Tuple, face_index: int | np.ndarray) -> int | List[int]

Locate a point in the mesh faces.

Parameters:

Name	Type	Description	Default
point	Point | ndarray | list | Tuple	The point to check.	required
face_index	int | ndarray	The index of the mesh face.	required

Returns:

Name	Type	Description
indexes	int | list[int]	index if the face that the point is located in, or a list of indexes if the point is on the edge of multiple faces.

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

def locate_point(
    self, point: Point | np.ndarray | list | Tuple, face_index: int | np.ndarray
) -> int | List[int]:
    """Locate a point in the mesh faces.

    Args:
        point:
            The point to check.
        face_index:
            The index of the mesh face.

    Returns:
        indexes (int|list[int]):
            index if the face that the point is located in, or a list of indexes if the point is on the edge of
            multiple faces.
    """
    if not isinstance(point, Point) and isinstance(
        point, (list, tuple, np.ndarray)
    ):
        point = Point(point)
    else:
        raise TypeError(
            "point must be a Point object, a list, a tuple, or an np.ndarray of coordinates"
        )

    index_list = []
    for ind in face_index:
        face = self.get_face_by_index(ind)
        face_polygon = Polygon(face)
        if face_polygon.contains(point):
            return ind
        else:
            # create a closed line string from the face coordinates
            closed_line = LineString(np.vstack([face, face[0]]))
            if closed_line.contains(point):
                index_list.append(ind)

    return index_list

resolve_next_face_from_edges(node, edges: Edges, verbose_index: int = None) -> int

Helper to resolve the next face index when traversing between two edges at a node.

Parameters:

Name	Type	Description	Default
node	int	The node index.	required
edges	Edges	The edges connecting the node, containing left and right edge indices.	required
verbose_index	int	Step index for verbose output.	None

Returns:

Name	Type	Description
next_face_index	int	The next face index.

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

def resolve_next_face_from_edges(
    self, node, edges: Edges, verbose_index: int = None
) -> int:
    """
    Helper to resolve the next face index when traversing between two edges at a node.

    Args:
        node (int): The node index.
        edges (Edges):
            The edges connecting the node, containing left and right edge indices.
        verbose_index (int, optional):
            Step index for verbose output.

    Returns:
        next_face_index (int):
            The next face index.
    """
    left_faces = self.edge_face_connectivity[edges.left, :]
    right_faces = self.edge_face_connectivity[edges.right, :]

    if left_faces[0] in right_faces and left_faces[1] in right_faces:
        fn1 = self.face_node[left_faces[0]]
        fe1 = self.face_edge_connectivity[left_faces[0]]

        if self.verbose and verbose_index is not None:
            print(
                f"{verbose_index}: those edges are shared by two faces: {left_faces}"
            )
            print(f"{verbose_index}: face {left_faces[0]} has nodes: {fn1}")
            print(f"{verbose_index}: face {left_faces[0]} has edges: {fe1}")

        # nodes of the face should be listed in clockwise order
        # 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] == edges.right:
            next_face_index = left_faces[0]
        else:
            next_face_index = left_faces[1]

    elif left_faces[0] in right_faces:
        next_face_index = left_faces[0]
    elif left_faces[1] in right_faces:
        next_face_index = left_faces[1]
    else:
        raise ValueError(
            f"Shouldn't come here .... left edge {edges.left}"
            f" and right edge {edges.right} don't share any face"
        )
    return next_face_index

RiverSegment dataclass

Represents a segment of a river line.

Parameters:

Name	Type	Description	Default
min_relative_distance	float	The relative distance along the previous segment where the last intersection occurred. Used to filter intersections along the current segment.	required
current_point	ndarray	A 1D array of shape (2,) containing the x, y coordinates of the current point of the line segment.	required
previous_point	ndarray	A 1D array of shape (2,) containing the x, y coordinates of the previous point of the line segment.	required

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

@dataclass
class RiverSegment:
    """Represents a segment of a river line.

    Args:
        min_relative_distance (float):
            The relative distance along the previous segment where the last intersection occurred. Used to filter
            intersections along the current segment.
        current_point (np.ndarray):
            A 1D array of shape (2,) containing the x, y coordinates of the current point of the line segment.
        previous_point (np.ndarray):
            A 1D array of shape (2,) containing the x, y coordinates of the previous point of the line segment.
    """

    index: int
    min_relative_distance: float
    current_point: np.ndarray
    previous_point: np.ndarray
    distances: np.ndarray = field(default_factory=lambda: np.zeros(0))
    edges: np.ndarray = field(default_factory=lambda: np.zeros(0, dtype=np.int64))
    nodes: np.ndarray = field(default_factory=lambda: np.zeros(0, dtype=np.int64))

    def is_length_zero(self) -> bool:
        """Check if the segment has zero length."""
        return np.array_equal(self.current_point, self.previous_point)

    @property
    def theta(self):
        """Calculate the angle of the segment in radians."""
        theta = math.atan2(
            self.current_point[1] - self.previous_point[1],
            self.current_point[0] - self.previous_point[0],
        )
        return theta

    def select_first_intersection(self):
        """Select the first crossing from a set of edges and their associated distances.

        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
        """
        # get the minimum distance from the distances array
        min_distance_indices = self.distances == np.amin(self.distances)
        self.distances = self.distances[min_distance_indices]
        self.edges = self.edges[min_distance_indices]
        self.nodes = self.nodes[min_distance_indices]

theta property

Calculate the angle of the segment in radians.

is_length_zero() -> bool

Check if the segment has zero length.

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

def is_length_zero(self) -> bool:
    """Check if the segment has zero length."""
    return np.array_equal(self.current_point, self.previous_point)

select_first_intersection()

Select the first crossing from a set of edges and their associated distances.

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

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

def select_first_intersection(self):
    """Select the first crossing from a set of edges and their associated distances.

    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
    """
    # get the minimum distance from the distances array
    min_distance_indices = self.distances == np.amin(self.distances)
    self.distances = self.distances[min_distance_indices]
    self.edges = self.edges[min_distance_indices]
    self.nodes = self.nodes[min_distance_indices]

Mesh Processing

dfastbe.bank_erosion.mesh.processor

module for processing mesh-related operations.

IntersectionState dataclass

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

@dataclass
class IntersectionState:
    coords = np.ndarray
    face_indexes = np.ndarray
    point_index = 0
    current_face_index: int
    vertex_index: np.ndarray | int

    def __init__(self, shape: Tuple[int, int] = (100, 2), verbose: bool = False):
        """Initialize the intersection results with given shape."""
        self.coords = np.zeros(shape, dtype=np.float64)
        self.face_indexes = np.zeros(shape[0], dtype=np.int64)
        self.point_index = 0
        self.current_face_index: int = -1
        self.verbose = verbose

    def update(self, current_bank_point, shape_length: bool = None):
        """Finalize a segment

        Enlarge arrays if needed, set coordinates and index, and increment ind.
        """
        shape_length = self.point_index + 1 if shape_length is None else shape_length

        if self.point_index == self.coords.shape[0]:
            # last coordinate reached, so enlarge arrays
            self.coords = enlarge(self.coords, (shape_length, 2))
            self.face_indexes = enlarge(self.face_indexes, (shape_length,))

        self.coords[self.point_index] = current_bank_point

        if self.current_face_index == -2:
            self.face_indexes[self.point_index] = self.vertex_index[0]
        else:
            self.face_indexes[self.point_index] = self.current_face_index
        self.point_index += 1

    def _log_mesh_transition(self, log_status: Status):
        """Helper to print mesh transition information for debugging."""
        index_str = "outside" if self.current_face_index == -1 else self.current_face_index
        if self.current_face_index == -2:
            index_str = f"edge between {self.vertex_index}"

        log_status.print(index_str)

    def _update_main_attributes(self, log_status):
        if self.verbose:
            log_status.transition_type = "node"
            self._log_mesh_transition(log_status)

        face_indexes = log_status.face_index
        if isinstance(face_indexes, (int, np.integer)):
            self.current_face_index = face_indexes
        elif hasattr(face_indexes, "__len__") and len(face_indexes) == 1:
            self.current_face_index = face_indexes[0]
        else:
            self.current_face_index = -2
            self.vertex_index = face_indexes

    def update_index_and_log(self, status: Status, edge, faces):
        """
        Helper to update mesh index and log transitions for intersect_line_mesh.
        """
        if status.face_index is not None:
            self._update_main_attributes(status)
            return

        status.transition_type = "edge"
        status.transition_index = edge

        for i, face in enumerate(faces):
            if face == self.current_face_index:
                other_face = faces[1 - i]
                if self.verbose:
                    status.face_index = other_face
                    self._log_mesh_transition(status)
                self.current_face_index = other_face
                return

        raise ValueError(
            f"Shouldn't come here .... index {self.current_face_index} differs from both faces "
            f"{faces[0]} and {faces[1]} associated with slicing edge {edge}"
        )

__init__(shape: Tuple[int, int] = (100, 2), verbose: bool = False)

Initialize the intersection results with given shape.

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

def __init__(self, shape: Tuple[int, int] = (100, 2), verbose: bool = False):
    """Initialize the intersection results with given shape."""
    self.coords = np.zeros(shape, dtype=np.float64)
    self.face_indexes = np.zeros(shape[0], dtype=np.int64)
    self.point_index = 0
    self.current_face_index: int = -1
    self.verbose = verbose

update(current_bank_point, shape_length: bool = None)

Finalize a segment

Enlarge arrays if needed, set coordinates and index, and increment ind.

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

def update(self, current_bank_point, shape_length: bool = None):
    """Finalize a segment

    Enlarge arrays if needed, set coordinates and index, and increment ind.
    """
    shape_length = self.point_index + 1 if shape_length is None else shape_length

    if self.point_index == self.coords.shape[0]:
        # last coordinate reached, so enlarge arrays
        self.coords = enlarge(self.coords, (shape_length, 2))
        self.face_indexes = enlarge(self.face_indexes, (shape_length,))

    self.coords[self.point_index] = current_bank_point

    if self.current_face_index == -2:
        self.face_indexes[self.point_index] = self.vertex_index[0]
    else:
        self.face_indexes[self.point_index] = self.current_face_index
    self.point_index += 1

update_index_and_log(status: Status, edge, faces)

Helper to update mesh index and log transitions for intersect_line_mesh.

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

def update_index_and_log(self, status: Status, edge, faces):
    """
    Helper to update mesh index and log transitions for intersect_line_mesh.
    """
    if status.face_index is not None:
        self._update_main_attributes(status)
        return

    status.transition_type = "edge"
    status.transition_index = edge

    for i, face in enumerate(faces):
        if face == self.current_face_index:
            other_face = faces[1 - i]
            if self.verbose:
                status.face_index = other_face
                self._log_mesh_transition(status)
            self.current_face_index = other_face
            return

    raise ValueError(
        f"Shouldn't come here .... index {self.current_face_index} differs from both faces "
        f"{faces[0]} and {faces[1]} associated with slicing edge {edge}"
    )

MeshProcessor

Class to process bank lines and intersect them with a mesh.

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

class MeshProcessor:
    """Class to process bank lines and intersect them with a mesh."""

    def __init__(self, river_data: ErosionRiverData, mesh_data: MeshData, verbose: bool = False):
        """Constructor for MeshProcessor."""
        self.bank_lines = river_data.bank_lines
        self.mesh_data = mesh_data
        self.river_data = river_data
        self.wrapper = MeshWrapper(mesh_data, verbose=verbose)

    def get_fairway_data(self, river_axis: LineGeometry) -> FairwayData:
        log_text("chainage_to_fairway")
        # intersect fairway and mesh
        fairway_intersection_coords, fairway_face_indices = self.wrapper.intersect_with_coords(river_axis.as_array())

        if self.river_data.debug:
            arr = (fairway_intersection_coords[:-1] + fairway_intersection_coords[1:]) / 2
            line_geom = LineGeometry(arr, crs=river_axis.crs)
            line_geom.to_file(
                file_name=f"{str(self.river_data.output_dir)}/fairway_face_indices.shp",
                data={"iface": fairway_face_indices},
            )

        return FairwayData(fairway_face_indices, fairway_intersection_coords)

    def get_bank_data(self) -> BankData:
        """Intersect bank lines with a mesh and return bank data.

        Returns:
            BankData object containing bank line coordinates, face indices, and other bank-related data.
        """
        n_bank_lines = len(self.bank_lines)

        bank_line_coords = []
        bank_face_indices = []
        for bank_index in range(n_bank_lines):
            line_coords = np.array(self.bank_lines.geometry[bank_index].coords)
            log_text("bank_nodes", data={"ib": bank_index + 1, "n": len(line_coords)})

            coords_along_bank, face_indices = self.wrapper.intersect_with_coords(line_coords)
            bank_line_coords.append(coords_along_bank)
            bank_face_indices.append(face_indices)

        bank_order, data = self._link_lines_to_stations(bank_line_coords, bank_face_indices)

        return BankData.from_column_arrays(
            data,
            SingleBank,
            bank_lines=self.bank_lines,
            n_bank_lines=n_bank_lines,
            bank_order=bank_order,
        )

    def _link_lines_to_stations(self, bank_line_coords, bank_face_indices):
        # linking bank lines to chainage
        log_text("chainage_to_banks")
        river_center_line = self.river_data.river_center_line.as_array()
        n_bank_lines = len(self.bank_lines)

        bank_chainage_midpoints = [None] * n_bank_lines
        is_right_bank = [True] * n_bank_lines
        for bank_index, coords in enumerate(bank_line_coords):
            segment_mid_points = LineGeometry((coords[:-1, :] + coords[1:, :]) / 2)
            chainage_mid_points = segment_mid_points.intersect_with_line(
                river_center_line
            )

            # check if the bank line is defined from low chainage to high chainage
            if chainage_mid_points[0] > chainage_mid_points[-1]:
                # if not, flip the bank line and all associated data
                chainage_mid_points = chainage_mid_points[::-1]
                bank_line_coords[bank_index] = bank_line_coords[bank_index][::-1, :]
                bank_face_indices[bank_index] = bank_face_indices[bank_index][::-1]

            bank_chainage_midpoints[bank_index] = chainage_mid_points

            # check if the bank line is a left or right bank
            # when looking from low-to-high chainage
            is_right_bank[bank_index] = on_right_side(
                coords, river_center_line[:, :2]
            )
            if is_right_bank[bank_index]:
                log_text("right_side_bank", data={"ib": bank_index + 1})
            else:
                log_text("left_side_bank", data={"ib": bank_index + 1})

        bank_order = tuple("right" if val else "left" for val in is_right_bank)
        data = {
            'is_right_bank': is_right_bank,
            'bank_line_coords': bank_line_coords,
            'bank_face_indices': bank_face_indices,
            'bank_chainage_midpoints': bank_chainage_midpoints
        }
        return bank_order, data

__init__(river_data: ErosionRiverData, mesh_data: MeshData, verbose: bool = False)

Constructor for MeshProcessor.

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

def __init__(self, river_data: ErosionRiverData, mesh_data: MeshData, verbose: bool = False):
    """Constructor for MeshProcessor."""
    self.bank_lines = river_data.bank_lines
    self.mesh_data = mesh_data
    self.river_data = river_data
    self.wrapper = MeshWrapper(mesh_data, verbose=verbose)

get_bank_data() -> BankData

Intersect bank lines with a mesh and return bank data.

Returns:

Type	Description
BankData	BankData object containing bank line coordinates, face indices, and other bank-related data.

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

def get_bank_data(self) -> BankData:
    """Intersect bank lines with a mesh and return bank data.

    Returns:
        BankData object containing bank line coordinates, face indices, and other bank-related data.
    """
    n_bank_lines = len(self.bank_lines)

    bank_line_coords = []
    bank_face_indices = []
    for bank_index in range(n_bank_lines):
        line_coords = np.array(self.bank_lines.geometry[bank_index].coords)
        log_text("bank_nodes", data={"ib": bank_index + 1, "n": len(line_coords)})

        coords_along_bank, face_indices = self.wrapper.intersect_with_coords(line_coords)
        bank_line_coords.append(coords_along_bank)
        bank_face_indices.append(face_indices)

    bank_order, data = self._link_lines_to_stations(bank_line_coords, bank_face_indices)

    return BankData.from_column_arrays(
        data,
        SingleBank,
        bank_lines=self.bank_lines,
        n_bank_lines=n_bank_lines,
        bank_order=bank_order,
    )

MeshWrapper

A class for processing mesh-related operations.

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

class MeshWrapper:
    """A class for processing mesh-related operations."""

    def __init__(
        self, mesh_data: MeshData, d_thresh: float = 0.001, verbose: bool = False
    ):
        self.mesh_data = mesh_data
        self.d_thresh = d_thresh
        self.verbose = verbose

    def _read_target_object(self, xy_coords):
        self.given_coords = xy_coords
        self.intersection_state = IntersectionState(xy_coords.shape, self.verbose)

    def _handle_first_point(self, current_bank_point: np.ndarray):
        # get the point on the mesh face that is closest to the first bank point
        dx = self.mesh_data.x_face_coords - current_bank_point[0]
        dy = self.mesh_data.y_face_coords - current_bank_point[1]
        closest_cell_ind = np.nonzero(
            ~(
                (dx < 0).all(axis=1)
                | (dx > 0).all(axis=1)
                | (dy < 0).all(axis=1)
                | (dy > 0).all(axis=1)
            )
        )[0]
        index, vindex = self._find_starting_face(closest_cell_ind)
        self.intersection_state.current_face_index = index
        self.intersection_state.vertex_index = vindex
        self.intersection_state.update(current_bank_point)

    def _find_starting_face(self, face_indexes: np.ndarray):
        """Find the starting face for a bank line segment.

        This function determines the face index and vertex indices of the mesh
        that the first point of a bank line segment is associated with.

        Args:
            face_indexes (np.ndarray): Array of possible cell indices.
        """
        if len(face_indexes) == 0 and self.verbose:
            print("starting outside mesh")

        edges_indexes = self.mesh_data.locate_point(self.given_coords[0], face_indexes)
        if not isinstance(edges_indexes, list):
            # A single index was returned
            index = edges_indexes
            vindex = None
        else:
            # A list of indices is expected
            if edges_indexes:
                if self.verbose:
                    print(f"starting on edge of {edges_indexes}")

                index = -2 if len(edges_indexes) > 1 else edges_indexes[0]
                vindex = edges_indexes if len(edges_indexes) > 1 else None
            else:
                if self.verbose:
                    print("starting outside mesh")
                index = -1
                vindex = None

        return index, vindex

    def _log_slice_status(self, j, prev_b, bpj):
        if prev_b > 0:
            print(f"{j}: -- no further slices along this segment --")
        else:
            print(f"{j}: -- no slices along this segment --")

        if self.intersection_state.current_face_index >= 0:
            pnt = Point(bpj)
            polygon_k = self.mesh_data.get_face_by_index(self.intersection_state.current_face_index, as_polygon=True)

            if not polygon_k.contains(pnt):
                raise ValueError(
                    f"{j}: ERROR: point actually not contained within {self.intersection_state.current_face_index}!"
                )

    def _process_node_transition(self, segment: RiverSegment, node):
        """Process the transition at a node when a segment ends or continues."""
        finished = False

        if self.verbose:
            print(
                f"{segment.index}: moving via node {node} on edges {segment.edges} at {segment.distances[0]}"
            )

        # figure out where we will be heading afterwards ...
        if segment.distances[0] < 1.0:
            # segment passes through node and enter non-neighbouring cell ...
            # direction of current segment from bpj1 to bpj
            theta = segment.theta
        else:
            if (
                np.isclose(segment.distances[0], 1.0, rtol=RTOL, atol=ATOL)
                and segment.index == len(self.given_coords) - 1
            ):
                # catch case of last segment
                if self.verbose:
                    print(f"{segment.index}: last point ends in a node")

                self.intersection_state.update(segment.current_point)
                theta = 0.0
                finished = True
            else:
                # this segment ends in the node, so check next segment ...
                # direction of next segment from bpj to bp[j+1]
                theta = math.atan2(
                    self.given_coords[segment.index + 1][1]
                    - segment.current_point[1],
                    self.given_coords[segment.index + 1][0]
                    - segment.current_point[0],
                )
        next_face_index = None
        if not finished:
            next_face_index = self.mesh_data.resolve_next_face_by_direction(
                theta, node, segment.index
            )
        return False, next_face_index

    def _slice_by_node_or_edge(
        self, segment: RiverSegment, node, edge, faces
    ):
        finished = False
        next_face_index = None
        if node >= 0:
            # if we slice at a node ...
            finished, next_face_index = self._process_node_transition(segment, node)

        elif segment.distances[0] == 1:
            # ending at slice point, so ending on an edge ...
            if self.verbose:
                print(
                    f"{segment.index}: ending on edge {edge} at {segment.distances[0]}"
                )
            # figure out where we will be heading afterwards ...
            if segment.index == len(self.given_coords) - 1:
                # catch case of last segment
                if self.verbose:
                    print(f"{segment.index}: last point ends on an edge")
                self.intersection_state.update(segment.current_point)
                finished = True
            else:
                next_point = [self.given_coords[segment.index + 1][0], self.given_coords[segment.index + 1][1]]
                next_face_index = self.determine_next_face_on_edge(segment, next_point, edge, faces)

        return finished, next_face_index

    def _process_bank_segment(self, segment: RiverSegment):

        while True:
            if self.intersection_state.current_face_index == -2:
                index_src = self._resolve_ambiguous_edge_transition(segment, self.intersection_state.vertex_index)

                if len(index_src) == 1:
                    self.intersection_state.current_face_index = index_src[0]
                    self.intersection_state.vertex_index = index_src[0:1]
                else:
                    self.intersection_state.current_face_index = -2

            elif segment.is_length_zero():
                # segment has zero length
                break
            else:
                segment.distances, segment.edges, segment.nodes = (
                    self.mesh_data.find_segment_intersections(
                        self.intersection_state.current_face_index,
                        segment,
                    )
                )

            if len(segment.edges) == 0:
                # rest of segment associated with same face
                shape_length = self.intersection_state.point_index * SHAPE_MULTIPLIER
                self.intersection_state.update(
                    segment.current_point, shape_length=shape_length
                )
                break

            if len(segment.edges) > 1:
                segment.select_first_intersection()

            # slice location identified ...   (number of edges should be 1)
            node = segment.nodes[0]
            edge = segment.edges[0]
            faces = self.mesh_data.edge_face_connectivity[edge]
            segment.min_relative_distance = segment.distances[0]

            finished, face_index = self._slice_by_node_or_edge(
                segment,
                node,
                edge,
                faces,
            )
            if finished:
                break

            status = Status(**{
                "step": segment.index,
                "transition_index": node,
                "face_index": face_index,
                "prev_b": segment.min_relative_distance,
            })
            self.intersection_state.update_index_and_log(status, edge, faces)
            segment_x = (
                    segment.previous_point
                    + segment.min_relative_distance
                    * (segment.current_point - segment.previous_point)
            )
            shape_length = self.intersection_state.point_index * SHAPE_MULTIPLIER
            self.intersection_state.update(segment_x, shape_length=shape_length)
            if segment.min_relative_distance == 1:
                break

    def _resolve_ambiguous_edge_transition(self, segment: RiverSegment, vindex):
        """Resolve ambiguous edge transitions when a line segment is on the edge of multiple mesh faces."""
        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 = self.mesh_data.find_segment_intersections(i, segment)
            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)

        segment.edges, id_edges = np.unique(edges, return_index=True)
        segment.distances = b[id_edges]
        segment.nodes = nodes[id_edges]
        index_src = index_src[id_edges]

        return index_src

    def resolve_next_face_by_direction(
            self, theta: float, node, verbose_index: int = None
    ):
        """Helper to resolve the next face index based on the direction theta at a node."""

        if self.verbose:
            print(f"{verbose_index}: moving in direction theta = {theta}")

        edges = self.mesh_data.find_edges(theta, node, verbose_index)

        if self.verbose:
            print(f"{verbose_index}: the edge to the left is edge {edges.left}")
            print(f"{verbose_index}: the edge to the right is edge {edges.right}")

        if edges.left == edges.right:
            if self.verbose:
                print(f"{verbose_index}: continue along edge {edges.left}")

            next_face_index = self.mesh_data.edge_face_connectivity[edges.left, :]
        else:
            if self.verbose:
                print(
                    f"{verbose_index}: continue between edges {edges.left}"
                    f" on the left and {edges.right} on the right"
                )
            next_face_index = self.mesh_data.resolve_next_face_from_edges(
                node, edges, verbose_index
            )
        return next_face_index

    def determine_next_face_on_edge(
        self, segment: RiverSegment, next_point: List[float], edge, faces,
    ):
        """Determine the next face to continue along an edge based on the segment direction."""
        theta = math.atan2(
            next_point[1] - segment.current_point[1],
            next_point[0] - segment.current_point[0],
            )
        if self.verbose:
            print(f"{segment.index}: moving in direction theta = {theta}")

        theta_edge = self.mesh_data.calculate_edge_angle(edge)
        if theta == theta_edge or theta == -theta_edge:
            if self.verbose:
                print(f"{segment.index}: continue along edge {edge}")
            next_face_index = faces
        else:
            # check whether the (extended) segment slices any edge of faces[0]
            fe1 = self.mesh_data.face_edge_connectivity[faces[0]]
            reversed_segment = RiverSegment(
                index=segment.index,
                previous_point=segment.current_point,
                current_point=next_point,
                min_relative_distance=0,
            )
            _, _, edges = self.mesh_data.calculate_edge_intersections(
                fe1,
                reversed_segment,
                False,
            )
            # 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]
            # if there are no slices for faces[0], we continue towards faces[1]
            next_face_index = faces[0] if len(edges) > 0 else faces[1]
        return next_face_index

    def intersect_with_coords(self, given_coords: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """Intersects a coords 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.

        Returns:
            Tuple[np.ndarray, np.ndarray]:
                A tuple containing:
                - `coords` (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.
        """
        self._read_target_object(given_coords)
        for point_index, current_bank_point in enumerate(given_coords):
            if self.verbose:
                print(
                    f"Current location: {current_bank_point[0]}, {current_bank_point[1]}"
                )
            if point_index == 0:
                self._handle_first_point(current_bank_point)
            else:
                segment = RiverSegment(
                    index=point_index,
                    min_relative_distance=0,
                    current_point=current_bank_point,
                    previous_point=given_coords[point_index - 1],
                )
                self._process_bank_segment(segment)

        # clip to actual length (idx refers to segments, so we can ignore the last value)
        self.intersection_state.coords = self.intersection_state.coords[: self.intersection_state.point_index]
        self.intersection_state.face_indexes = self.intersection_state.face_indexes[: self.intersection_state.point_index - 1]

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

        # since index refers to segments, don't return the first one
        return self.intersection_state.coords, self.intersection_state.face_indexes

determine_next_face_on_edge(segment: RiverSegment, next_point: List[float], edge, faces)

Determine the next face to continue along an edge based on the segment direction.

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

def determine_next_face_on_edge(
    self, segment: RiverSegment, next_point: List[float], edge, faces,
):
    """Determine the next face to continue along an edge based on the segment direction."""
    theta = math.atan2(
        next_point[1] - segment.current_point[1],
        next_point[0] - segment.current_point[0],
        )
    if self.verbose:
        print(f"{segment.index}: moving in direction theta = {theta}")

    theta_edge = self.mesh_data.calculate_edge_angle(edge)
    if theta == theta_edge or theta == -theta_edge:
        if self.verbose:
            print(f"{segment.index}: continue along edge {edge}")
        next_face_index = faces
    else:
        # check whether the (extended) segment slices any edge of faces[0]
        fe1 = self.mesh_data.face_edge_connectivity[faces[0]]
        reversed_segment = RiverSegment(
            index=segment.index,
            previous_point=segment.current_point,
            current_point=next_point,
            min_relative_distance=0,
        )
        _, _, edges = self.mesh_data.calculate_edge_intersections(
            fe1,
            reversed_segment,
            False,
        )
        # 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]
        # if there are no slices for faces[0], we continue towards faces[1]
        next_face_index = faces[0] if len(edges) > 0 else faces[1]
    return next_face_index

intersect_with_coords(given_coords: np.ndarray) -> Tuple[np.ndarray, np.ndarray]

Intersects a coords 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.

Returns:

Type	Description
Tuple[ndarray, ndarray]	Tuple[np.ndarray, np.ndarray]: A tuple containing: - coords (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_with_coords(self, given_coords: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """Intersects a coords 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.

    Returns:
        Tuple[np.ndarray, np.ndarray]:
            A tuple containing:
            - `coords` (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.
    """
    self._read_target_object(given_coords)
    for point_index, current_bank_point in enumerate(given_coords):
        if self.verbose:
            print(
                f"Current location: {current_bank_point[0]}, {current_bank_point[1]}"
            )
        if point_index == 0:
            self._handle_first_point(current_bank_point)
        else:
            segment = RiverSegment(
                index=point_index,
                min_relative_distance=0,
                current_point=current_bank_point,
                previous_point=given_coords[point_index - 1],
            )
            self._process_bank_segment(segment)

    # clip to actual length (idx refers to segments, so we can ignore the last value)
    self.intersection_state.coords = self.intersection_state.coords[: self.intersection_state.point_index]
    self.intersection_state.face_indexes = self.intersection_state.face_indexes[: self.intersection_state.point_index - 1]

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

    # since index refers to segments, don't return the first one
    return self.intersection_state.coords, self.intersection_state.face_indexes

resolve_next_face_by_direction(theta: float, node, verbose_index: int = None)

Helper to resolve the next face index based on the direction theta at a node.

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

def resolve_next_face_by_direction(
        self, theta: float, node, verbose_index: int = None
):
    """Helper to resolve the next face index based on the direction theta at a node."""

    if self.verbose:
        print(f"{verbose_index}: moving in direction theta = {theta}")

    edges = self.mesh_data.find_edges(theta, node, verbose_index)

    if self.verbose:
        print(f"{verbose_index}: the edge to the left is edge {edges.left}")
        print(f"{verbose_index}: the edge to the right is edge {edges.right}")

    if edges.left == edges.right:
        if self.verbose:
            print(f"{verbose_index}: continue along edge {edges.left}")

        next_face_index = self.mesh_data.edge_face_connectivity[edges.left, :]
    else:
        if self.verbose:
            print(
                f"{verbose_index}: continue between edges {edges.left}"
                f" on the left and {edges.right} on the right"
            )
        next_face_index = self.mesh_data.resolve_next_face_from_edges(
            node, edges, verbose_index
        )
    return next_face_index

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.

Data Models

The Bank Erosion module uses several data models to represent inputs, calculation parameters, and results:

Calculation Data Models

dfastbe.bank_erosion.data_models.calculation

Erosion-related data structures.

BankData dataclass

Bases: BaseBank[SingleBank]

Class to hold bank-related data.

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]

bank_chainage_midpoints: List[np.ndarray] property

Get the chainage midpoints of the bank lines.

bank_line_coords: List[np.ndarray] property

Get the coordinates of the bank lines.

height: List[np.ndarray] property

Get the bank height.

is_right_bank: List[bool] property

Get the bank direction.

num_stations_per_bank: List[int] property

Get the number of stations per bank.

BaseBank dataclass

Bases: Generic[GenericType]

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

@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]

bank_protection_level: List[np.ndarray] property

Get the bank protection level.

tauc: List[np.ndarray] property

Get the critical bank shear stress values.

ErosionResults dataclass

Class to hold erosion results.

Parameters:

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)

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.

Input Data Models

dfastbe.bank_erosion.data_models.inputs

BankLinesResultsError

Bases: Exception

Custom exception for BankLine results errors.

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 = 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

    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_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_center_line_arr = self.river_center_line.as_array()
        river_axis_km = river_axis.intersect_with_line(river_center_line_arr)

        # clip river axis to reach of interest (get the closest point to the first and last station)
        i1 = np.argmin(
            ((river_center_line_arr[0, :2] - river_axis_numpy) ** 2).sum(axis=1)
        )
        i2 = np.argmin(
            ((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 = LineGeometry(river_axis_numpy, crs=self.config_file.crs)
        river_axis.add_data(data={"stations": river_axis_km})
        return river_axis

process_river_axis_by_center_line() -> 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.

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

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_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_center_line_arr = self.river_center_line.as_array()
    river_axis_km = river_axis.intersect_with_line(river_center_line_arr)

    # clip river axis to reach of interest (get the closest point to the first and last station)
    i1 = np.argmin(
        ((river_center_line_arr[0, :2] - river_axis_numpy) ** 2).sum(axis=1)
    )
    i2 = np.argmin(
        ((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 = LineGeometry(river_axis_numpy, crs=self.config_file.crs)
    river_axis.add_data(data={"stations": river_axis_km})
    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, verbose: bool = False) -> 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,
            verbose=verbose,
        )

    @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(verbose: bool = False) -> 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, verbose: bool = False) -> 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,
        verbose=verbose,
    )

ShipsParameters dataclass

Data for ships going through the fairway for bank erosion simulation.

Parameters:

Name	Type	Description	Default
config_file	ConfigFile	Configuration file containing parameters.	required
velocity	List[ndarray]	Ship velocities for each bank.	required
number	List[ndarray]	Number of ships for each bank.	required
num_waves	List[ndarray]	Number of waves per ship for each bank.	required
draught	List[ndarray]	Draught of ships for each bank.	required
type	List[ndarray]	Type of ships for each bank.	required
slope	List[ndarray]	Slope values for each bank.	required
reed	List[ndarray]	Reed values for each bank.	required

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

@dataclass
class ShipsParameters:
    """Data for ships going through the fairway for bank erosion simulation.

    Args:
        config_file (ConfigFile):
            Configuration file containing parameters.
        velocity (List[np.ndarray]):
            Ship velocities for each bank.
        number (List[np.ndarray]):
            Number of ships for each bank.
        num_waves (List[np.ndarray]):
            Number of waves per ship for each bank.
        draught (List[np.ndarray]):
            Draught of ships for each bank.
        type (List[np.ndarray]):
            Type of ships for each bank.
        slope (List[np.ndarray]):
            Slope values for each bank.
        reed (List[np.ndarray]):
            Reed values for each bank.
    """

    config_file: ConfigFile
    velocity: List[np.ndarray]
    number: List[np.ndarray]
    num_waves: List[np.ndarray]
    draught: List[np.ndarray]
    type: List[np.ndarray]
    slope: List[np.ndarray]
    reed: List[np.ndarray]

    REED_DAMPING_COEFFICIENT = 8.5e-4
    REED_DAMPING_EXPONENT = 0.8

    @classmethod
    def get_ship_data(
        cls, num_stations_per_bank: List[int], config_file: ConfigFile
    ) -> "ShipsParameters":
        """Get ship parameters from the configuration file.

        Args:
            num_stations_per_bank (List[int]):
                The number of stations per bank.
            config_file (ConfigFile):
                Configuration file containing parameters.

        Returns:
            ShipsParameters: An instance of ShipsParameters with parameters read from the config file.
        """

        param_defs = cls._get_initial_parameter_definitions()
        param_resolved = cls._get_parameters(
            config_file, num_stations_per_bank, param_defs
        )

        param_dict = {
            "velocity": param_resolved["vship"],
            "number": param_resolved["nship"],
            "num_waves": param_resolved["nwave"],
            "draught": param_resolved["draught"],
            "type": param_resolved["shiptype"],
            "slope": param_resolved["slope"],
            "reed": param_resolved["reed"],
        }

        return cls(config_file, **param_dict)

    @staticmethod
    def _get_parameters(
        config_file: ConfigFile,
        num_stations_per_bank: List[int],
        param_defs: List[Parameters],
        level_i_str: str = "",
    ) -> Dict[str, np.ndarray]:
        """Resolve a list of Parameters with values from the configuration file."""
        param_resolved = {}

        for param in param_defs:
            param_resolved[f"{param.name.lower()}"] = config_file.get_parameter(
                "Erosion",
                f"{param.name}{level_i_str}",
                num_stations_per_bank,
                default=param.default,
                valid=param.valid,
                onefile=param.onefile,
                positive=param.positive,
                ext=param.ext,
            )
        return param_resolved

    @classmethod
    def _get_initial_parameter_definitions(cls) -> List[Parameters]:
        """Get parameter definitions for discharge parameters.

        Returns:
            List[namedtuple]: List of parameter definitions.
        """
        return [
            Parameters("VShip", None, None, True, True, None),
            Parameters("NShip", None, None, True, True, None),
            Parameters("NWave", 5, None, True, True, None),
            Parameters("Draught", None, None, True, True, None),
            Parameters("ShipType", None, [1, 2, 3], True, None, None),
            Parameters("Slope", 20, None, None, True, "slp"),
            Parameters("Reed", 0, None, None, True, "rdd"),
        ]

    def _get_discharge_parameter_definitions(self) -> List[Parameters]:
        """Get parameter definitions for discharge parameters.

        Returns:
            List[namedtuple]: List of parameter definitions.
        """
        return [
            Parameters("VShip", self.velocity, None, None, None, None),
            Parameters("NShip", self.number, None, None, None, None),
            Parameters("NWave", self.num_waves, None, None, None, None),
            Parameters("Draught", self.draught, None, None, None, None),
            Parameters("ShipType", self.type, [1, 2, 3], True, None, None),
            Parameters("Slope", self.slope, None, None, True, "slp"),
            Parameters("Reed", self.reed, None, None, True, "rdd"),
        ]

    @staticmethod
    def _calculate_ship_derived_parameters(
        slope_values: List[np.ndarray], reed_values: List[np.ndarray]
    ) -> Tuple[List[np.ndarray], List[np.ndarray]]:
        """Calculate derived parameters from slope and reed values.

        Args:
            slope_values (List[np.ndarray]): Slope values for each bank.
            reed_values (List[np.ndarray]): Reed values for each bank.

        Returns:
            Tuple[List[np.ndarray], List[np.ndarray]]: Calculated mu_slope and mu_reed values.
        """
        mu_slope = []
        mu_reed = []

        for ps, pr in zip(slope_values, reed_values):
            # Calculate mu_slope (inverse of slope for non-zero values)
            mus = ps.copy()
            mus[mus > 0] = 1.0 / mus[mus > 0]
            mu_slope.append(mus)

            # Calculate mu_reed (empirical damping coefficient)
            mu_reed.append(
                ShipsParameters.REED_DAMPING_COEFFICIENT
                * pr**ShipsParameters.REED_DAMPING_EXPONENT
            )

        return mu_slope, mu_reed

    def read_discharge_parameters(
        self,
        level_i: int,
        num_stations_per_bank: List[int],
    ) -> SingleLevelParameters:
        """Read Discharge level parameters.

        Read all discharge-specific input arrays for level_i.

        Args:
            level_i (int):
                The index of the discharge level.
            num_stations_per_bank (List[int]):
                The number of stations per bank.

        Returns:
            SingleLevelParameters: The discharge level parameters.
        """
        param_defs = self._get_discharge_parameter_definitions()
        param_resolved = self._get_parameters(
            self.config_file, num_stations_per_bank, param_defs, f"{level_i + 1}"
        )

        mu_slope, mu_reed = self._calculate_ship_derived_parameters(
            param_resolved["slope"], param_resolved["reed"]
        )

        return SingleLevelParameters.from_column_arrays(
            {
                "id": level_i,
                "ship_velocity": param_resolved["vship"],
                "num_ship": param_resolved["nship"],
                "num_waves_per_ship": param_resolved["nwave"],
                "ship_draught": param_resolved["draught"],
                "ship_type": param_resolved["shiptype"],
                "par_slope": param_resolved["slope"],
                "par_reed": param_resolved["reed"],
                "mu_slope": mu_slope,
                "mu_reed": mu_reed,
            },
            SingleParameters,
        )

get_ship_data(num_stations_per_bank: List[int], config_file: ConfigFile) -> ShipsParameters classmethod

Get ship parameters from the configuration file.

Parameters:

Name	Type	Description	Default
num_stations_per_bank	List[int]	The number of stations per bank.	required
config_file	ConfigFile	Configuration file containing parameters.	required

Returns:

Name	Type	Description
ShipsParameters	ShipsParameters	An instance of ShipsParameters with parameters read from the config file.

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

@classmethod
def get_ship_data(
    cls, num_stations_per_bank: List[int], config_file: ConfigFile
) -> "ShipsParameters":
    """Get ship parameters from the configuration file.

    Args:
        num_stations_per_bank (List[int]):
            The number of stations per bank.
        config_file (ConfigFile):
            Configuration file containing parameters.

    Returns:
        ShipsParameters: An instance of ShipsParameters with parameters read from the config file.
    """

    param_defs = cls._get_initial_parameter_definitions()
    param_resolved = cls._get_parameters(
        config_file, num_stations_per_bank, param_defs
    )

    param_dict = {
        "velocity": param_resolved["vship"],
        "number": param_resolved["nship"],
        "num_waves": param_resolved["nwave"],
        "draught": param_resolved["draught"],
        "type": param_resolved["shiptype"],
        "slope": param_resolved["slope"],
        "reed": param_resolved["reed"],
    }

    return cls(config_file, **param_dict)

read_discharge_parameters(level_i: int, num_stations_per_bank: List[int]) -> SingleLevelParameters

Read Discharge level parameters.

Read all discharge-specific input arrays for level_i.

Parameters:

Name	Type	Description	Default
level_i	int	The index of the discharge level.	required
num_stations_per_bank	List[int]	The number of stations per bank.	required

Returns:

Name	Type	Description
SingleLevelParameters	SingleLevelParameters	The discharge level parameters.

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

def read_discharge_parameters(
    self,
    level_i: int,
    num_stations_per_bank: List[int],
) -> SingleLevelParameters:
    """Read Discharge level parameters.

    Read all discharge-specific input arrays for level_i.

    Args:
        level_i (int):
            The index of the discharge level.
        num_stations_per_bank (List[int]):
            The number of stations per bank.

    Returns:
        SingleLevelParameters: The discharge level parameters.
    """
    param_defs = self._get_discharge_parameter_definitions()
    param_resolved = self._get_parameters(
        self.config_file, num_stations_per_bank, param_defs, f"{level_i + 1}"
    )

    mu_slope, mu_reed = self._calculate_ship_derived_parameters(
        param_resolved["slope"], param_resolved["reed"]
    )

    return SingleLevelParameters.from_column_arrays(
        {
            "id": level_i,
            "ship_velocity": param_resolved["vship"],
            "num_ship": param_resolved["nship"],
            "num_waves_per_ship": param_resolved["nwave"],
            "ship_draught": param_resolved["draught"],
            "ship_type": param_resolved["shiptype"],
            "par_slope": param_resolved["slope"],
            "par_reed": param_resolved["reed"],
            "mu_slope": mu_slope,
            "mu_reed": mu_reed,
        },
        SingleParameters,
    )

For more details, see Bank Erosion Input Data Models.

Workflow

The typical workflow for bank erosion calculation is:

1. Initialize the Erosion object with a configuration file
2. Call the run method to start the erosion calculation process
3. The run method orchestrates the entire process:
4. Processes the river axis
5. Gets fairway data
6. Calculates bank-fairway distance
7. Prepares initial conditions
8. Processes discharge levels
9. Computes erosion per level
10. Post-processes results
11. Writes output files
12. 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.