smashbox.model.mesh

View Source

  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(
 96                f"</> Error: file {param.outlets_database} not found"
 97            )
 98
 99        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
100        Input_outletsID = param.outletsID
101
102        if param.outletsID is not None:
103
104            if len(param.outletsID) > 0:
105                stations_calage = (
106                    stations_calage.set_index(
107                        param.outlets_database_fields["id"]
108                    )
109                    .loc[param.outletsID]
110                    .reset_index()
111                )
112
113            # check and print outlets not in database
114            for index, id_outlets in enumerate(Input_outletsID):
115
116                if (
117                    id_outlets
118                    not in stations_calage[
119                        param.outlets_database_fields["id"]
120                    ].values
121                ):
122
123                    print(
124                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
125                    )
126
127            if query is not None:
128                stations_calage = stations_calage.query(query)
129
130            if bbox is not None:
131                stations_calage = stations_calage[
132                    (
133                        stations_calage[
134                            param.outlets_database_fields["coord_x"]
135                        ]
136                        >= param.bbox["left"]
137                    )
138                    & (
139                        stations_calage[
140                            param.outlets_database_fields["coord_x"]
141                        ]
142                        <= param.bbox["right"]
143                    )
144                    & (
145                        stations_calage[
146                            param.outlets_database_fields["coord_y"]
147                        ]
148                        >= param.bbox["bottom"]
149                    )
150                    & (
151                        stations_calage[
152                            param.outlets_database_fields["coord_y"]
153                        ]
154                        <= param.bbox["top"]
155                    )
156                ]
157
158            # check and print outlets exluded outside bbox
159            for index, id_outlets in enumerate(Input_outletsID):
160                if (
161                    id_outlets
162                    not in stations_calage[
163                        param.outlets_database_fields["id"]
164                    ].values
165                ):
166                    print(
167                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
168                    )
169
170            duplicated_stations = stations_calage.duplicated(
171                subset=[param.outlets_database_fields["id"]], keep=False
172            )
173            if len(stations_calage[duplicated_stations]) > 0:
174                print(
175                    "Duplicated stations found in database... Last dupicates are removed:"
176                )
177                print(stations_calage[duplicated_stations])
178                stations_calage = stations_calage.drop_duplicates(
179                    subset=[param.outlets_database_fields["id"]],
180                    keep="first",
181                    ignore_index=True,
182                )
183
184            if len(stations_calage) == 0:
185                print(
186                    f"</> Error: outlets {param.outletsID} not found in"
187                    "{param.outlets_database}"
188                )
189                raise ValueError(
190                    f"</> Error: outlets {param.outletsID} not found in"
191                    "{param.outlets_database}"
192                )
193
194            columns = {
195                "coord_x": param.outlets_database_fields["coord_x"],
196                "coord_y": param.outlets_database_fields["coord_y"],
197                "area": param.outlets_database_fields["area"],
198                "id": param.outlets_database_fields["id"],
199            }
200
201            # first build
202            self.mesh = smash.factory.generate_mesh(
203                flwdir_path=param.flowdir,
204                bbox=bbox,
205                x=np.array(stations_calage[columns["coord_x"]][:]),
206                y=np.array(stations_calage[columns["coord_y"]][:]),
207                area=np.array(
208                    stations_calage[columns["area"]][:] * 1e6
209                ),  # Convert km² to m²
210                code=np.array(stations_calage[columns["id"]][:]),
211                epsg=param.epsg,
212                shp_path=param.outlets_shapefile,
213                max_depth=max_depth,
214                area_error_th=area_error_th,
215            )
216
217            if lacuna_threshold is not None:
218                setup = copy.deepcopy(self._setup.setup)
219                setup.update(
220                    {
221                        "read_prcp": False,
222                        "read_pet": False,
223                        "read_qobs": True,
224                        "adjust_interception": False,
225                        "compute_mean_atmos": False,
226                    }
227                )
228
229                # filter gauge if lacuna exceed a thresholds
230                model = smash.Model(setup, self.mesh)
231                qobs = model.response_data.q
232
233                valid_gauge = model.mesh.code[
234                    np.where(
235                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
236                        >= (100.0 - lacuna_threshold)
237                    )
238                ]
239
240                unvalid_gauge = model.mesh.code[
241                    np.where(
242                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
243                        > lacuna_threshold
244                    )
245                ]
246
247                if len(unvalid_gauge) > 0:
248                    print(
249                        f"</> Remove gauges from the mesh where total lacuna "
250                        f"between {setup['start_time']} and {setup['end_time']} "
251                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
252                    )
253
254                stations_calage = stations_calage.loc[
255                    stations_calage[columns["id"]].isin(valid_gauge)
256                ]
257
258                if len(stations_calage) == 0:
259                    print(
260                        "</> Warnings, no outlets/gauge will be added to the mesh !"
261                    )
262
263                # rebuild mesh
264                self.mesh = smash.factory.generate_mesh(
265                    flwdir_path=param.flowdir,
266                    bbox=bbox,
267                    x=np.array(stations_calage[columns["coord_x"]][:]),
268                    y=np.array(stations_calage[columns["coord_y"]][:]),
269                    area=np.array(
270                        stations_calage[columns["area"]][:] * 1e6
271                    ),  # Convert km² to m²
272                    code=np.array(stations_calage[columns["id"]][:]),
273                    epsg=param.epsg,
274                    shp_path=param.outlets_shapefile,
275                    max_depth=max_depth,
276                    area_error_th=area_error_th,
277                )
278
279            self.load_catchment_polygon(
280                param=param, outlets_db=stations_calage
281            )
282            # ~ if param.outlets_shapefile is not None:
283            # ~ print("</> Outlets shapefile detected. Loading outlets ...")
284            # ~ if param.outlets_database_fields["id_shapefile"] != "None":
285            # ~ col_id = param.outlets_database_fields["id_shapefile"]
286            # ~ code_shape_file = []
287            # ~ for i in range(len(self.mesh["code"])):
288            # ~ sta = self.mesh["code"][i]
289            # ~ code_shape_file.extend(
290            # ~ stations_calage.loc[
291            # ~ stations_calage[columns["id"]] == sta, col_id
292            # ~ ].to_list()
293            # ~ )
294            # ~ code_shape_file = np.array(code_shape_file)
295            # ~ else:
296            # ~ code_shape_file = self.mesh["code"]
297
298            # ~ catchment_polygon = gpd.read_file(param.outlets_shapefile)
299            # ~ self.catchment_polygon = catchment_polygon.loc[
300            # ~ catchment_polygon.code.isin(code_shape_file)
301            # ~ ]
302            # ~ del catchment_polygon
303
304        else:
305            stations_calage = pd.DataFrame(None)
306
307            if bbox is None:
308                raise ValueError(
309                    "Bbox is None. If no outlets provided, the bbox must be defined."
310                )
311
312            self.mesh = smash.factory.generate_mesh(
313                flwdir_path=param.flowdir,
314                bbox=bbox,
315                epsg=param.epsg,
316                max_depth=max_depth,
317            )
318
319    def load_catchment_polygon(
320        self,
321        param: smashboxparam | None = None,
322        outlets_db: pd.DataFrame | None = None,
323    ):
324
325        if param.outlets_shapefile is not None:
326
327            if os.path.exists(param.outlets_database):
328                if outlets_db is None:
329                    outlets_db = pd.read_csv(param.outlets_database)
330            else:
331                raise ValueError(
332                    f"</> Error: file {param.outlets_database} not found"
333                )
334
335            print("</> Outlets shapefile detected. Loading outlets ...")
336            if len(param.outlets_database_fields["id_shapefile"]) > 0:
337                col_id = param.outlets_database_fields["id_shapefile"]
338                code_shape_file = []
339                for i in range(len(self.mesh["code"])):
340                    sta = self.mesh["code"][i]
341                    code_shape_file.extend(
342                        outlets_db.loc[
343                            outlets_db[param.outlets_database_fields["id"]]
344                            == sta,
345                            col_id,
346                        ].to_list()
347                    )
348                code_shape_file = np.array(code_shape_file)
349            else:
350                code_shape_file = self.mesh["code"]
351
352            catchment_polygon = gpd.read_file(param.outlets_shapefile)
353            self.catchment_polygon = catchment_polygon.loc[
354                catchment_polygon.code.isin(code_shape_file)
355            ]
356            del catchment_polygon

