smashbox.tutorial_readme

  1# -*- coding: utf-8 -*-
  2
  3if __name__ == "__main__":
  4
  5    import smashbox
  6    import numpy as np
  7
  8    # import pyhdf5_handler
  9    # import smash
 10
 11    # setup,mesh=smash.factory.load_dataset('Cance')
 12    # model=smash.Model(setup,mesh)
 13
 14    bbox = {
 15        "left": 925000.0,
 16        "bottom": 6230000.0,
 17        "right": 974000.0,
 18        "top": 6255000.0,
 19    }
 20
 21    sb = smashbox.SmashBox()
 22
 23    sb.myparam.set_param("bbox", bbox)
 24    sb.myparam.set_param(
 25        "smash_parameters", "/home/maxime/DEV/smashbox/smashbox/asset/params"
 26    )
 27    sb.myparam.set_param(
 28        "outlets_database",
 29        "/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_sites.csv",
 30    )
 31
 32    sb.newmodel("graffas_zone")
 33    sb.graffas_zone.mysetup.load_setup("setup_rhax_gr4_dt3600")
 34
 35    sb.graffas_zone.mysetup.update_setup(
 36        {
 37            "pet_directory": "/home/maxime/DATA/ETP-SFR-FRA-INTERA_L93",
 38            "prcp_directory": "/home/maxime/DATA/PLUIE",
 39            "qobs_directory": "/home/maxime/DATA/QOBS_SITE_60M",
 40        }
 41    )  # Change value of the Smash setup
 42
 43    sb.graffas_zone.generate_mesh(
 44        query="(SURF>20) & (INFLUENCE=='Influence nulle ou faible')"
 45    )  # direct method to build the mesh (recommended)
 46
 47    sb.graffas_zone.myplot.plot_catchment_surface_consistency(
 48        fig_settings={
 49            "figname": "../images/mesh_surface_consistency.png",
 50            "xsize": 4,
 51            "ysize": 4,
 52            "font_ratio": 1,
 53        },
 54        ax_settings={"title_fontsize": 12},
 55    )
 56    sb.graffas_zone.myplot.plot_catchment_surface_error(
 57        fig_settings={
 58            "figname": "../images/mesh_surface_error.png",
 59            "xsize": 4,
 60            "ysize": 4,
 61            "font_ratio": 1,
 62        },
 63        ax_settings={"title_fontsize": 12},
 64    )
 65
 66    sb.graffas_zone.myplot.plot_mesh(
 67        fig_settings={
 68            "figname": "../images/mesh.png",
 69            "xsize": 6,
 70            "ysize": 4,
 71        },
 72        ax_settings={"title_fontsize": 14},
 73    )
 74
 75    nrow = sb.graffas_zone.mymesh.mesh["nrow"]
 76    ncol = sb.graffas_zone.mymesh.mesh["ncol"]
 77    chunk_size = 18
 78    nb_chunk = 21
 79    graffas_prcp = np.zeros(shape=(nrow, ncol, nb_chunk * chunk_size))
 80    fctr = 100
 81    for i in range(nb_chunk):
 82        graffas_prcp[:, :, i * chunk_size : (i + 1) * chunk_size] = (
 83            fctr
 84            * np.arange(i * chunk_size, (i + 1) * chunk_size)
 85            / (chunk_size * nb_chunk)
 86            * np.cos(np.arange(0, chunk_size * 10, 10) * np.pi / 180)
 87        )
 88    graffas_prcp = np.where(graffas_prcp < 0, 0.0, graffas_prcp)
 89
 90    sb.graffas_zone.model_warmup(warmup=365)
 91
 92    sb.graffas_zone.atmos_data_connector(
 93        input_prcp=graffas_prcp, input_dt=3600.0
 94    )
 95
 96    sb.graffas_zone.model()
 97
 98    sb.graffas_zone.myplot.target = "mysmashmodel"
 99    sb.graffas_zone.myplot.plot_parameters(
100        param="cp",
101        fig_settings={
102            "figname": "../images/param_cp.png",
103            "xsize": 5,
104            "ysize": 5,
105        },
106        ax_settings={"font_ratio": 1},
107    )
108
109    sb.graffas_zone.myplot.multiplot_parameters(
110        fig_settings={
111            "figname": "../images/multiplot_param.png",
112            "font_ratio": 1,
113            "xsize": 7,
114            "ysize": 6,
115        },
116        ax_settings={"title_fontsize": 10, "font_ratio": 0.6},
117    )
118
119    sb.graffas_zone.forward_run(
120        invert_states=True, return_options={"q_domain": True}
121    )
122
123    rand_q = np.random.rand(
124        *sb.graffas_zone.mysmashmodel.smash.response_data.q.shape
125    )
126    sb.graffas_zone.mysmashmodel.smash.response_data.q[
127        :
128    ] = sb.graffas_zone.mysmashmodel.smash.response.q[
129        :
130    ] + sb.graffas_zone.mysmashmodel.smash.response.q[
131        :
132    ] * (
133        rand_q[:] - 0.5
134    ) * np.repeat(
135        np.arange(rand_q.shape[0]).reshape((rand_q.shape[0], 1)),
136        rand_q.shape[1],
137        axis=1,
138    )
139
140    sb.graffas_zone.mysmashmodel.mystats.fspatial_stats()
141    sb.graffas_zone.mysmashmodel.mystats.foutlets_stats()
142    sb.graffas_zone.mysmashmodel.mystats.fmisfit_stats(use_smash_metrics=False)
143    sb.graffas_zone.mysmashmodel.mystats.fmisfit_stats(use_smash_metrics=True)
144
145    # test des metrics
146    # -------------------------------------------------------
147    from smash.fcore import _mwd_metrics as smash_metrics
148
149    metric_list = [
150        "nse",
151        "kge",
152        "lgrm",
153        "mae",
154        "mape",
155        "mse",
156        "nnse",
157        "rmse",
158        "se",
159    ]
160
161    for metric_name in metric_list:
162        # metric_name = "nse"
163        fmetric = getattr(smash_metrics, metric_name)
164        metric = (
165            np.zeros(shape=(len(sb.graffas_zone.mysmashmodel.smash.mesh.code)))
166            + np.nan
167        )
168        for i in range(len(sb.graffas_zone.mysmashmodel.smash.mesh.code)):
169            qobs = sb.graffas_zone.mysmashmodel.smash.response_data.q[i, :]
170            qsim = sb.graffas_zone.mysmashmodel.smash.response.q[i, :]
171            metric[i] = fmetric(qobs, qsim)
172
173        if hasattr(
174            sb.graffas_zone.mysmashmodel.mystats.misfit_stats.results,
175            metric_name,
176        ):
177            smashbox_metric = getattr(
178                sb.graffas_zone.mysmashmodel.mystats.misfit_stats.results,
179                metric_name,
180            )
181            print(f"diff on metric {metric_name}:")
182            print(smashbox_metric - metric)
183        else:
184            print(f"no metric {metric_name} found")
185
186    # -------------------------------------------------------
187
188    sb.graffas_zone.mysmashmodel.mystats.fquantile_stats(
189        chunk_size=3,
190        estimate_method="MLE",
191        ncpu=4,
192        fit="gumbel",
193        compute_uncertainties=False,
194    )
195
196    sb.graffas_zone.myplot.plot_outlets_quantile(
197        fig_settings={
198            "figname": "../images/quantile.png",
199            "xsize": 10,
200            "ysize": 10,
201        },
202    )
203    sb.graffas_zone.myplot.plot_outlets_quantile(quantile_obs=False)
204
205    sb.graffas_zone.mysmashmodel.mystats.fmaxima_stats(
206        chunk_size=3, cumulated_maxima=True
207    )
208
209    # import matplotlib.pyplot as plt
210
211    # var = "spatial_quantile_outlets"
212
213    # fig, ax = plt.subplots()
214
215    # var_outlets = getattr(es.graffas_zone.mystats.quantile_stats, var)
216
217    # for i in range(var_outlets.shape[0]):
218    #     ax.plot(var_outlets[i, :], "x", markersize=5)
219
220    # for i in range(len(es.graffas_zone.mymesh.mesh["code"])):
221    #     coords = sb.graffas_zone.mymesh.mesh["gauge_pos"][i]
222    #     var_outlets = getattr(
223    #         sb.graffas_zone.mystats.quantile_stats, var.removesuffix("_outlets")
224    #     )[coords[0], coords[1], :]
225
226    #     ax.plot(var_outlets[:, 1], "o", markersize=4)
227
228    # save the smashbox container:
229    sb.graffas_zone.save_model_container(path="../smashbox_saved")
230
231    import smashbox
232
233    sb_in = smashbox.SmashBox()
234    sb_in.load_containers(path="../smashbox_saved")
235
236    dict_es = sb.graffas_zone.save_model_container_hdf5()
237
238    sb.graffas_zone.myplot.plot_spatial_stats(
239        fig_settings={
240            "figname": "../images/spatial_stats_max.png",
241            "xsize": 5,
242            "ysize": 5,
243        },
244        ax_settings={"font_ratio": 0.8},
245    )
246
247    sb.graffas_zone.myplot.plot_xy_quantile(
248        X=20,
249        Y=25,
250        fig_settings={
251            "figname": "../images/quantile_XY.png",
252            "xsize": 4,
253            "ysize": 4,
254        },
255        ax_settings={"font_ratio": 1},
256    )
257
258    sb.graffas_zone.myplot.plot_outlets_quantile(
259        fig_settings={
260            "figname": "../images/quantile_outlets.png",
261            "xsize": 6,
262            "ysize": 8,
263        },
264        ax_settings={"font_ratio": 0.8},
265    )
266
267    sb.graffas_zone.myplot.plot_spatial_quantile(
268        T=2,
269        duration=1,
270        fig_settings={
271            "figname": "../images/spatial_quantile.png",
272            "xsize": 4,
273            "ysize": 4,
274        },
275        ax_settings={"font_ratio": 0.8},
276    )
277
278    sb.graffas_zone.myplot.multiplot_spatial_quantile(
279        duration=1,
280        fig_settings={
281            "figname": "../images/multiplot_spatial_quantile.png",
282            "xsize": 6,
283            "ysize": 6,
284        },
285        ax_settings={"font_ratio": 0.6},
286    )
287
288    sb.graffas_zone.myplot.plot_hydrograph(
289        fig_settings={
290            "figname": "../images/hydrogram.png",
291        },
292    )
293
294    sb.graffas_zone.myplot.plot_misfit(
295        misfit="nse",
296        fig_settings={
297            "figname": "../images/nse2.png",
298            "xsize": 4,
299            "ysize": 3,
300        },
301    )
302
303    sb.graffas_zone.myplot.multiplot_misfit(
304        fig_settings={
305            "figname": "../images/multiplot_misfit2.png",
306            "xsize": 6,
307            "ysize": 6,
308        },
309        ax_settings={"font_ratio": 0.4},
310    )
311
312    sb.graffas_zone.myplot.plot_misfit_map(
313        misfit="nse",
314        fig_settings={
315            "figname": "../images/misfit_map.png",
316            "xsize": 8,
317            "ysize": 6,
318        },
319        ax_settings={"font_ratio": 1},
320    )