1if __name__ == "__main__":
2
3 import smashbox
4
5 from smashbox import tools
6
7 # @tools.tools.autocast_args
8 # def test(
9 # a: int | None = None,
10 # b: None | int = None,
11 # c: int | None = 1,
12 # d: int = 2,
13 # e: None = None,
14 # f: None | dict | tuple = None,
15 # g=1,
16 # ):
17
18 # print(a, b, c, d, e)
19
20 # test(a=1.0, b=1.0, c=1.0, d=1.0, e=None, f=(2,), g=2)
21
22 # Main object creation
23 sb = smashbox.SmashBox()
24
25 # Main input parameters of main object: class param, attribute 'myparam'
26 bbox = {"left": 925000.0, "bottom": 6230000.0, "right": 974000.0, "top": 6255000.0}
27 # sb.currentmodel #current active model
28 sb.myparam.list_param() # list the parameters of the current active model
29 sb.myparam.set_param(
30 "outlets_database",
31 "/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_stations_example.csv",
32 )
33 sb.myparam.set_param("outletsID", ["Y4604020", "Y4624010", "Y4002010", "Y4105210"])
34 sb.myparam.set_param("bbox", bbox)
35
36 sb.myparam.set_param(
37 "smash_parameters", "/home/maxime/DEV/smashbox/smashbox/asset/params"
38 ) # change the value of myparam used byRealCollobrier
39
40 # New model creation
41 sb.newmodel(
42 "RealCollobrier"
43 ) # create a new model named RealCollobrier and copy the previous parameters used for mymodel1 (default behaviour, but it can turned off)
44
45 # create an second model with different parameters
46 sb.myparam.param.setup_file = "setup_rhax_gr5_dt3600.yaml"
47
48 # Param associated to the model (hidden) has been copied
49 sb.RealCollobrier._myparam.list_param() # Notice that myparam is also stored in each model object. So you can access to the model param inside the model attribute.
50
51 # Setup utilities
52 sb.RealCollobrier.mysetup.setup # print the setup to the screen
53 sb.RealCollobrier.mysetup.update_setup(
54 {
55 "hydrological_module": "gr4",
56 "end_time": "2014-02-01 12:00",
57 "start_time": "2014-01-01 01:00",
58 "pet_directory": "/home/maxime/DATA/ETP-SFR-FRA-INTERA_L93",
59 "prcp_directory": "/home/maxime/DATA/PLUIE",
60 "qobs_directory": "/home/maxime/DATA/QOBS_60min",
61 }
62 ) # Change value of the Smash setup
63 sb.RealCollobrier.mysetup.setup # print the setup to the screen
64
65 sb.RealCollobrier.generate_mesh() # direct method to build the mesh (recommended)
66
67 sb.RealCollobrier.model()
68
69 sb.RealCollobrier.mysmashmodel.set_rr_parameters("cp", 500.0)
70 sb.RealCollobrier.mysmashmodel.set_rr_parameters("ct", 200.0)
71 sb.RealCollobrier.mysmashmodel.set_rr_parameters("llr", 50.0)
72 sb.RealCollobrier.mysmashmodel.set_rr_parameters("kexc", 0.0)
73
74 sb.RealCollobrier.optimize(
75 cost_options={"end_warmup": "2014-01-10 01:00", "gauge": "Y4604020"},
76 mapping="uniform",
77 optimizer="sbs",
78 optimize_options={
79 "bounds": {
80 "cp": [1, 1000],
81 "ct": [1, 1000],
82 "kexc": [-50, 50],
83 "llr": [1, 1000],
84 },
85 "termination_crit": {"maxiter": 10},
86 },
87 )
88
89 sb.RealCollobrier.optimize(
90 cost_options={
91 "end_warmup": "2014-01-10 01:00",
92 "gauge": ["Y4604020", "Y4624010"],
93 },
94 mapping="distributed",
95 optimizer="lbfgsb",
96 optimize_options={
97 "bounds": {
98 "cp": [1, 1000],
99 "ct": [1, 1000],
100 "kexc": [-50, 50],
101 "llr": [1, 1000],
102 },
103 "termination_crit": {"maxiter": 10},
104 },
105 )
106
107 # Stat utilities
108
109 # plot utilities
110 sb.RealCollobrier.myplot.plot_mesh()
111
112 import smashbox
113 import smash
114 from smashbox.plot import plot
115
116 # Testing graffas connector
117 import numpy as np
118
119 bbox = {"left": 925000.0, "bottom": 6230000.0, "right": 974000.0, "top": 6255000.0}
120
121 # graffas_prcp = np.zeros(shape=(145, 170, 1300))
122 # rng = np.random.default_rng()
123 # graffas_prcp = (
124 # graffas_prcp
125 # + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
126 # )
127 # graffas_pet = np.zeros(shape=(142, 166, 17520)) + 1.0
128
129 # Main object creation
130 sb = smashbox.SmashBox()
131
132 sb.myparam.set_param("bbox", bbox)
133 sb.myparam.set_param(
134 "smash_parameters", "/home/maxime/DEV/smashbox/smashbox/asset/params"
135 )
136 sb.myparam.set_param(
137 "outlets_database",
138 "/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_sites.csv",
139 )
140 # sb.myparam.set_param(
141 # "outlets_shapefile", "/home/maxime/DATA/BNBVlight/BNBV_light_bassins.shp"
142 # )
143 # sb.myparam.set_param("enhanced_smash_input_data", True)
144
145 sb.newmodel("graffas_zone")
146 sb.graffas_zone.mysetup.load_setup("setup_rhax_gr4_dt3600")
147
148 sb.graffas_zone.mysetup.update_setup(
149 {
150 "pet_directory": "/home/maxime/DATA/ETP-SFR-FRA-INTERA_L93",
151 "prcp_directory": "/home/maxime/DATA/PLUIE",
152 "qobs_directory": "/home/maxime/DATA/QOBS_60min",
153 }
154 ) # Change value of the Smash setup
155
156 sb.graffas_zone.generate_mesh(
157 lacuna_threshold=20.0
158 ) # direct method to build the mesh (recommended)
159
160 sb.graffas_zone.myplot.plot_catchment_surface_consistency(
161 fig_settings={
162 "figname": "../images/mesh_surface_consistency.png",
163 "xsize": 5,
164 "ysize": 5,
165 "font_ratio": 2,
166 },
167 ax_settings={"title_fontsize": 12},
168 )
169 sb.graffas_zone.myplot.plot_catchment_surface_error(
170 fig_settings={
171 "figname": "../images/mesh_surface_error.png",
172 "xsize": 5,
173 "ysize": 5,
174 "font_ratio": 2,
175 },
176 ax_settings={"title_fontsize": 12},
177 )
178
179 sb.graffas_zone.myplot.plot_mesh(
180 fig_settings={
181 "figname": "../images/mesh.png",
182 "xsize": 10,
183 "ysize": 10,
184 },
185 ax_settings={"title_fontsize": 20},
186 )
187
188 nrow = sb.graffas_zone.mymesh.mesh["nrow"]
189 ncol = sb.graffas_zone.mymesh.mesh["ncol"]
190 chunk_size = 18
191 nb_chunk = 21
192 graffas_prcp = np.zeros(shape=(nrow, ncol, nb_chunk * chunk_size))
193 fctr = 100
194 for i in range(nb_chunk):
195 graffas_prcp[:, :, i * chunk_size : (i + 1) * chunk_size] = (
196 fctr
197 * np.arange(i * chunk_size, (i + 1) * chunk_size)
198 / (chunk_size * nb_chunk)
199 * np.cos(np.arange(0, chunk_size * 10, 10) * np.pi / 180)
200 )
201 graffas_prcp = np.where(graffas_prcp < 0, 0.0, graffas_prcp)
202
203 sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp)
204 # sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp, input_pet=graffas_pet)
205
206 sb.graffas_zone.model()
207
208 sb.graffas_zone.myplot.plot_parameters(
209 param="cp",
210 fig_settings={
211 "figname": "../images/param_cp.png",
212 "xsize": 10,
213 "ysize": 10,
214 },
215 ax_settings={"title_fontsize": 12},
216 )
217 sb.graffas_zone.myplot.multiplot_parameters(
218 fig_settings={
219 "figname": "../images/multiplot_param.png",
220 "xsize": 10,
221 "ysize": 10,
222 "font_ratio": 1,
223 },
224 ax_settings={"title_fontsize": 10, "font_ratio": 1},
225 )
226
227 # sb.graffas_zone.generate_model(
228 # return_options={"q_domain": True}
229 # ) # build and run the model with smash.Model()
230
231 sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
232
233 # sb.graffas_zone.myplot.plot_parameters(param="cp")
234 # sb.graffas_zone.myplot.multiplot_parameters()
235 # import matplotlib.pyplot as plt
236 # from matplotlib.colors import LogNorm, SymLogNorm
237
238 # q = sb.graffas_zone.extra_smash_results.q_domain[:, :, 500]
239 # q = np.where(es.graffas_zone.mysmashmodel.mesh.active_cell == 0, np.nan, q)
240 # plt.imshow(q, norm=SymLogNorm(1e-4))
241
242 rand_q = np.random.rand(*sb.graffas_zone.mysmashmodel.response_data.q.shape)
243 sb.graffas_zone.mysmashmodel.response_data.q[
244 :
245 ] = sb.graffas_zone.mysmashmodel.response.q[
246 :
247 ] + sb.graffas_zone.mysmashmodel.response.q[
248 :
249 ] * (
250 rand_q[:] - 0.5
251 ) * np.repeat(
252 np.arange(rand_q.shape[0]).reshape((rand_q.shape[0], 1)), rand_q.shape[1], axis=1
253 )
254
255 sb.graffas_zone.mystats.fspatial_stats()
256 sb.graffas_zone.mystats.foutlets_stats()
257 sb.graffas_zone.mystats.fmisfit_stats()
258
259 sb.graffas_zone.mystats.fquantile_stats(
260 chunk_size=3,
261 estimate_method="MLE",
262 ncpu=6,
263 fit="gumbel",
264 compute_uncertainties=True,
265 )
266
267 # sb.graffas_zone.mystats.misfit.rmse()
268 # sb.graffas_zone.mystats.misfit.nse(outlets_name=["MedEst_528", "MedEst_617"])
269 # res = sb.graffas_zone.mystats.stats_misfit()
270
271 from smashbox.stats import stats
272
273 # results = stats.spatial_quantiles(
274 # sb.graffas_zone.extra_smash_results.q_domain,
275 # chunk_size=3,
276 # quantile_duration=48,
277 # estimate_method="MLE",
278 # )
279
280 # results = stats.spatial_quantiles(
281 # sb.graffas_zone.extra_smash_results.q_domain,
282 # chunk_size=3,
283 # quantile_duration=72,
284 # estimate_method="MLE",
285 # )
286
287 # array = stats.time_resample_array(
288 # array=es.graffas_zone.extra_smash_results.q_domain,
289 # quantile_duration=1,
290 # model_time_step=3600,
291 # t_axis=2,
292 # )
293 # array.shape
294
295 # stats.compute_maxima(
296 # array=array,
297 # t_axis=2,
298 # nb_minimum_chunks=4,
299 # chunk_size=2,
300 # quantile_duration=48,
301 # )
302
303 # array = stats.time_resample_array(
304 # array=es.graffas_zone.extra_smash_results.q_domain,
305 # quantile_duration=48,
306 # model_time_step=3600,
307 # t_axis=2,
308 # )
309 # array.shape
310
311 # maxima = stats.compute_maxima(
312 # array=array,
313 # t_axis=2,
314 # nb_minimum_chunks=5,
315 # chunk_size=2,
316 # quantile_duration=48,
317 # )
318
319 # sb.graffas_zone.mystats.fmaxima_stats(chunk_size=4, cumulated_maxima=False)
320
321 # sb.graffas_zone.mystats.fstats_quantile_2(
322 # chunk_size=10, estimate_method="MLE", ncpu=6, fit="gumbel"
323 # )
324
325 # sb.graffas_zone.mystats.fmaxima_stats(chunk_size=2, cumulated_maxima=True)
326
327 # sb.graffas_zone.mystats.quantile_stats.spatial_cumulated_maxima[100, 10, :, 0]
328
329 sb.graffas_zone.myplot.plot_spatial_stats(
330 fig_settings={
331 "figname": "../images/spatial_stats_max.png",
332 "xsize": 10,
333 "ysize": 10,
334 },
335 ax_settings={"font_ratio": 1.0},
336 )
337
338 # results = stats.fit_quantile(
339 # maxima=es.graffas_zone.mystats.quantile_stats.spatial_cumulated_maxima[
340 # :, :, :, 0
341 # ],
342 # t_axis=2,
343 # return_periods=[2, 5, 10, 20, 50, 100],
344 # fit="gumbel",
345 # estimate_method="MLE",
346 # quantile_duration=1,
347 # ncpu=6,
348 # )
349 # results = stats.parametric_bootstrap_uncertainties(q_results=results, ncpu=6)
350
351 # fig, ax = smashbox.plot.plot.plot_xy_quantile(
352 # results,
353 # X=100,
354 # Y=10,
355 # fig_settings={"figname": "../images/quantile_XY.png", "font_ratio": 2},
356 # ax_settings={"legend_fontsize": 20},
357 # )
358
359 sb.graffas_zone.myplot.plot_xy_quantile(
360 X=20,
361 Y=25,
362 fig_settings={"figname": "../images/quantile_XY.png", "font_ratio": 2},
363 ax_settings={"legend_fontsize": 20},
364 )
365
366 sb.graffas_zone.myplot.plot_outlets_quantile(
367 fig_settings={
368 "figname": "../images/quantile_outlets.png",
369 "xsize": 10,
370 "ysize": 10,
371 },
372 ax_settings={"font_ratio": 0.5},
373 )
374
375 sb.graffas_zone.myplot.plot_spatial_quantile(
376 T=2,
377 duration=1,
378 fig_settings={
379 "figname": "../images/spatial_quantile.png",
380 "xsize": 10,
381 "ysize": 10,
382 },
383 ax_settings={"font_ratio": 1.5},
384 )
385
386 sb.graffas_zone.myplot.multiplot_spatial_quantile(
387 duration=1,
388 fig_settings={
389 "figname": "../images/multiplot_spatial_quantile.png",
390 "xsize": 10,
391 "ysize": 10,
392 },
393 ax_settings={"font_ratio": 1.0},
394 )
395
396 # fig, ax = plot.plot_chro(
397 # sb.graffas_zone.mysmashmodel.response.q,
398 # ax_settings={"title": "my title", "title_fontsize": 42},
399 # plot_settings={"label": "discharge at gauge 0"},
400 # )
401
402 sb.graffas_zone.myplot.plot_hydrograph(
403 fig_settings={
404 "figname": "../images/hydrogram.png",
405 "xsize": 10,
406 "ysize": 10,
407 },
408 )
409
410 sb.graffas_zone.myplot.plot_misfit(
411 misfit="nse",
412 fig_settings={
413 "figname": "../images/nse.png",
414 "xsize": 10,
415 "ysize": 7,
416 },
417 )
418
419 sb.graffas_zone.myplot.multiplot_misfit(
420 fig_settings={
421 "figname": "../images/multiplot_misfit.png",
422 "xsize": 10,
423 "ysize": 10,
424 },
425 )
426
427 sb.graffas_zone.myplot.plot_misfit_map(
428 misfit="nse",
429 fig_settings={"figname": "../images/misfit_map.png", "xsize": 10, "ysize": 10},
430 ax_settings={"font_ratio": 1},
431 )
432
433 # sb.graffas_zone.export_parameters()
434
435 import smashbox
436 import smash
437 from smashbox.plot import plot
438
439 # Testing graffas connector
440 import numpy as np
441
442 bbox = {"left": 925000.0, "bottom": 6230000.0, "right": 974000.0, "top": 6255000.0}
443
444 sb = smashbox.SmashBox()
445
446 sb.myparam.set_param("bbox", bbox)
447 sb.myparam.set_param("outlets_database", "db_sites")
448 sb.myparam.set_param(
449 "outlets_shapefile",
450 "/home/maxime/DEV/smashbox/smashbox/asset/shapefile/Hydro_Bassins.shp",
451 )
452
453 sb.myparam.set_param(
454 "outlets_database_fields",
455 {
456 "coord_x": "X_L93",
457 "coord_y": "Y_L93",
458 "area": "SURF",
459 "id": "CODE_SITE",
460 },
461 )
462
463 sb.newmodel("rex")
464
465 sb.rex.mysetup.update_setup(
466 {
467 "pet_directory": "/home/maxime/DATA/ETP-SFR-FRA-INTERA_L93",
468 "prcp_directory": "/home/maxime/DATA/PLUIE",
469 "qobs_directory": "/home/maxime/DATA/QOBS_SITE_60M",
470 "start_time": "2014-01-01 00:00",
471 "end_time": "2014-02-01 00:00",
472 }
473 )
474
475 sb.rex.generate_mesh(query="(SURF>20)")
476 sb.rex.myplot.plot_mesh()
477
478 graffas_prcp = np.zeros(shape=(25, 49, 365 * 24 * 5)) + 1
479 rng = np.random.default_rng()
480 graffas_prcp = (
481 graffas_prcp
482 + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
483 )
484
485 sb.rex.atmos_data_connector(input_prcp=graffas_prcp[:, :, :])
486 sb.rex.forward_run(return_options={"q_domain": True})
487
488 sb.rex.mystats.fmaxima_stats(t_axis=2.0, chunk_size=365, cumulated_maxima=True)
489
490 graffas_prcp = np.zeros(shape=(25, 49, 365 * 24 * 5)) + 1
491 rng = np.random.default_rng()
492 graffas_prcp = (
493 graffas_prcp
494 + rng.multinomial(20, [0.2, 0.8], size=graffas_prcp.shape)[:, :, :, 0]
495 )
496
497 sb.rex.atmos_data_connector(input_prcp=graffas_prcp[:, :, :])
498 sb.rex.forward_run(invert_states=True, return_options={"q_domain": True})
499
500 sb.rex.mystats.fmaxima_stats(chunk_size=365, cumulated_maxima=True)
501
502 sb.rex.mystats.fquantile_stats(
503 chunk_size=365,
504 estimate_method="MLE",
505 ncpu=6,
506 fit="gumbel",
507 compute_uncertainties=True,
508 )
509
510 sb.rex.mystats.fquantile_stats2(
511 chunk_size=365, estimate_method="MLE", ncpu=6, fit="gumbel"
512 )
513
514 sb.rex.myplot.plot_outlets_quantile()
515 sb.rex.myplot.multiplot_spatial_quantile()
516
517 array = stats.time_resample_array(
518 array=sb.rex.extra_smash_results.q_domain,
519 quantile_duration=72,
520 model_time_step=3600,
521 t_axis=2,
522 )
523 array.shape
524
525 maxima = stats.compute_maxima(
526 array=array,
527 t_axis=2,
528 nb_minimum_chunks=4,
529 chunk_size=365,
530 quantile_duration=72,
531 )
532
533 sb.rex.mystats.fquantile_stats(
534 chunk_size=365, estimate_method="MLE", ncpu=6, fit="gumbel"
535 )