class mesh: View Source

 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(
 97                f"</> Error: file {param.outlets_database} not found"
 98            )
 99
100        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
101        Input_outletsID = param.outletsID
102
103        if param.outletsID is not None:
104
105            if len(param.outletsID) > 0:
106                stations_calage = (
107                    stations_calage.set_index(
108                        param.outlets_database_fields["id"]
109                    )
110                    .loc[param.outletsID]
111                    .reset_index()
112                )
113
114            # check and print outlets not in database
115            for index, id_outlets in enumerate(Input_outletsID):
116
117                if (
118                    id_outlets
119                    not in stations_calage[
120                        param.outlets_database_fields["id"]
121                    ].values
122                ):
123
124                    print(
125                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
126                    )
127
128            if query is not None:
129                stations_calage = stations_calage.query(query)
130
131            if bbox is not None:
132                stations_calage = stations_calage[
133                    (
134                        stations_calage[
135                            param.outlets_database_fields["coord_x"]
136                        ]
137                        >= param.bbox["left"]
138                    )
139                    & (
140                        stations_calage[
141                            param.outlets_database_fields["coord_x"]
142                        ]
143                        <= param.bbox["right"]
144                    )
145                    & (
146                        stations_calage[
147                            param.outlets_database_fields["coord_y"]
148                        ]
149                        >= param.bbox["bottom"]
150                    )
151                    & (
152                        stations_calage[
153                            param.outlets_database_fields["coord_y"]
154                        ]
155                        <= param.bbox["top"]
156                    )
157                ]
158
159            # check and print outlets exluded outside bbox
160            for index, id_outlets in enumerate(Input_outletsID):
161                if (
162                    id_outlets
163                    not in stations_calage[
164                        param.outlets_database_fields["id"]
165                    ].values
166                ):
167                    print(
168                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
169                    )
170
171            duplicated_stations = stations_calage.duplicated(
172                subset=[param.outlets_database_fields["id"]], keep=False
173            )
174            if len(stations_calage[duplicated_stations]) > 0:
175                print(
176                    "Duplicated stations found in database... Last dupicates are removed:"
177                )
178                print(stations_calage[duplicated_stations])
179                stations_calage = stations_calage.drop_duplicates(
180                    subset=[param.outlets_database_fields["id"]],
181                    keep="first",
182                    ignore_index=True,
183                )
184
185            if len(stations_calage) == 0:
186                print(
187                    f"</> Error: outlets {param.outletsID} not found in"
188                    "{param.outlets_database}"
189                )
190                raise ValueError(
191                    f"</> Error: outlets {param.outletsID} not found in"
192                    "{param.outlets_database}"
193                )
194
195            columns = {
196                "coord_x": param.outlets_database_fields["coord_x"],
197                "coord_y": param.outlets_database_fields["coord_y"],
198                "area": param.outlets_database_fields["area"],
199                "id": param.outlets_database_fields["id"],
200            }
201
202            # first build
203            self.mesh = smash.factory.generate_mesh(
204                flwdir_path=param.flowdir,
205                bbox=bbox,
206                x=np.array(stations_calage[columns["coord_x"]][:]),
207                y=np.array(stations_calage[columns["coord_y"]][:]),
208                area=np.array(
209                    stations_calage[columns["area"]][:] * 1e6
210                ),  # Convert km² to m²
211                code=np.array(stations_calage[columns["id"]][:]),
212                epsg=param.epsg,
213                shp_path=param.outlets_shapefile,
214                max_depth=max_depth,
215                area_error_th=area_error_th,
216            )
217
218            if lacuna_threshold is not None:
219                setup = copy.deepcopy(self._setup.setup)
220                setup.update(
221                    {
222                        "read_prcp": False,
223                        "read_pet": False,
224                        "read_qobs": True,
225                        "adjust_interception": False,
226                        "compute_mean_atmos": False,
227                    }
228                )
229
230                # filter gauge if lacuna exceed a thresholds
231                model = smash.Model(setup, self.mesh)
232                qobs = model.response_data.q
233
234                valid_gauge = model.mesh.code[
235                    np.where(
236                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
237                        >= (100.0 - lacuna_threshold)
238                    )
239                ]
240
241                unvalid_gauge = model.mesh.code[
242                    np.where(
243                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
244                        > lacuna_threshold
245                    )
246                ]
247
248                if len(unvalid_gauge) > 0:
249                    print(
250                        f"</> Remove gauges from the mesh where total lacuna "
251                        f"between {setup['start_time']} and {setup['end_time']} "
252                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
253                    )
254
255                stations_calage = stations_calage.loc[
256                    stations_calage[columns["id"]].isin(valid_gauge)
257                ]
258
259                if len(stations_calage) == 0:
260                    print(
261                        "</> Warnings, no outlets/gauge will be added to the mesh !"
262                    )
263
264                # rebuild mesh
265                self.mesh = smash.factory.generate_mesh(
266                    flwdir_path=param.flowdir,
267                    bbox=bbox,
268                    x=np.array(stations_calage[columns["coord_x"]][:]),
269                    y=np.array(stations_calage[columns["coord_y"]][:]),
270                    area=np.array(
271                        stations_calage[columns["area"]][:] * 1e6
272                    ),  # Convert km² to m²
273                    code=np.array(stations_calage[columns["id"]][:]),
274                    epsg=param.epsg,
275                    shp_path=param.outlets_shapefile,
276                    max_depth=max_depth,
277                    area_error_th=area_error_th,
278                )
279
280            self.load_catchment_polygon(
281                param=param, outlets_db=stations_calage
282            )
283            # ~ if param.outlets_shapefile is not None:
284            # ~ print("</> Outlets shapefile detected. Loading outlets ...")
285            # ~ if param.outlets_database_fields["id_shapefile"] != "None":
286            # ~ col_id = param.outlets_database_fields["id_shapefile"]
287            # ~ code_shape_file = []
288            # ~ for i in range(len(self.mesh["code"])):
289            # ~ sta = self.mesh["code"][i]
290            # ~ code_shape_file.extend(
291            # ~ stations_calage.loc[
292            # ~ stations_calage[columns["id"]] == sta, col_id
293            # ~ ].to_list()
294            # ~ )
295            # ~ code_shape_file = np.array(code_shape_file)
296            # ~ else:
297            # ~ code_shape_file = self.mesh["code"]
298
299            # ~ catchment_polygon = gpd.read_file(param.outlets_shapefile)
300            # ~ self.catchment_polygon = catchment_polygon.loc[
301            # ~ catchment_polygon.code.isin(code_shape_file)
302            # ~ ]
303            # ~ del catchment_polygon
304
305        else:
306            stations_calage = pd.DataFrame(None)
307
308            if bbox is None:
309                raise ValueError(
310                    "Bbox is None. If no outlets provided, the bbox must be defined."
311                )
312
313            self.mesh = smash.factory.generate_mesh(
314                flwdir_path=param.flowdir,
315                bbox=bbox,
316                epsg=param.epsg,
317                max_depth=max_depth,
318            )
319
320    def load_catchment_polygon(
321        self,
322        param: smashboxparam | None = None,
323        outlets_db: pd.DataFrame | None = None,
324    ):
325
326        if param.outlets_shapefile is not None:
327
328            if os.path.exists(param.outlets_database):
329                if outlets_db is None:
330                    outlets_db = pd.read_csv(param.outlets_database)
331            else:
332                raise ValueError(
333                    f"</> Error: file {param.outlets_database} not found"
334                )
335
336            print("</> Outlets shapefile detected. Loading outlets ...")
337            if len(param.outlets_database_fields["id_shapefile"]) > 0:
338                col_id = param.outlets_database_fields["id_shapefile"]
339                code_shape_file = []
340                for i in range(len(self.mesh["code"])):
341                    sta = self.mesh["code"][i]
342                    code_shape_file.extend(
343                        outlets_db.loc[
344                            outlets_db[param.outlets_database_fields["id"]]
345                            == sta,
346                            col_id,
347                        ].to_list()
348                    )
349                code_shape_file = np.array(code_shape_file)
350            else:
351                code_shape_file = self.mesh["code"]
352
353            catchment_polygon = gpd.read_file(param.outlets_shapefile)
354            self.catchment_polygon = catchment_polygon.loc[
355                catchment_polygon.code.isin(code_shape_file)
356            ]
357            del catchment_polygon

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

