smashbox.model.mesh 

  1import smash
  2import pandas as pd
  3import geopandas as gpd
  4import os
  5import numpy as np
  6from smashbox.init.param import smashboxparam
  7import copy
  8
  9
 10class mesh:
 11    """Class mesh(). This class has three functions and one attributes to generate, load,
 12    write and store the mesh used by Smash."""
 13
 14    def __init__(self, setup):
 15
 16        self._setup = setup
 17
 18        self.mesh = None
 19        self.catchment_polygon = None
 20
 21    def write_mesh(self, filename: os.PathLike | None = None):
 22        """
 23        Write the mesh of the Smash model unsing hdf5 format.
 24        :param filename: Path to the file where to write the mesh, defaults to None
 25        :type filename: os.PathLike | None, optional
 26
 27        """
 28        if filename is not None:
 29            smash.io.save_mesh(self.mesh, path=filename)
 30        else:
 31            raise ValueError(f"Output filename '{filename}' is None.")
 32
 33    def load_mesh(self, filename: os.PathLike | None = None):
 34        """
 35        Read a mesh for Smash stored with the hdf5 format.
 36        :param filename: path to the hdf5 file
 37        :type filename: TYPE
 38
 39        """
 40
 41        if os.path.exists(filename):
 42            self.mesh = smash.io.read_mesh(filename)
 43        else:
 44            raise ValueError(f"{filename} does not exist.")
 45
 46    def generate_mesh(
 47        self,
 48        param: smashboxparam | None = None,
 49        query: str | None = None,
 50        max_depth: float = 1.0,
 51        area_error_th: None | float = None,
 52        lacuna_threshold: None | float = None,
 53    ):
 54        """
 55        :param param: Class param.smashboxparam(), store main smashbox parameters
 56        :type param: param.smashboxparam()
 57        :param query: Any pandas dataframe query as string: '(SURF>20) & (SURF<100)'. This query
 58        must be build using the field (column name) in the outlet database.
 59        https://pandas.pydata.org/docs/user_guide/indexing.html#the-query-method
 60        :type query: str
 61        :max_depth: The maximum depth accepted by the algorithm to find the catchment outlet.
 62            A **max_depth** of 1 means that the algorithm will search among the
 63            combinations in
 64            (``row - 1``, ``row``, ``row + 1``; ``col - 1``, ``col``, ``col + 1``),
 65            the coordinates that minimize
 66            the relative error between the given catchment area and the modeled
 67            catchment area calculated from the
 68            flow directions file.
 69        :type `int`, default 1
 70        :param area_error_th: Tolerance error during the positionning of the outlets. If the Error `(Ssim-Sobs)/Sobs > area_error_th`, the outlet will be excluded.
 71        :type area_error_th: float
 72        :param lacuna_threshold: Lacuna threshold for the discharges in percent. All gauge where the proportion of lacuna exceed this threshold will be removed.
 73        :type: float | Nonetype
 74
 75        """
 76
 77        if param.bbox is not None:
 78            bbox = [
 79                param.bbox["left"],
 80                param.bbox["right"],
 81                param.bbox["bottom"],
 82                param.bbox["top"],
 83            ]
 84        else:
 85            bbox = None
 86
 87        if not os.path.exists(param.outlets_database):
 88            param.outlets_database = os.path.join(
 89                param.asset_dir, "outlets", param.outlets_database
 90            )
 91
 92        if os.path.exists(param.outlets_database):
 93            stations_calage = pd.read_csv(param.outlets_database)
 94        else:
 95            raise ValueError(f"</> Error: file {param.outlets_database} not found")
 96
 97        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
 98        Input_outletsID = param.outletsID
 99
