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