smashbox 1.2
A python library developed by Hydris-hydrologie (https://www.hydris-hydrologie.fr/) and INRAE to provide an technical environment for performing hydrological simulation with Smash. This library is preloaded with data (France only), comes with preconfigured setup, own some extra features such as the in-memory atmospheric data connector, has statistical and plot capabilities. Plot and stats functions can be used as a separate module with your own Smash model object.
Keywords: hydro-smash, plot, smash, statistics

# Presentation

SmashBox is a python library developed by Inrae and Hydris-hydrologie (https://recover.paca.hub.inrae.fr/qui-sommes-nous/nos-equipes/rhax and https://www.hydris-hydrologie.fr/) to provide a Python environment for performing hydrological simulation with Smash. This library is preloaded with data (France only), comes with preconfigured setup, own some extra features such as the in-memory atmospheric data connector, has statistical and plot capabilities. Most of the features present in this module can be used with a separate Smash model object.

This module can also be used from R with the Reticulate package.

# Installation

Warning: SmashBox use the upcoming version of SMASH 1.2. Currently, to use SMashBox, you need to install SMASH manually in your Python environment from the "main" github branch (the most up-to-date version of SMASH).

To install SmashBox, you should use a python virtual environment.

```shell
python3 -m venv smashbox_env
source smashbox_env/bin/activate
pip install smashbox
```
Ensure you Python environment is loaded and load SmashBox in Python as follow:

```python
import smashbox
```

with R language you must load smashbox with the Reticulate package as follow:

```R
library(reticulate)
use_virtualenv("smashbox_env")
smashbox <- import("smashbox")
```
# API documentation  

The API documentation is hosted at https://maximejay.codeberg.page/smashbox.html. This documentation is auto-generated using pdoc (https://pdoc.dev/docs/pdoc.html).  

```bash
pdoc smashbox/ --math -o ./html  
```

To generate the doc for several modules, use for instance:

```bash
pdoc ./pyhdf5_handler/pyhdf5_handler/ smashbox/smashbox/ --math -o ./html 
```

To publish the doc on codeberg pages, you must flatten the directory. All statics pages must be at the root directory. A script can be used for this:
```bash
#modify the path to the documentation folder inside the script !
bash flatten_doc.sh
```
The script flatten_doc.sh can be found at https://codeberg.org/maximejay/pages.

# Getting started

Python and R tutorials:
----------

Before anything, few Python and R tutorials are available in the documentation section Submodules/tutorials.

Initialisation:
---------------

Here we create an smashbox object `es`. This object `es` is a container for performing hydrological simulation with smash.

```python
import smashbox
es = smashbox.SmashBox()
```

`es` has 3 methods `smashbox.init.smashbox.SmashBox.help`, `smashbox.init.smashbox.SmashBox.newmodel`, `smashbox.init.smashbox.SmashBox.delmodel` and one attribute `myparam` which refer to the class `smashbox.init.param.param`.

The smashbox parameters:
-------------------------

The class `smashbox.init.param.param` own the main top-level parameters of SmashBox. These parameters controls the data used to build an hydrological model and behaviour of SmashBox. These parameters can be managed from the attribute `myparam` of the container `es`. The attribute `myparam` refer to the class `eaysysmash.init.param.param` which contains different methods to manipulate the attribute of the sub-class `eaysysmash.init.param.smashboxparam`.

```python
sb.myparam.list_param()
asset_dir=/home/maxime/.smashbox/asset
outlets_database=/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_stations_example.csv
setup_file=/home/maxime/.smashbox/asset/setup/setup_rhax_gr4_dt3600.yaml
flowdir=/home/maxime/.smashbox/asset/flwdir/flowdir_fr_1000m.tif
bbox=None
epsg=2154
outletsID=[]
outlets_shapefile=None
smash_parameters=/home/maxime/DEV/smashbox/smashbox/asset/params
smash_parameters_dt=None
outlets_database_fields={'coord_x': 'X_L93', 'coord_y': 'Y_L93', 'area': 'SURF', 'id': 'CODE'}
enhanced_smash_input_data=False
```

The different attributes listed above are also accessible from the attribute `sb.myparam.param` which own attributes of the class `smashbox.init.param.smashboxparam`. Every parameters of this class are documented here:`smashbox.init.param.smashboxparam`.

The assets directory `asset_dir` contains the data provided by SmashBox. These data are copied in your personal directory when you load the module for the first time. SmashBox comes with few data: four outlets database for France, the flow directions at a resolution of 1km²  and few preconfigured Smash model setup. 
Notice that data in your personal asset directory is the default data used by SmashBox. So, you can edit these data and copy other personal data in this directory to quickly retrieve it each time you use SmashBox.

Every parameters can be modified using the setters or by changing directly the value of the attribute. Both following commands are equivalent:

```python
sb.myparam.set_param(
       "outlets_database",
       "/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_stations_example.csv",
    )

sb.myparam.param.outlets_database = "/home/maxime/DEV/smashbox/smashbox/asset/outlets/db_stations_example.csv"
```

Building an new model environment for performing hydrological simulations:
--------------------------------------------------------------------------

A new SmashBox model environnement can be created using the function `smashbox.init.smashbox.SmashBox.newmodel`. In this example we call our environment "graffas_zone":

```python
sb.newmodel("graffas_zone")
```

This environment is created at the top level of the object `es` and is accessible via the attribute `graffas_zone`. This environment contains functions and class-attribute to setup, build and run a Smash model, and then make plots and compute statistics. Notice that the top-level parameters `sb.myparam` are fixed and copied into the environment under the attribute `sb.graffas_zone._myparam` (hidden) and should not be modified.


Generate the Smash mesh:
------------------------

The Smash mesh is generated using the method `smashbox.mesh.mesh.generate_mesh`. But one can use the wrapper `smashbox.model.model.generate_mesh` to build the mesh directly from the model container `graffas_zone`:

```python
sb.graffas_zone.generate_mesh()
```
 The function `generate_mesh` has the option `query` to filter the outlet in the database before including them inside the mesh. See the documentation `smashbox.model.model.model.generate_mesh`.

```python
sb.graffas_zone.generate_mesh(query="('SURF' > 20 and 'SURF' < 200)", lacuna_threshold=20.) #example with filtering surface outlets and filtering with percentile of lacuna on the simulated window given by the setup (start_time and end_time).
```

The Smash setup:
----------------

The Smash model setup is stored in `mysetup.setup` which is an attribute of the class `smashbox.model.setup.setup`. From this class you can change, load, write, set, get, list all properties of the Smash setup. For instance, one can list available setup, select the one we prefer and adjust some values of the setup using the function `update_setup`:

```pyhton
sb.graffas_zone.mysetup.list_available_setup()
├── setup_rhax_gr4_dt3600.yaml
├── setup_rhax_gr4_dt900.yaml
├── setup_rhax_gr5_dt3600.yaml
└── setup_rhax_gr5_dt900.yaml

sb.graffas_zone.mysetup.load_setup("setup_rhax_gr4_dt3600.yaml")

sb.graffas_zone.mysetup.update_setup(
       {
           "pet_directory": "/nas/RHAX/DONNEES/ETP/GRILLES/ETP-SFR-FRA-INTERA_L93/",
           "prcp_directory": "/nas/RHAX/DONNEES/PLUIE/SPATIAL/COMEP/L93_1km/60M/",
           "qobs_directory": "/nas/RHAX/DONNEES/DEBITS/Extraction_sites_SCHAPI_20250226/QTVAR_to_QM/QM_60M",
       }
   )
```

Simulation with Smash :
-----------------------

The creation of the Smash model object is analogous with Smash. When the Smash setup and the mesh are ready, one just need to call the function `smashbox.model.model.model.model`:

```python
sb.graffas_zone.model()
```

Notice: if `sb.myparam.param.enhanced_smash_input_data=True`, enhanced mode is active. That means new beta functions to generate the model are used. These functions essentially extent the capabilities to handle atmospherical data (continuous pet, inter-annual pet timezone, date pattern specification, ...). Refer to the documentation of `smashbox.init.param.param`and at section `smashbox.html#extended-feature-for-reading-atmospheric-data`.

The Smash model object is created and stored behind the attribute `sb.graffas_zone.mysmashmodel`. 

Before going further, it is recommended to warm the Smash model. The function `smashbox.model.model.model.model_warmup` has been made for this. In the following example, a warmup is performed during 365 days before the start-time defined in the Smash model setup. These model states will be used for the forward run.

```python
sb.graffas_zone.model_warmup(warmup=365)
```

Then it is easy to invoke a forward run of the model:

```python
sb.graffas_zone.forward_run()
```

Notice that if the model does not exist, `forward_run` will automatically create it. Moreover, five options can be passed to the function`smashbox.model.model.model.forward_run`:
* invert_state: bool, default is False. If True, the final model states is used to start the simulation.
* warmup: int | None, default is None. Number of days of model warmup before the `start_time`.
* cost_options: Dictionary containing computation cost options for simulated and observed responses (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html)
* common_options: Dictionary containing common options with two elements, ncpu and verbose (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html)
* return_options: A dictionary with return options as explained in the Smash documentation (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html). In that case, the returned results are stored behind the attribute  `sb.graffas_zone.extra_smash_results`.


# In memory atmospheric data connector

Instead of reading atmospheric precipitation and evapotranspiration on the disk (Smash handle geotiff file), SmashBox provide a simple way to transfer a matrix with three dimension (X direction, Y direction and time) to the Smash model directly.

Let's create an numpy.array and generates periodic increasing and arbitrary rainfall pattern during `(21*18)` time-steps:

```python
nrow = sb.graffas_zone.mymesh.mesh["nrow"]
ncol = sb.graffas_zone.mymesh.mesh["ncol"]
chunk_size = 18
nb_chunk = 21
graffas_prcp = np.zeros(shape=(nrow, ncol, nb_chunk * chunk_size))
fctr = 100
for i in range(nb_chunk):
    graffas_prcp[:, :, i * chunk_size : (i + 1) * chunk_size] = (
        fctr
        * np.arange(i * chunk_size, (i + 1) * chunk_size)
        / (chunk_size * nb_chunk)
        * np.cos(np.arange(0, chunk_size * 10, 10) * np.pi / 180)
    )
graffas_prcp = np.where(graffas_prcp < 0, 0.0, graffas_prcp)
```

Then one can connect this matrix to Smash using the atmos_data_connector `smashbox.model.model.model.atmos_data_connector`:

```python
sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp, input_dt=3600.)
```
If the domain and the resolution of the input matrix differ from the domain and the resolution of the Smash model then the input matrix is cropped and resampled. Three method can be used:
* "rasterio1" : appear to be the slowest method ,`smashbox.model.atmos_data_connector.read_input_data`
* "rasterio2" : more efficient, `smashbox.model.atmos_data_connector.read_input_data2`
* "home_made_with_scipy_zoom" : default method, efficient, with resampling before cropping (minimal shifting in x and y direction), `smashbox.model.atmos_data_connector.read_input_data3`.

The bounding box and the epsg code of the domain of the input matrix should be also passed to the function to ensure best results. If the bounding box is not provided, the bounding box of the Smash mesh is used to position the input matrix.
If the input matrix has exactly the same 2D shape of the Smash mesh, the resampling and the crop operation are skipped for faster processing.

If the time-step of the precipitation differ from the model time-step, it is recommended to provide the time-step of the input precipitation. The matrix will downscaled or upscaled automatically. Use the keyword `input_dt=3600` for the atmos_data_connector.

Then we can call the forward run as follow:

```python
sb.graffas_zone.forward_run(return_options={"q_domain": True})
```

Also, if several consecutive simulation are planned with different but successive atmospheric precipitation, one ca invoke forward_run as follow:

```python
chunk_size=10
nb_chunk=int(graffas_prcp.shape[2]/chunk_size)
qdomain={}
for i in range(nb_chunck):
    sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp[:,:,i*chunk_size:(i+1)*chunk_size])
    sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
    #do some staff here...
```
In that case, the final states is used as initial states (no need to warm the model at each simulation) and the discharges for each pixel of the domain is stored in `sb.graffas_zone.extra_smash_results`. Analysis such as quantile computation can be performed for each chunk.

For deeper knowledge on this subject follow the R tutorial `smashbox.tutorials.R_Graffas_tutorial`.

# Load, calibrate and export Smash parameters:

## Default parameters

Default parameters for the asset directory are loaded. These parameters comes from the French regionalisation 2025 (J. Demargne). These parameters has been calibrated for a times-step of 900 seconds. Thus we must use a Smash model with time-step 900s. Otherwise, the parameters won't be valid.

## Converting the parameters for different simulation time-step

According to the phd of  A. Ficchi 2017, GR4 prameters can be théorically converted for different time-step. SmashBox provide a function to convert the parameter to the model time-step. Assume our model run at an  hourly time-step (we use dt=3600 in the Smash setup), the default loaded parameters will be automatically adjusted using this function. We recommend to do that just after the model creation.

```python
sb.graffas_zone.transform_parameters(original_dt=900)
```
## Optimizing the parameters

If the above solution is not satisfying, it is still possible to calibrate the Smash parameters using Smash. Moreover easySmash provide a wrapped function of smash.optimze in order to quickly optimize the parameters in a chosen period and use it for an other period. Just run:

```python
sb.graffas_zone.optimize_parameters(start_time="2010-01-01 00:00",end_time="2014-01-01 00:00", cost_options={end_warmup:"2011-01-01 00:00"})
```
Any arguments of smash.optimize can be passed to this function.

## Export and import parameters

SmashBox and Smash provide functions manually to export and import parameters:

```python
sb.graffas_zone.export_parameters(path="path/to/my/param")
sb.graffas_zone.import_parameters(path="path/to/my/param")

```

# Extended feature for reading atmospheric data:

SmashBox has an enhanced version of the smash.model() method. These functions essentially extent the capabilities to handle atmospherics data (continuous pet, inter-annual pet timezone, date pattern specification, ...) and bring more options and flexibility. Here the list of the new features:
* read atmospheric data from different source which can substitute to the first one if data are unavailable.
* Merge all reading function in one.
* Configure the pattern of the date to search in the filename. Use common date formats. Handle the occurrence number (starting from 0) with %n at the end of the date pattern. ex: %Y%m%d%H%1.
* Improve logs: new log clearly warn user about which files have been read and which files are missing.
* Fix a mistake during the reading of the evapotranspration. In SMASH the evapotranspiration from the previous day is read instead of the current one.
* Read several data source by priority: each kind of data may have different source. If one is missing, the model will read the second one, etc...
* Partially handle the reading of the continuous evapotranspiration in an operational context. To do that, a sim-link of the etp data (delayed by 1 day) is created named with the date of the current day.
* Handle time zone to shift the desegregation curve of the PET during the day.
* Improve the speed of the reading. Technically, an index of the dates and the corresponding data files is created. To efficiency run through the long list of data files and performs a regex search to match a date, a simple search algorithm is build on top of  the main loop to avoid performing a regex on thousand files. This improvement is noticeable for model running on few time-step. It is particularly important when running smash in an operational context.

The setup file must be designed as follow:

```yaml
adjust_interception: false
compute_mean_atmos: true
continuous_pet: #define if the pet data are contiuous or not
  1: true
  2: false
daily_interannual_pet: true #always use that if some pet data are not continuous
dt: 900
hydrological_module: grc
pet_date_pattern: '%Y%m%d' # the date pattern to search in the pet filename (for continuous pet)
pet_directories: # sevreal directories can be defined. The identifier correspond to the priority of the data source:  1 = higher priority
    1 : "/home/maxime/DATA/REUNION/ETPJ_continue"
    2 : "/home/maxime/DATA/REUNION/ETPJ_interannuelle_250_lnZ"
prcp_access: '%Y/%m/%d'
prcp_conversion_factor: 1.0
prcp_date_pattern: '%Y%m%d%H%M' # the date pattern to search in the prcp filename. an identifier %N can be added to target the right date on the filename, ex '%Y%m%d%H%M%2'
prcp_directories: # sevreral directories can be defined. The identifier correspond to the priority of the data source:  1 = higher priority
    1 : '/home/maxime/DATA/REUNION/PLUIE_REUNION/tests_tr/ANTILOPE_J1'
    2 : '/home/maxime/DATA/REUNION/PLUIE_REUNION/tests_tr/ANTILOPE_TR'
    3 : '/home/maxime/DATA/REUNION/PLUIE_REUNION/tests_tr/GRID_LESOL_500m'
qobs_directory: ''
read_pet: true
read_prcp: true
read_qobs: false
routing_module: lr
timezone: "Indian/Reunion" #define the timezone to adjust the hourly pet computation.
```

In order to activate this feature, you need to properly initialise the smashbox top-level parameters class `smashbox.init.param.smashboxparam` before to create the container for your simulation:

```python
sb.myparam.param.enhanced_smash_input_data=True
```
Then everything will work as explained above. Notice that the enhanced read input data features are fully compatible with the official Smash setup.

# Statistical capabilities
 
The statistical features are handled by the submodule `smashbox.stats`. `smashbox.stats.mystats.mystats` features some methods ready to be applied on any smashbox object. `smashbox.stats.stats` owns generic functions that can be used outside the SmashBox object.

Misfit criterion at the outlets
------------------------------

If discharges observations are available, misfit criterion can be computed at each gauge. The available misfit criterion are:
* nse : Nash Sutcliff criteria $nse= 1 - \frac{\sum_{i=1}^{i=N} (Q^i_{sim}-Q^i_{obs})^2 }{ \sum_{i=1}^{i=N} (Q^i_{obs}-\overline{Q_{obs}})^2 }$
* nnse : Normalized Nash Sutcliffe criteria (Nossent and Bauwens, 2012) $nnse = \frac{1}{2-NSE} \ \in[0:1]$
* mse : The Mean Squared Error, $mse=\frac{1}{n}\sum^n_{k=1}(f_k - o_k)^2$
* rmse : Root Mean Square Error, $rmse=\sqrt{\frac{1}{n}\sum^n_{k=1}(f_k - o_k)^2}$
* nrmse : Normalized Root Mean Square Error, $nrmse=\frac{rmse}{\overline{Q_{obs}}}$
* se: Square Error, $se=\sum{(Q_{obs}-Q_{sim})^2}$
* mae: The Mean Absolute Error: $mae=\frac{1}{N}\sum_{i=0}^{N}{abs(Qsim^i-Qobs^i)}$
* mape: The Mean Absolute Percentage Error : $mape=\frac{1}{N}\sum_{i=0}{N}{abs(\frac{Qsim^i-Qobs^i}{Qobs^i})}$
* lgrm: The Logarithmic Error : $lgrm=\sum$_{i=0}{N}Qobs^i \times ln(\frac{Qsim^i}{Qobs^i})^2$
* kge: The Kling Gupta Efficiency criteria.

To compute all these misfits, just run `sb.graffas_zone.mystats.fmisfit_stats()` function, and all results will be stored in `sb.graffas_zone.mystats.misfit_stats`:

```python
sb.graffas_zone.mystats.fmisfit_stats()
sb.graffas_zone.mystats.misfit_stats.results.nse
```

Results are stored in attribute `sb.graffas_zone.mystats.misfit_stats.results`. Criterion can be computed separately using the misfit functions in `sb.graffas_zone.mystats.misfit_stats`: 2 separates functions can be used for each statistic can be used. These functions are all named by the statistics name. But some have the prefix `sm_` which stand for `smash_`. For those functions the statistics are computed using the smash library. We advice to use function prefixed with `sm_` because the Smash library is more robust and has been fully tested.

Outlet discharges statistics
----------------------------

Basic statistics can be computed on the observed and simulated discharge at every outlets. These statistic are the mean, the median, the 20% and 80% percentile of the discharge, the maximum and the minimum. The function `sb.graffas_zone.mystats.foutlets_stats()` compute all of these statistics and the results are stored in `sb.graffas_zone.mystats.outlets_stats`:

```python
sb.graffas_zone.mystats.foutlets_stats()
sb.graffas_zone.mystats.outlets_stats.results_sim.maximum
sb.graffas_zone.mystats.outlets_stats.results_obs.maximum
```

Spatial discharges statistics
-----------------------------

Basic spatial statistics can be computed on the simulated discharge at every active pixel of the domain. These statistic are the mean, the median, the 20% and 80% percentile of the discharge, the maximum and the minimum. This require to ran the forward model with return option `{q_domain:True}`. The function `sb.graffas_zone.mystats.fspatial_stats()` compute all of these statistics and the results are stored in `sb.graffas_zone.mystats.spatial_stats`:

```python
sb.graffas_zone.mystats.fspatial_stats()
sb.graffas_zone.mystats.spatial_stats.results.maximum
```

Spatial discharges quantile computation
---------------------------------------

SmashBox has the capabilities to compute the discharges quantile spatially for various return period. The following process is applied:

* Perform a forward run with return option `{q_domain:True}` so that the spatial discharges is saved for all time-steps.
* Resample the spatial discharge over time if needed (depending the duration of the quantile)
* Compute the annual maxima (annual chunck size) for each active pixel of the domain (notice that the chunck size can be adjusted but default is 365 days)
* Adjust an extrem law, such as Gumbel or GEV, on the maxima values. This part use scipy algorithm.
* Compute the quantile for the given return period

To perform this computation one use `smashbox.stats.mystats.fstats_quantile` function. In the following example, the length of the chronicle is about 400 hours (16 days) and the chunk_size is set to 4 days (instead of 365 days) to have at least four chunks to compute the quantile (in other word, the return period 'unit' is 4 days). By default the estimated return period T are `[2, 5, 10, 20, 50, 100]` (unit depend on the chunk size), and the quantile duration are `[1, 2, 3, 4, 6, 12, 24, 48, 72]` hours.

```python
sb.graffas_zone.mystats.fquantile_stats(
       chunk_size=4, estimate_method="MLE", ncpu=6, fit="gumbel"
   )
```

The results of the quantiles are then stored for each duration behind the attribute `sb.graffas_zone.mystats.quantile_stats`:

```python
sb.graffas_zone.mystats.quantile_stats.__dict__.keys()
Out[9]: dict_keys(['spatial_cumulated_maxima', 'spatial_cumulated_maxima_outlets', 'spatial_maxima', 'spatial_maxima_outlets', 'spatial_quantile', 'spatial_quantile_outlets', 'Quantile_1h', 'Quantile_2h', 'Quantile_3h', 'Quantile_4h', 'Quantile_6h', 'Quantile_12h', 'Quantile_24h', 'Quantile_48h', 'Quantile_72h'])
```

The attributes `spatial_quantile` and `spatial_maxima` store respectively the spatial quantile and the spatial maxima in an 4D arrays (respectively of shape `(nbx, nby, duration, T)` and `(nbx, nby, nb_chunk, duration)`) 

The attributes `spatial_quantile_outlets` and `spatial_maxima_outlets` store respectively the quantile and the maxima at every outlets in an 3D arrays (respectively of shape `(nb_outlets, duration, T)` and `(nb_outlets, nb_chunk, duration)`)

The attributes `spatial_cumulated_maxima` `spatial_cumulated_maxima_outlets` and  store the maximum discharges (maxima) respectively spatially and for outlets only, but it accumulated these values after each call of `sb.graffas_zone.mystats.fmaxima_stats`. Shapes are `(nbx, nby, nb_chunk*nb_simulation, duration)` and `(nb_outlets, nb_chunk*nb_simulation, duration)`. `nb_simulation` represent the number of smash run and the number of call of `sb.graffas_zone.mystats.fmaxima_stats`.

Remark: despite parallel computing is used (on the pixel), this computation may consume an high quantity of memory and can take a long time. Depending the size of the grid and for how much duration the quantile are computed, the function `sb.graffas_zone.mystats.fstats_quantile2` perform the computation in parallel for each duration. The second method could be more efficient in some case.

For extremely long rainfall chronic the forward run can be split to save memory. If our goal is to compute the discharges quantile, it is possible to only compute the maxima after each forward run. These maxima will be automatically stacked (with option `cumulated_maxima=True`) and the Gumbel or GEV adjustment can be achieved for each duration after all:

```python
from smashbox import stats
chunk_size=4
nb_chunk=int(graffas_prcp.shape[2]/chunk_size)
qdomain={}
for i in range(nb_chunck):
    sb.graffas_zone.atmos_data_connector(input_prcp=graffas_prcp[:,:,i*chunk_size:(i+1)*chunk_size])
    sb.graffas_zone.forward_run(invert_states=True, return_options={"q_domain": True})
    sb.graffas_zone.mystats.fmaxima_stats(chunk_size=4, cumulated_maxima=True)

results={}
for i, duration in enumerate([1, 2, 3, 4, 6, 12, 24, 48, 72]):
    results["duration"] = stats.fit_quantile(
        maxima=sb.graffas_zone.mystats.quantile_stats.spatial_cumulated_maxima[:, :, :, i],
        t_axis=2,
        return_periods=[2, 5, 10, 20, 50, 100],
        fit="gumbel",
        estimate_method="MLE",
        quantile_duration=duration,
        ncpu=6,
    )
```

# Multi-model statistics (Beta)

Several model containers can be created at the root of the SmashBox object. Each model container owns statistic (quantile, misfit, spatial criteria, ...). SmashBox has a feature to aggregate all model statistics and compute the mean, the median, the variance the min, the max and the distribution. These statistics will be computed for each model where the bounding-box match the one in `sb.myparam` attribute. The class `smashbox.init.multimodel_statistics.compute_multimodel_statistics` is instanced when a model container is created.

```python
sb.multimodel_statistics.compute_multimodel_statistics().
```


# Export results

We can save the smashbox container to an HDF5 database:

```python
sb.graffas_zone.save_model_container(path_to_hdf5="my_smashbox_data.hdf5")
```

We can also export the data of the smashbox container to a python dictionary if `path_to_hdf5è is left empty:

```python
dict_sb=sb.graffas_zone.save_model_container()
```
This dictionary will be eventually converted to a R list object and can be saved in a Rdata.


# Plotting capabilities

SmashBox has some plotting capabilities. The plot are produced with the Matplotlib library. We use the mpld3 librairy to render the plot in html and visualize it directly in your favorite navigator. Every plot can be configured and saved on the disk as any format that matplotlib support. Plotting features are handled by the submodule `smashbox.plot`.

Graphics settings
-----------------

Every plotting functions of SmashBox have options to configure the plot. These options can be passed to any function as a dictionary:

* `fig_settings`: figure settings are related to the class `smashbox.plot.plot.fig_properties`. These options configure the figure output. With a dictionary, one can define and change any attribute of the class `smashbox.plot.plot.fig_properties`.
* `ax_settings`: Ax settings refers to the class `smashbox.plot.plot.ax_properties`. These options are useful to change the appaerance of the graphics (changing the matplotlib ax properties). With a dictionary, one can define and change any attribute of the class `smashbox.plot.plot.ax_properties`.
* `plot_settings`: plot settings refers to the class `smashbox.plot.plot.plot_properties`. Some function accept `plot_settings` options. `plot_settings` serves to configure the line of the matplotlib plot function (colour, line style, line width).

If no output is configured, i.e `fig_settings.figname=None`, the figure will popup to the screen using the matplotlib UI. However, from R language, reticulate is not able to display the figure using matplotlib UI. As workaround, SmashBox use mpld3 to output the figure in html and display it directly in your navigator.

Plotting the mesh
-----------------

With the function `smashbox.plot.myplot.myplot.plot_mesh`, SmashBox is able to plot the mesh of the Smash model in order to visualize the hydrographic network, the catchment and the localisation of the different outlets. In the following example, the figure will be saved in the file mesh.png, its size will measure 20*20 inches and the global font will be multiply by 2.

```python
sb.graffas_zone.myplot.plot_mesh(
       fig_settings={"figname": "mesh.png", "xsize": 20, "ysize": 20, "font_ratio": 2},
       ax_settings={"title_fontsize": 32},
   )
```

![image info](./images/mesh.png)

Two other function helps to characterize the catchment surface error between the mesh and the real surface. These functions are `smashbox.plot.myplot.myplot.plot_catchment_surface_consistency` and `smashbox.plot.myplot.myplot.plot_catchment_surface_error`.

```python
sb.graffas_zone.myplot.plot_catchment_surface_consistency(
       fig_settings={"figname": "mesh_surface_consistency.png", "xsize": 20, "ysize": 20, "font_ratio": 2},
       ax_settings={"title_fontsize": 32},
   )
sb.graffas_zone.myplot.plot_catchment_surface_error(
       fig_settings={"figname": "mesh_surface_error.png", "xsize": 20, "ysize": 20, "font_ratio": 2},
       ax_settings={"title_fontsize": 32},
   )
```
![image info](./images/mesh_surface_consistency.png)
![image info](./images/mesh_surface_error.png)

Plotting the parameters
-----------------------

It is easy to plot the model parameters. Just use both functions `smashbox.plot.myplot.myplot.plot_parameters` and `smashbox.plot.myplot.myplot.multiplot_parameters`.

```python
sb.graffas_zone.myplot.plot_parameters(
    param="cp",
    fig_settings={
        "figname": "param_cp.png",
        "xsize": 10,
        "ysize": 10,
    },
    ax_settings={"title_fontsize": 12},
)
sb.graffas_zone.myplot.multiplot_parameters(
    fig_settings={
        "figname": "multiplot_param.png",
        "xsize": 20,
        "ysize": 20,
        "font_ratio": 1,
    },
    ax_settings={"title_fontsize": 10, "font_ratio": 2},
)
```

![image info](./images/param_cp.png)
![image info](./images/multiplot_param.png)

Plotting the quantile
---------------------

SmashBox allow to plot the discharges quantile in many different ways. Four plot are possible.

The function `smashbox.plot.myplot.myplot.plot_xy_quantile` allow you to plot the quantile adjustment at any coordinate of the domain:

```python
sb.graffas_zone.myplot.plot_xy_quantile(
       X=100, Y=10, fig_settings={"font_ratio": 2}, ax_settings={"legend_fontsize": 20}
    )
```

![image info](./images/quantile_XY.png)

It is also possible to plot the quantile at every outlets of the domain using the function `smashbox.plot.myplot.myplot.plot_outlets_quantile`. In the following example the figure will display in html in your navigator:

```python
sb.graffas_zone.myplot.plot_outlets_quantile(
    html_show=True,
    fig_settings={
        "figname": "quantile_outlets.png",
        "xsize": 10,
        "ysize": 10,
    },
)
``` 

![image info](./images/quantile_outlets.png)

The function `smashbox.plot.myplot.myplot.plot_spatial_quantile` allows to plot the spatial discharges quantile on a map for a given return period and duration.

```python
sb.graffas_zone.myplot.plot_spatial_quantile(
    T=2,
    duration=1,
    fig_settings={
        "figname": "spatial_quantile.png",
        "xsize": 10,
        "ysize": 10,
    },
    ax_settings={"font_ratio": 1.5},
)
```

![image info](./images/spatial_quantile.png)

The latest function `smashbox.plot.myplot.myplot.multiplot_spatial_quantile` allow to multiplot the map of the quantile for a given duration and for all return period.

```python
sb.graffas_zone.myplot.multiplot_spatial_quantile(
    duration=1,
    fig_settings={
        "figname": "../images/multiplot_spatial_quantile.png",
        "xsize": 10,
        "ysize": 10,
    },
    ax_settings={"font_ratio": 0.8},
)
```

![image info](./images/multiplot_spatial_quantile.png)

Plotting the statistics
-----------------------

Any spatial statistic can be plotted on a map using the function `smashbox.plot.myplot.myplot.plot_spatial_stats` as follow:

```python
sb.graffas_zone.myplot.plot_spatial_stats(
    fig_settings={
        "figname": "../images/spatial_stats_max.png",
        "xsize": 10,
        "ysize": 10,
    },
    ax_settings={"font_ratio": 1.0},
)
```

![image info](./images/spatial_stats_max.png)


Plotting the misfit criterion
-----------------------------

The misift criterion can be plotted using the function `smashbox.plot.myplot.myplot.plot_misfit`. It is a bar plot of the misfit scores at every outlets.

```python
sb.graffas_zone.myplot.plot_misfit(
    misfit="nse",
    fig_settings={
        "figname": "../images/nse.png",
        "xsize": 10,
        "ysize": 7,
    },
)
```
![image info](./images/nse.png)

Multiplot of every misfit criteria can be plotted directly using the function `smashbox.plot.myplot.myplot.plot_misfit`:

```python
sb.graffas_zone.myplot.multiplot_misfit(
    fig_settings={
        "figname": "../images/multiplot_misfit.png",
        "xsize": 10,
        "ysize": 10,    
    },
)
```

![image info](./images/multiplot_misfit.png)

Finally a map of the misfit criteria at the outlets is plotted using the function `smashbox.plot.myplot.myplot.plot_misfit_map`:

```python
fig, ax = sb.graffas_zone.myplot.plot_misfit_map(
    misfit="nnse",
    fig_settings={"figname": "../images/misfit_map.png", "xsize": 10, "ysize": 10},
    ax_settings={"font_ratio": 1},
)
```

![image info](./images/misfit_map.png)

Plotting hydrogram
------------------

Hydrogram of the discharges can be plotted for different outlets with the function `smashbox.plot.myplot.myplot.plot_hydrograph`:

```python
sb.graffas_zone.myplot.plot_hydrograph(
    fig_settings={
        "figname": "../images/hydrogram.png",
        "xsize": 10,
        "ysize": 10,
    },
)
```

![image info](./images/hydrogram.png)

# SmashBox tree structure object:

SmashBox comes with many other methods than described above. One can display the tree structure of the SmashBox object to better understand how it is structured and help to discover the existing methods associated to each object.

```python
smashbox.tools.tools.build_object_tree(es,"es")

Object sb (SmashBox)
├── Method: delmodel()
├── Object graffas_zone (model)
│   ├── Method: atmos_data_connector()
│   ├── Method: export_parameters()
│   ├── Unknown attribute type: extra_smash_results = <class 
│   │   'smash.core.simulation.run.run.ForwardRun'>
│   ├── Method: forward_run()
│   ├── Method: generate_mesh()
│   ├── Method: import_parameters()
│   ├── Method: model()
│   ├── Method: model_warmup()
│   ├── Object mymesh (mesh)
│   │   ├── Attribute: catchment_polygon = None
│   │   ├── Method: generate_mesh()
│   │   ├── Method: load_mesh()
│   │   ├── Attribut:  mesh = dict_keys(['xres', 'yres', 'xmin', 'ymax', 'epsg',
│   │   │   'nrow', 'ncol', 'dx', 'dy', 'flwdir', 'flwdst', 'flwacc', 'npar', 
│   │   │   'ncpar', 'cscpar', 'cpar_to_rowcol', 'flwpar', 'nac', 'active_cell',
│   │   │   'ng', 'gauge_pos', 'code', 'area', 'area_dln'])
│   │   └── Method: write_mesh()
│   ├── Object myplot (myplot)
│   │   ├── Method: multiplot_misfit()
│   │   ├── Method: multiplot_parameters()
│   │   ├── Method: multiplot_spatial_quantile()
│   │   ├── Method: plot_catchment_surface_consistency()
│   │   ├── Method: plot_catchment_surface_error()
│   │   ├── Method: plot_hydrograph()
│   │   ├── Method: plot_mesh()
│   │   ├── Method: plot_misfit()
│   │   ├── Method: plot_misfit_map()
│   │   ├── Method: plot_outlet_stats()
│   │   ├── Method: plot_outlets_quantile()
│   │   ├── Method: plot_parameters()
│   │   ├── Method: plot_spatial_quantile()
│   │   ├── Method: plot_spatial_stats()
│   │   └── Method: plot_xy_quantile()
│   ├── Object mysetup (setup)
│   │   ├── Method: get_setup()
│   │   ├── Method: list_available_setup()
│   │   ├── Method: load_setup()
│   │   ├── Method: set_setup()
│   │   ├── Attribut:  setup = dict_keys(['adjust_interception', 
│   │   │   'daily_interannual_pet', 'dt', 'end_time', 'snow_module', 
│   │   │   'hydrological_module', 'pet_directory', 'prcp_conversion_factor', 
│   │   │   'prcp_directory', 'qobs_directory', 'read_pet', 'read_prcp', 
│   │   │   'read_qobs', 'routing_module', 'start_time'])
│   │   ├── Attribute: setup_file = 
│   │   │   '/home/maxime/.smashbox/asset/setup/setup_rhax_gr4_dt3600.yaml'
│   │   ├── Method: update_setup()
│   │   └── Method: write_setup()
│   ├── Unknown attribute type: mysmashmodel = <class 
│   │   'smash.core.model.model.Model'>
│   ├── Object mystats (mystats)
│   │   ├── Method: fmaxima_stats()
│   │   ├── Method: fmisfit_stats()
│   │   ├── Method: foutlets_stats()
│   │   ├── Method: fquantile_stats()
│   │   ├── Method: fquantile_stats2()
│   │   ├── Method: fspatial_stats()
│   │   ├── Object misfit_stats (misfit_stats)
│   │   │   ├── Method: kge()
│   │   │   ├── Method: lgrm()
│   │   │   ├── Method: mae()
│   │   │   ├── Method: mape()
│   │   │   ├── Method: mse()
│   │   │   ├── Method: nnse()
│   │   │   ├── Method: nrmse()
│   │   │   ├── Method: nse()
│   │   │   ├── Object results (misfit_results)
│   │   │   │   ├── Attribute: kge = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: lgrm = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: mae = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: mape = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: mse = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: nnse = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: nrmse = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: nse = <class 'numpy.ndarray'>
│   │   │   │   ├── Attribute: rmse = <class 'numpy.ndarray'>
│   │   │   │   └── Attribute: se = <class 'numpy.ndarray'>
│   │   │   ├── Method: rmse()
│   │   │   ├── Method: se()
│   │   │   ├── Method: sm_kge()
│   │   │   ├── Method: sm_lgrm()
│   │   │   ├── Method: sm_mae()
│   │   │   ├── Method: sm_mape()
│   │