100        if param.outletsID is not None:
101            if len(param.outletsID) > 0:
102                stations_calage = stations_calage.loc[
103                    stations_calage[param.outlets_database_fields["id"]].isin(
104                        param.outletsID
105                    )
106                ]
107
108            # check and print outlets not in database
109            for index, id_outlets in enumerate(Input_outletsID):
110
111                if (
112                    id_outlets
113                    not in stations_calage[param.outlets_database_fields["id"]].values
114                ):
115
116                    print(
117                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
118                    )
119
120            if query is not None:
121                stations_calage = stations_calage.query(query)
122
123            if bbox is not None:
124                stations_calage = stations_calage[
125                    (
126                        stations_calage[param.outlets_database_fields["coord_x"]]
127                        >= param.bbox["left"]
128                    )
129                    & (
130                        stations_calage[param.outlets_database_fields["coord_x"]]
131                        <= param.bbox["right"]
132                    )
133                    & (
134                        stations_calage[param.outlets_database_fields["coord_y"]]
135                        >= param.bbox["bottom"]
136                    )
137                    & (
138                        stations_calage[param.outlets_database_fields["coord_y"]]
139                        <= param.bbox["top"]
140                    )
141                ]
142
143            # check and print outlets exluded outside bbox
144            for index, id_outlets in enumerate(Input_outletsID):
145                if (
146                    id_outlets
147                    not in stations_calage[param.outlets_database_fields["id"]].values
148                ):
149                    print(
150                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
151                    )
152
153            if len(stations_calage) == 0:
154                print(
155                    f"</> Error: outlets {param.outletsID} not found in"
156                    "{param.outlets_database}"
157                )
158                raise ValueError(
159                    f"</> Error: outlets {param.outletsID} not found in"
160                    "{param.outlets_database}"
161                )
162
163            columns = {
164                "coord_x": param.outlets_database_fields["coord_x"],
165                "coord_y": param.outlets_database_fields["coord_y"],
166                "area": param.outlets_database_fields["area"],
167                "id": param.outlets_database_fields["id"],
168            }
169
170            # fist build
171            self.mesh = smash.factory.generate_mesh(
172                flwdir_path=param.flowdir,
173                bbox=bbox,
174                x=np.array(stations_calage[columns["coord_x"]][:]),
175                y=np.array(stations_calage[columns["coord_y"]][:]),
176                area=np.array(
177                    stations_calage[columns["area"]][:] * 1e6
178                ),  # Convert km² to m²
179                code=np.array(stations_calage[columns["id"]][:]),
180                epsg=param.epsg,
181                shp_path=param.outlets_shapefile,
182                max_depth=max_depth,
183                area_error_th=area_error_th,
184            )
185
186            if lacuna_threshold is not None:
187                setup = copy.deepcopy(self._setup.setup)
188                setup.update(
189                    {
190                        "read_prcp": False,
191                        "read_pet": False,
192                        "read_qobs": True,
193                        "adjust_interception": False,
194                        "compute_mean_atmos": False,
195                    }
196                )
197
198                # filter gauge if lacuna exceed a thresholds
199                model = smash.Model(setup, self.mesh)
200                qobs = model.response_data.q
201
202                valid_gauge = model.mesh.code[
203                    np.where(
204                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
205                        >= (100.0 - lacuna_threshold)
206                    )
207                ]
208
209                unvalid_gauge = model.mesh.code[
210                    np.where(
211                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
212                        > lacuna_threshold
213                    )
214                ]
215
216                if len(unvalid_gauge) > 0:
217                    print(
218                        f"</> Remove gauges from the mesh where total lacuna "
219                        f"between {setup['start_time']} and {setup['end_time']} "
220                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
221                    )
222
223                stations_calage = stations_calage.loc[
224                    stations_calage[columns["id"]].isin(valid_gauge)
225                ]
226
227                if len(stations_calage) == 0:
228                    print("</> Warnings, no outlets/gauge will be added to the mesh !")
229
230                # rebuild mesh
231                self.mesh = smash.factory.generate_mesh(
232                    flwdir_path=param.flowdir,
233                    bbox=bbox,
234                    x=np.array(stations_calage[columns["coord_x"]][:]),
235                    y=np.array(stations_calage[columns["coord_y"]][:]),
236                    area=np.array(
237                        stations_calage[columns["area"]][:] * 1e6
238                    ),  # Convert km² to m²
239                    code=np.array(stations_calage[columns["id"]][:]),
240                    epsg=param.epsg,
241                    shp_path=param.outlets_shapefile,
242                    max_depth=max_depth,
243                    area_error_th=area_error_th,
244                )
245
246            if param.outlets_shapefile is not None:
247                print("</> Outlets shapefile detected. Loading outlets ...")
248                if param.outlets_database_fields["id_shapefile"] != "None":
249                    col_id = param.outlets_database_fields["id_shapefile"]
250                    code_shape_file = []
251                    for i in range(len(self.mesh["code"])):
252                        sta = self.mesh["code"][i]
253                        code_shape_file.extend(
254                            stations_calage.loc[
255                                stations_calage[columns["id"]] == sta, col_id
256                            ].to_list()
257                        )
258                    code_shape_file = np.array(code_shape_file)
259                else:
260                    code_shape_file = self.mesh["code"]
261
262                catchment_polygon = gpd.read_file(param.outlets_shapefile)
263                self.catchment_polygon = catchment_polygon.loc[
264                    catchment_polygon.code.isin(code_shape_file)
265                ]
266                del catchment_polygon
267
268        else:
269            stations_calage = pd.DataFrame(None)
270
271            if bbox is None:
272                raise ValueError(
273                    "Bbox is None. If no outlets provided, the bbox must be defined."
274                )
275
276            self.mesh = smash.factory.generate_mesh(
277                flwdir_path=param.flowdir, bbox=bbox, epsg=param.epsg, max_depth=max_depth
278            )
class mesh: 
 11class mesh:
 12    """Class mesh(). This class has three functions and one attributes to generate, load,
 13    write and store the mesh used by Smash."""
 14
 15    def __init__(self, setup):
 16
 17        self._setup = setup
 18
 19        self.mesh = None
 20        self.catchment_polygon = None
 21
 22    def write_mesh(self, filename: os.PathLike | None = None):
 23        """
 24        Write the mesh of the Smash model unsing hdf5 format.
 25        :param filename: Path to the file where to write the mesh, defaults to None
 26        :type filename: os.PathLike | None, optional
 27
 28        """
 29        if filename is not None:
 30            smash.io.save_mesh(self.mesh, path=filename)
 31        else:
 32            raise ValueError(f"Output filename '{filename}' is None.")
 33
 34    def load_mesh(self, filename: os.PathLike | None = None):
 35        """
 36        Read a mesh for Smash stored with the hdf5 format.
 37        :param filename: path to the hdf5 file
 38        :type filename: TYPE
 39
 40        """
 41
 42        if os.path.exists(filename):
 43            self.mesh = smash.io.read_mesh(filename)
 44        else:
 45            raise ValueError(f"{filename} does not exist.")
 46
 47    def generate_mesh(
 48        self,
 49        param: smashboxparam | None = None,
 50        query: str | None = None,
 51        max_depth: float = 1.0,
 52        area_error_th: None | float = None,
 53        lacuna_threshold: None | float = None,
 54    ):
 55        """
 56        :param param: Class param.smashboxparam(), store main smashbox parameters
 57        :type param: param.smashboxparam()
 58        :param query: Any pandas dataframe query as string: '(SURF>20) & (SURF<100)'. This query
 59        must be build using the field (column name) in the outlet database.
 60        https://pandas.pydata.org/docs/user_guide/indexing.html#the-query-method
 61        :type query: str
 62        :max_depth: The maximum depth accepted by the algorithm to find the catchment outlet.
 63            A **max_depth** of 1 means that the algorithm will search among the
 64            combinations in
 65            (``row - 1``, ``row``, ``row + 1``; ``col - 1``, ``col``, ``col + 1``),
 66            the coordinates that minimize
 67            the relative error between the given catchment area and the modeled
 68            catchment area calculated from the
 69            flow directions file.
 70        :type `int`, default 1
 71        :param area_error_th: Tolerance error during the positionning of the outlets. If the Error `(Ssim-Sobs)/Sobs > area_error_th`, the outlet will be excluded.
 72        :type area_error_th: float
 73        :param lacuna_threshold: Lacuna threshold for the discharges in percent. All gauge where the proportion of lacuna exceed this threshold will be removed.
 74        :type: float | Nonetype
 75
 76        """
 77
 78        if param.bbox is not None:
 79            bbox = [
 80                param.bbox["left"],
 81                param.bbox["right"],
 82                param.bbox["bottom"],
 83                param.bbox["top"],
 84            ]
 85        else:
 86            bbox = None
 87
 88        if not os.path.exists(param.outlets_database):
 89            param.outlets_database = os.path.join(
 90                param.asset_dir, "outlets", param.outlets_database
 91            )
 92
 93        if os.path.exists(param.outlets_database):
 94            stations_calage = pd.read_csv(param.outlets_database)
 95        else:
 96            raise ValueError(f"</> Error: file {param.outlets_database} not found")
 97
 98        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
 99        Input_outletsID = param.outletsID