mesh(setup) View Source

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): View Source

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): View Source

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): View Source

 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(
 97                f"</> Error: file {param.outlets_database} not found"
 98            )
 99
100        # Pointeur ou copy of param.outletsID ? ici pointeur, ca veux dire que les stations enlevé du mesh sont aussi enlevé de param
101        Input_outletsID = param.outletsID
102
103        if param.outletsID is not None:
104
105            if len(param.outletsID) > 0:
106                stations_calage = (
107                    stations_calage.set_index(
108                        param.outlets_database_fields["id"]
109                    )
110                    .loc[param.outletsID]
111                    .reset_index()
112                )
113
114            # check and print outlets not in database
115            for index, id_outlets in enumerate(Input_outletsID):
116
117                if (
118                    id_outlets
119                    not in stations_calage[
120                        param.outlets_database_fields["id"]
121                    ].values
122                ):
123
124                    print(
125                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` does not exist in the outlets database. This outlets will not be included in the mesh..."
126                    )
127
128            if query is not None:
129                stations_calage = stations_calage.query(query)
130
131            if bbox is not None:
132                stations_calage = stations_calage[
133                    (
134                        stations_calage[
135                            param.outlets_database_fields["coord_x"]
136                        ]
137                        >= param.bbox["left"]
138                    )
139                    & (
140                        stations_calage[
141                            param.outlets_database_fields["coord_x"]
142                        ]
143                        <= param.bbox["right"]
144                    )
145                    & (
146                        stations_calage[
147                            param.outlets_database_fields["coord_y"]
148                        ]
149                        >= param.bbox["bottom"]
150                    )
151                    & (
152                        stations_calage[
153                            param.outlets_database_fields["coord_y"]
154                        ]
155                        <= param.bbox["top"]
156                    )
157                ]
158
159            # check and print outlets exluded outside bbox
160            for index, id_outlets in enumerate(Input_outletsID):
161                if (
162                    id_outlets
163                    not in stations_calage[
164                        param.outlets_database_fields["id"]
165                    ].values
166                ):
167                    print(
168                        f"</> Warning: Outlets `{Input_outletsID.pop(index)}` outside the boundingbox {bbox}. This outlets will not be included in the mesh..."
169                    )
170
171            duplicated_stations = stations_calage.duplicated(
172                subset=[param.outlets_database_fields["id"]], keep=False
173            )
174            if len(stations_calage[duplicated_stations]) > 0:
175                print(
176                    "Duplicated stations found in database... Last dupicates are removed:"
177                )
178                print(stations_calage[duplicated_stations])
179                stations_calage = stations_calage.drop_duplicates(
180                    subset=[param.outlets_database_fields["id"]],
181                    keep="first",
182                    ignore_index=True,
183                )
184
185            if len(stations_calage) == 0:
186                print(
187                    f"</> Error: outlets {param.outletsID} not found in"
188                    "{param.outlets_database}"
189                )
190                raise ValueError(
191                    f"</> Error: outlets {param.outletsID} not found in"
192                    "{param.outlets_database}"
193                )
194
195            columns = {
196                "coord_x": param.outlets_database_fields["coord_x"],
197                "coord_y": param.outlets_database_fields["coord_y"],
198                "area": param.outlets_database_fields["area"],
199                "id": param.outlets_database_fields["id"],
200            }
201
202            # first build
203            self.mesh = smash.factory.generate_mesh(
204                flwdir_path=param.flowdir,
205                bbox=bbox,
206                x=np.array(stations_calage[columns["coord_x"]][:]),
207                y=np.array(stations_calage[columns["coord_y"]][:]),
208                area=np.array(
209                    stations_calage[columns["area"]][:] * 1e6
210                ),  # Convert km² to m²
211                code=np.array(stations_calage[columns["id"]][:]),
212                epsg=param.epsg,
213                shp_path=param.outlets_shapefile,
214                max_depth=max_depth,
215                area_error_th=area_error_th,
216            )
217
218            if lacuna_threshold is not None:
219                setup = copy.deepcopy(self._setup.setup)
220                setup.update(
221                    {
222                        "read_prcp": False,
223                        "read_pet": False,
224                        "read_qobs": True,
225                        "adjust_interception": False,
226                        "compute_mean_atmos": False,
227                    }
228                )
229
230                # filter gauge if lacuna exceed a thresholds
231                model = smash.Model(setup, self.mesh)
232                qobs = model.response_data.q
233
234                valid_gauge = model.mesh.code[
235                    np.where(
236                        np.sum(qobs > 0.0, axis=1) / qobs.shape[1] * 100.0
237                        >= (100.0 - lacuna_threshold)
238                    )
239                ]
240
241                unvalid_gauge = model.mesh.code[
242                    np.where(
243                        np.sum(qobs < 0.0, axis=1) / qobs.shape[1] * 100.0
244                        > lacuna_threshold
245                    )
246                ]
247
248                if len(unvalid_gauge) > 0:
249                    print(
250                        f"</> Remove gauges from the mesh where total lacuna "
251                        f"between {setup['start_time']} and {setup['end_time']} "
252                        f"exceed {lacuna_threshold}%: {unvalid_gauge}"
253                    )
254
255                stations_calage = stations_calage.loc[
256                    stations_calage[columns["id"]].isin(valid_gauge)
257                ]
258
259                if len(stations_calage) == 0:
260                    print(
261                        "</> Warnings, no outlets/gauge will be added to the mesh !"
262                    )
263
264                # rebuild mesh
265                self.mesh = smash.factory.generate_mesh(
266                    flwdir_path=param.flowdir,
267                    bbox=bbox,
268                    x=np.array(stations_calage[columns["coord_x"]][:]),
269                    y=np.array(stations_calage[columns["coord_y"]][:]),
270                    area=np.array(
271                        stations_calage[columns["area"]][:] * 1e6
272                    ),  # Convert km² to m²
273                    code=np.array(stations_calage[columns["id"]][:]),
274                    epsg=param.epsg,
275                    shp_path=param.outlets_shapefile,
276                    max_depth=max_depth,
277                    area_error_th=area_error_th,
278                )
279
280            self.load_catchment_polygon(
281                param=param, outlets_db=stations_calage
282            )
283            # ~ if param.outlets_shapefile is not None:
284            # ~ print("</> Outlets shapefile detected. Loading outlets ...")
285            # ~ if param.outlets_database_fields["id_shapefile"] != "None":
286            # ~ col_id = param.outlets_database_fields["id_shapefile"]
287            # ~ code_shape_file = []
288            # ~ for i in range(len(self.mesh["code"])):
289            # ~ sta = self.mesh["code"][i]
290            # ~ code_shape_file.extend(
291            # ~ stations_calage.loc[
292            # ~ stations_calage[columns["id"]] == sta, col_id
293            # ~ ].to_list()
294            # ~ )
295            # ~ code_shape_file = np.array(code_shape_file)
296            # ~ else:
297            # ~ code_shape_file = self.mesh["code"]
298
299            # ~ catchment_polygon = gpd.read_file(param.outlets_shapefile)
300            # ~ self.catchment_polygon = catchment_polygon.loc[
301            # ~ catchment_polygon.code.isin(code_shape_file)
302            # ~ ]
303            # ~ del catchment_polygon
304
305        else:
306            stations_calage = pd.DataFrame(None)
307
308            if bbox is None:
309                raise ValueError(
310                    "Bbox is None. If no outlets provided, the bbox must be defined."
311                )
312
313            self.mesh = smash.factory.generate_mesh(
314                flwdir_path=param.flowdir,
315                bbox=bbox,
316                epsg=param.epsg,
317                max_depth=max_depth,
318            )

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.

def load_catchment_polygon( self, param: smashbox.init.param.smashboxparam | None = None, outlets_db: pandas.DataFrame | None = None): View Source

320    def load_catchment_polygon(
321        self,
322        param: smashboxparam | None = None,
323        outlets_db: pd.DataFrame | None = None,
324    ):
325
326        if param.outlets_shapefile is not None:
327
328            if os.path.exists(param.outlets_database):
329                if outlets_db is None:
330                    outlets_db = pd.read_csv(param.outlets_database)
331            else:
332                raise ValueError(
333                    f"</> Error: file {param.outlets_database} not found"
334                )
335
336            print("</> Outlets shapefile detected. Loading outlets ...")
337            if len(param.outlets_database_fields["id_shapefile"]) > 0:
338                col_id = param.outlets_database_fields["id_shapefile"]
339                code_shape_file = []
340                for i in range(len(self.mesh["code"])):
341                    sta = self.mesh["code"][i]
342                    code_shape_file.extend(
343                        outlets_db.loc[
344                            outlets_db[param.outlets_database_fields["id"]]
345                            == sta,
346                            col_id,
347                        ].to_list()
348                    )
349                code_shape_file = np.array(code_shape_file)
350            else:
351                code_shape_file = self.mesh["code"]
352
353            catchment_polygon = gpd.read_file(param.outlets_shapefile)
354            self.catchment_polygon = catchment_polygon.loc[
355                catchment_polygon.code.isin(code_shape_file)
356            ]
357            del catchment_polygon