.. _user_guide.quickstart.real_case_france:

==================
Real case - France
==================

A real case on the entire France is considered.

First, open a Python interface:

.. code-block:: none

    python3
    
-------
Imports
-------

.. ipython:: python
    
    import smash
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.colors import LogNorm, SymLogNorm

---------------------   
Model object creation
---------------------

Creating a :class:`.Model` requires two input arguments: ``setup`` and ``mesh``. For this case, it is possible to directly load the both input dictionnaries using the :meth:`smash.load_dataset` method.

.. ipython:: python

    setup, mesh = smash.load_dataset("France")

.. note::

    The :ref:`user_guide.quickstart.real_case_france` is very similar to the :ref:`user_guide.quickstart.real_case_cance` except that here we do not consider catchments but a domain to make a simulation.

Setup argument
**************
    
Compared to the :ref:`user_guide.quickstart.real_case_cance`, less options have been filled in the ``setup`` dictionary.

.. ipython:: python

    setup

The only options which had been added to the setup is ``save_qsim_domain``. By default, to avoid storing in memory the entire grid of simulated discharge (especially on non active cells) when working with catchments this option is set to False.
But in case of simulation on a domain, we need to store the simulated discharge on the whole domain.

.. note::
    We won't read any observed discharge because we are not working on catchment entities but domain.

.. _user_guide.quickstart.real_case_france.mesh_argument:

Mesh argument
*************

Mesh composition
''''''''''''''''

.. ipython:: python

    mesh.keys()

Same as ``setup``, compared to the :ref:`user_guide.quickstart.real_case_cance`, less options have been filled in the ``mesh`` dictionary. All gauge specific options are not used in this case (i.e. ``code``, ``area``, ``gauge_pos`` ...).

We can still visualize the extent of the grid, the flow directions and flow accumulation.

.. ipython:: python

    mesh["nrow"], mesh["ncol"]

.. ipython:: python
    
    plt.imshow(mesh["flwdir"]);
    plt.colorbar(label="Flow direction (D8)");
    @savefig user_guide.quickstart.real_case_france.flwdir.png
    plt.title("Real case - France - Flow direction");


.. ipython:: python
    
    plt.imshow(mesh["flwacc"], norm=LogNorm());
    plt.colorbar(label="Flow accumulation (nb cells)");
    @savefig user_guide.quickstart.real_case_france.flwacc.png
    plt.title("Real case - France - Flow accumulation");

This mesh can also be automatically generated by providing to the :meth:`smash.generate_mesh` method the France flow directions and the bouding box ``(xmin, xmax, ymin, ymax)``.

.. ipython:: python

    flwdir = smash.load_dataset("flwdir")

    france_bbox = (100_000, 1_250_000, 6_050_000, 7_125_000)

    mesh2 = smash.generate_mesh(
        path=flwdir,
        bbox=france_bbox
    )

This ``mesh2`` created is a dictionnary which is identical to the ``mesh`` loaded with the :meth:`smash.load_dataset` method.

.. ipython:: python
    
    mesh2["nrow"], mesh2["ncol"]

Finally, create the :class:`smash.Model` object using the ``setup`` and ``mesh`` loaded.

.. ipython:: python

    model = smash.Model(setup, mesh)

    model

---
Run 
---

Forward run
***********

Make a forward run using the :meth:`.Model.run` method.

.. ipython:: python

    model.run(inplace=True);

We can visualize the simulated discharges after a forward run on the whole domain. Here for the last time step of simulation.

.. ipython:: python

    qsim = model.output.qsim_domain[..., -1]
    qsim = np.where(model.mesh.active_cell == 0, np.nan, qsim)

    plt.imshow(qsim, norm=SymLogNorm(1e-4));
    plt.colorbar(label="Discharge $(m^3/s)$");
    @savefig user_guide.quickstart.real_case_france.qsim.png
    plt.title("Real case - France - Discharge");

We can visualize the precipitation for the same time step. In addition to masking the non active cells, we mask the cells where there is no precipitation data (i.e. precipitation lower than 0).

.. ipython:: python

    prcp = model.input_data.prcp[..., -1]
    prcp = np.where(
        np.logical_or(model.mesh.active_cell == 0, prcp < 0),
        np.nan,
        prcp
    )

    plt.imshow(prcp);
    plt.colorbar(label="Precipitation (mm/h)");
    @savefig user_guide.quickstart.real_case_france.prcp.png
    plt.title("Real case - France - Precipitation");

.. ipython:: python
    :suppress:

    plt.close('all')