100
101        if param.outletsID is not None:
102            if len(param.outletsID) > 0:
103                stations_calage = stations_calage.loc[
104                    stations_calage[param.outlets_database_fields["id"]].isin(
105                        param.outletsID
106                    )
107                ]
108
109            # check and print outlets not in database
110            for index, id_outlets in enumerate(Input_outletsID):
111
112                if (
113                    id_outlets
114                    not in stations_calage[param.outlets_database_fields["id"]].values
115                ):
116
117                    print(
118                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
119                    )
120
121            if query is not None:
122                stations_calage = stations_calage.query(query)
123
124            if bbox is not None:
125                stations_calage = stations_calage[
126                    (
127                        stations_calage[param.outlets_database_fields["coord_x"]]
128                        >= param.bbox["left"]
129                    )
130                    & (
131                        stations_calage[param.outlets_database_fields["coord_x"]]
132                        <= param.bbox["right"]
133                    )
134                    & (
135                        stations_calage[param.outlets_database_fields["coord_y"]]
136                        >= param.bbox["bottom"]
137                    )
138                    & (
139                        stations_calage[param.outlets_database_fields["coord_y"]]
140                        <= param.bbox["top"]
141                    )
142                ]
143
144            # check and print outlets exluded outside bbox
145            for index, id_outlets in enumerate(Input_outletsID):
146                if (
147                    id_outlets
148                    not in stations_calage[param.outlets_database_fields["id"]].values
149                ):
150                    print(
151                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
152                    )
153
154            if len(stations_calage) == 0:
155                print(
156                    f"</> Error: outlets {param.outletsID} not found in"
157                    "{param.outlets_database}"
158                )
159                raise ValueError(
160                    f"</> Error: outlets {param.outletsID} not found in"
161                    "{param.outlets_database}"
162                )
163
164            columns = {
165                "coord_x": param.outlets_database_fields["coord_x"],
166                "coord_y": param.outlets_database_fields["coord_y"],
167                "area": param.outlets_database_fields["area"],
168                "id": param.outlets_database_fields["id"],
169            }
170
171            # fist build
172            self.mesh = smash.factory.generate_mesh(
173                flwdir_path=param.flowdir,
174                bbox=bbox,
175                x=np.array(stations_calage[columns["coord_x"]][:]),
176                y=np.array(stations_calage[columns["coord_y"]][:]),
177                area=np.array(
178                    stations_calage[columns["area"]][:] * 1e6
179                ),  # Convert km² to m²
180                code=np.array(stations_calage[columns["id"]][:]),
181                epsg=param.epsg,
182                shp_path=param.outlets_shapefile,
183                max_depth=max_depth,
184                area_error_th=area_error_th,
185            )
186
187            if lacuna_threshold is not None:
188                setup = copy.deepcopy(self._setup.setup)
189                setup.update(
190                    {
191                        "read_prcp": False,
192                        "read_pet": False,
193                        "read_qobs": True,
194                        "adjust_interception": False,
195                        "compute_mean_atmos": False,
196                    }
197                )
198
199                # filter gauge if lacuna exceed a thresholds
200                model = smash.Model(setup, self.mesh)
201                qobs = model.response_data.q
202
203                valid_gauge = model.mesh.code[
204                    np.where(
205                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
206                        >= (100.0 - lacuna_threshold)
207                    )
208                ]
209
210                unvalid_gauge = model.mesh.code[
211                    np.where(
212                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
213                        > lacuna_threshold
214                    )
215                ]
216
217                if len(unvalid_gauge) > 0:
218                    print(
219                        f"</> Remove gauges from the mesh where total lacuna "
220                        f"between {setup['start_time']} and {setup['end_time']} "
221                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
222                    )
223
224                stations_calage = stations_calage.loc[
225                    stations_calage[columns["id"]].isin(valid_gauge)
226                ]
227
228                if len(stations_calage) == 0:
229                    print("</> Warnings, no outlets/gauge will be added to the mesh !")
230
231                # rebuild mesh
232                self.mesh = smash.factory.generate_mesh(
233                    flwdir_path=param.flowdir,
234                    bbox=bbox,
235                    x=np.array(stations_calage[columns["coord_x"]][:]),
236                    y=np.array(stations_calage[columns["coord_y"]][:]),
237                    area=np.array(
238                        stations_calage[columns["area"]][:] * 1e6
239                    ),  # Convert km² to m²
240                    code=np.array(stations_calage[columns["id"]][:]),
241                    epsg=param.epsg,
242                    shp_path=param.outlets_shapefile,
243                    max_depth=max_depth,
244                    area_error_th=area_error_th,
245                )
246
247            if param.outlets_shapefile is not None:
248                print("</> Outlets shapefile detected. Loading outlets ...")
249                if param.outlets_database_fields["id_shapefile"] != "None":
250                    col_id = param.outlets_database_fields["id_shapefile"]
251                    code_shape_file = []
252                    for i in range(len(self.mesh["code"])):
253                        sta = self.mesh["code"][i]
254                        code_shape_file.extend(
255                            stations_calage.loc[
256                                stations_calage[columns["id"]] == sta, col_id
257                            ].to_list()
258                        )
259                    code_shape_file = np.array(code_shape_file)
260                else:
261                    code_shape_file = self.mesh["code"]
262
263                catchment_polygon = gpd.read_file(param.outlets_shapefile)
264                self.catchment_polygon = catchment_polygon.loc[
265                    catchment_polygon.code.isin(code_shape_file)
266                ]
267                del catchment_polygon
268
269        else:
270            stations_calage = pd.DataFrame(None)
271
272            if bbox is None:
273                raise ValueError(
274                    "Bbox is None. If no outlets provided, the bbox must be defined."
275                )
276
277            self.mesh = smash.factory.generate_mesh(
278                flwdir_path=param.flowdir, bbox=bbox, epsg=param.epsg, max_depth=max_depth
279            )

Class mesh(). This class has three functions and one attributes to generate, load, write and store the mesh used by Smash.

mesh(setup) 
15    def __init__(self, setup):
16
17        self._setup = setup
18
19        self.mesh = None
20        self.catchment_polygon = None
mesh
catchment_polygon
def write_mesh(self, filename: os.PathLike | None = None): 
22    def write_mesh(self, filename: os.PathLike | None = None):
23        """
24        Write the mesh of the Smash model unsing hdf5 format.
25        :param filename: Path to the file where to write the mesh, defaults to None
26        :type filename: os.PathLike | None, optional
27
28        """
29        if filename is not None:
30            smash.io.save_mesh(self.mesh, path=filename)
31        else:
32            raise ValueError(f"Output filename '{filename}' is None.")

Write the mesh of the Smash model unsing hdf5 format.

Parameters
  • filename: Path to the file where to write the mesh, defaults to None
def load_mesh(self, filename: os.PathLike | None = None): 
34    def load_mesh(self, filename: os.PathLike | None = None):
35        """
36        Read a mesh for Smash stored with the hdf5 format.
37        :param filename: path to the hdf5 file
38        :type filename: TYPE
39
40        """
41
42        if os.path.exists(filename):
43            self.mesh = smash.io.read_mesh(filename)
44        else:
45            raise ValueError(f"{filename} does not exist.")

Read a mesh for Smash stored with the hdf5 format.

Parameters
  • filename: path to the hdf5 file
def generate_mesh( self, param: smashbox.init.param.smashboxparam | None = None, query: str | None = None, max_depth: float = 1.0, area_error_th: None | float = None, lacuna_threshold: None | float = None): 
 47    def generate_mesh(
 48        self,
 49        param: smashboxparam | None = None,
 50        query: str | None = None,
 51        max_depth: float = 1.0,
 52        area_error_th: None | float = None,
 53        lacuna_threshold: None | float = None,
 54    ):
 55        """
 56        :param param: Class param.smashboxparam(), store main smashbox parameters
 57        :type param: param.smashboxparam()
 58        :param query: Any pandas dataframe query as string: '(SURF>20) & (SURF<100)'. This query
 59        must be build using the field (column name) in the outlet database.
 60        https://pandas.pydata.org/docs/user_guide/indexing.html#the-query-method
 61        :type query: str
 62        :max_depth: The maximum depth accepted by the algorithm to find the catchment outlet.
 63            A **max_depth** of 1 means that the algorithm will search among the
 64            combinations in
 65            (``row - 1``, ``row``, ``row + 1``; ``col - 1``, ``col``, ``col + 1``),
 66            the coordinates that minimize
 67            the relative error between the given catchment area and the modeled
 68            catchment area calculated from the
 69            flow directions file.
 70        :type `int`, default 1
 71        :param area_error_th: Tolerance error during the positionning of the outlets. If the Error `(Ssim-Sobs)/Sobs > area_error_th`, the outlet will be excluded.
 72        :type area_error_th: float
 73        :param lacuna_threshold: Lacuna threshold for the discharges in percent. All gauge where the proportion of lacuna exceed this threshold will be removed.
 74        :type: float | Nonetype
 75
 76        """
 77
 78        if param.bbox is not None:
 79            bbox = [
 80                param.bbox["left"],
 81                param.bbox["right"],
 82                param.bbox["bottom"],
 83                param.bbox["top"],
 84            ]
 85        else:
 86            bbox = None
 87
 88        if not os.path.exists(param.outlets_database):
 89            param.outlets_database = os.path.join(
 90                param.asset_dir, "outlets", param.outlets_database
 91            )
 92
 93        if os.path.exists(param.outlets_database):
 94            stations_calage = pd.read_csv(param.outlets_database)
 95        else:
 96            raise ValueError(f"</> Error: file {param.outlets_database} not found")
 97
 98        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
 99        Input_outletsID = param.outletsID
100
101        if param.outletsID is not None:
102            if len(param.outletsID) > 0:
103                stations_calage = stations_calage.loc[
104                    stations_calage[param.outlets_database_fields["id"]].isin(
105                        param.outletsID
106                    )
107                ]
108
109            # check and print outlets not in database
110            for index, id_outlets in enumerate(Input_outletsID):
111
112                if (
113                    id_outlets
114                    not in stations_calage[param.outlets_database_fields["id"]].values
115                ):
116
117                    print(
118                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
119                    )
120
121            if query is not None:
122                stations_calage = stations_calage.query(query)
123
124            if bbox is not None:
125                stations_calage = stations_calage[
126                    (
127                        stations_calage[param.outlets_database_fields["coord_x"]]
128                        >= param.bbox["left"]
129                    )
130                    & (
131                        stations_calage[param.outlets_database_fields["coord_x"]]
132                        <= param.bbox["right"]
133                    )
134                    & (
135                        stations_calage[param.outlets_database_fields["coord_y"]]
136                        >= param.bbox["bottom"]
137                    )
138                    & (
139                        stations_calage[param.outlets_database_fields["coord_y"]]
140                        <= param.bbox["top"]
141                    )
142                ]
143
144            # check and print outlets exluded outside bbox
145            for index, id_outlets in enumerate(Input_outletsID):
146                if (
147                    id_outlets
148                    not in stations_calage[param.outlets_database_fields["id"]].values
149                ):
150                    print(
151                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
152                    )
153
154            if len(stations_calage) == 0:
155                print(
156                    f"</> Error: outlets {param.outletsID} not found in"
157                    "{param.outlets_database}"
158                )
159                raise ValueError(
160                    f"</> Error: outlets {param.outletsID} not found in"
161                    "{param.outlets_database}"
162                )
163
164            columns = {
165                "coord_x": param.outlets_database_fields["coord_x"],
166                "coord_y": param.outlets_database_fields["coord_y"],
167                "area": param.outlets_database_fields["area"],
168                "id": param.outlets_database_fields["id"],
169            }
170
171            # fist build
172            self.mesh = smash.factory.generate_mesh(
173                flwdir_path=param.flowdir,
174                bbox=bbox,
175                x=np.array(stations_calage[columns["coord_x"]][:]),
176                y=np.array(stations_calage[columns["coord_y"]][:]),
177                area=np.array(
178                    stations_calage[columns["area"]][:] * 1e6
179                ),  # Convert km² to m²
180                code=np.array(stations_calage[columns["id"]][:]),
181                epsg=param.epsg,
182                shp_path=param.outlets_shapefile,
183                max_depth=max_depth,
184                area_error_th=area_error_th,
185            )
186
187            if lacuna_threshold is not None:
188                setup = copy.deepcopy(self._setup.setup)
189                setup.update(
190                    {
191                        "read_prcp": False,
192                        "read_pet": False,
193                        "read_qobs": True,
194                        "adjust_interception": False,
195                        "compute_mean_atmos": False,
196                    }
197                )
198
199                # filter gauge if lacuna exceed a thresholds
200                model = smash.Model(setup, self.mesh)
201                qobs = model.response_data.q
202
203                valid_gauge = model.mesh.code[
204                    np.where(
205                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
206                        >= (100.0 - lacuna_threshold)
207                    )
208                ]
209
210                unvalid_gauge = model.mesh.code[
211                    np.where(
212                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
213                        > lacuna_threshold
214                    )
215                ]
216
217                if len(unvalid_gauge) > 0:
218                    print(
219                        f"</> Remove gauges from the mesh where total lacuna "
220                        f"between {setup['start_time']} and {setup['end_time']} "
221                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
222                    )
223
224                stations_calage = stations_calage.loc[
225                    stations_calage[columns["id"]].isin(valid_gauge)
226                ]
227
228                if len(stations_calage) == 0:
229                    print("</> Warnings, no outlets/gauge will be added to the mesh !")
230
231                # rebuild mesh
232                self.mesh = smash.factory.generate_mesh(
233                    flwdir_path=param.flowdir,
234                    bbox=bbox,
235                    x=np.array(stations_calage[columns["coord_x"]][:]),
236                    y=np.array(stations_calage[columns["coord_y"]][:]),
237                    area=np.array(
238                        stations_calage[columns["area"]][:] * 1e6
239                    ),  # Convert km² to m²
240                    code=np.array(stations_calage[columns["id"]][:]),
241                    epsg=param.epsg,
242                    shp_path=param.outlets_shapefile,
243                    max_depth=max_depth,
244                    area_error_th=area_error_th,
245                )
246
247            if param.outlets_shapefile is not None:
248                print("</> Outlets shapefile detected. Loading outlets ...")
249                if param.outlets_database_fields["id_shapefile"] != "None":
250                    col_id = param.outlets_database_fields["id_shapefile"]
251                    code_shape_file = []
252                    for i in range(len(self.mesh["code"])):
253                        sta = self.mesh["code"][i]
254                        code_shape_file.extend(
255                            stations_calage.loc[
256                                stations_calage[columns["id"]] == sta, col_id
257                            ].to_list()
258                        )
259                    code_shape_file = np.array(code_shape_file)
260                else:
261                    code_shape_file = self.mesh["code"]
262
263                catchment_polygon = gpd.read_file(param.outlets_shapefile)
264                self.catchment_polygon = catchment_polygon.loc[
265                    catchment_polygon.code.isin(code_shape_file)
266                ]
267                del catchment_polygon
268
269        else:
270            stations_calage = pd.DataFrame(None)
271
272            if bbox is None:
273                raise ValueError(
274                    "Bbox is None. If no outlets provided, the bbox must be defined."
275                )
276
277            self.mesh = smash.factory.generate_mesh(
278                flwdir_path=param.flowdir, bbox=bbox, epsg=param.epsg, max_depth=max_depth
279            )
Parameters
  • param: Class param.smashboxparam(), store main smashbox parameters
  • query: Any pandas dataframe query as string: '(SURF>20) & (SURF<100)'. This query must be build using the field (column name) in the outlet database. https://pandas.pydata.org/docs/user_guide/indexing.html#the-query-method :max_depth: The maximum depth accepted by the algorithm to find the catchment outlet. A *max_depth* of 1 means that the algorithm will search among the combinations in (row - 1, row, row + 1; col - 1, col, col + 1), the coordinates that minimize the relative error between the given catchment area and the modeled catchment area calculated from the flow directions file. :type int, default 1
  • area_error_th: Tolerance error during the positionning of the outlets. If the Error (Ssim-Sobs)/Sobs > area_error_th, the outlet will be excluded.
  • lacuna_threshold: Lacuna threshold for the discharges in percent. All gauge where the proportion of lacuna exceed this threshold will be removed.