Precipitation indices#

This section aims to go into detail on how to compute and vizualize precipitation indices.

First, open a Python interface:

python3

Imports#

In [1]: import smash

In [2]: import numpy as np

In [3]: import matplotlib.pyplot as plt

Model object creation#

To compute precipitation indices, you need to create a smash.Model object. For this case, we will use the Cance dataset used in the User Guide section: Real case - Cance.

Load the setup and mesh dictionaries using the smash.load_dataset() method and create the smash.Model object.

In [4]: setup, mesh = smash.load_dataset("Cance")

In [5]: model = smash.Model(setup, mesh)

Precipitation indices computation#

Once the smash.Model is created. The precipitation indices, computed for each gauge and each time step, are available using the Model.prcp_indices() method.

In [6]: res = model.prcp_indices();

In [7]: res
Out[7]: 
  d1: array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
  d2: array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
 std: array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
  vg: array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

Note

NaN value means that there is no precipitation at this specific gauge and time step and therefore no precipitation indices.

The precipitation indices results are represented as a smash.PrcpIndicesResult object containning 4 attributes which are 4 different indices:

std : The precipitation spatial standard deviation,
d1 : The first scaled moment [Zoccatelli et al., 2011],
d2 : The second scaled moment [Zoccatelli et al., 2011],
vg : The vertical gap [Emmanuel et al., 2015] .

Each attributes (i.e. precipitation indices) of the smash.PrcpIndicesResult object is a numpy.ndarray of shape (number of gauge, number of time step)

In [8]: res.std
Out[8]: 
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

In [9]: res.std.shape
Out[9]: (3, 1440)

Precipitation indices description#

Precipitation spatial standard deviation (std)#

Simply the standard deviation.

Scaled moments (d1 and d2)#

The spatial scaled moments are described in [Zoccatelli et al., 2011] in the section 2 Spatial moments of catchment rainfall: definitions.

The first scaled moment \(\delta 1\) describes the distance of the centroid of catchment rainfall with respect to the average value of the flow distance (i.e. the catchment centroid). Values of \(\delta 1\) close to 1 reflect a rainfall distribution either concentrated close to the position of the catchment centroid or spatially homogeneous, with values less than one indicating that rainfall is distributed near the basin outlet, and values greater than one indicating that rainfall is distributed towards the catchment headwaters.

The second scaled moment \(\delta 2\) describes the dispersion of the rainfall-weighted flow distances about their mean value with respect to the dispersion of the flow distances. Values of \(\delta 2\) close to 1 reflect a uniform-like rainfall distribution, with values less than 1 indicating that rainfall is characterised by a unimodal distribution along the flow distance. Values greater than 1 are generally rare, and indicate cases of multimodal rainfall distributions.

Vertical gap (vg)#

The vertical gap is described in [Emmanuel et al., 2015] in the section 5.2 The proposed indexes.

VG values close to zero indicate a rainfall distribution over the catchment revealing weak spatial variability. The higher the VG value, the more concentrated the rainfall over a small part of the catchment.

Precipitation indices visualization#

Most of the precipitation indices are calculated based on flow distances. As a reminder and to facilitate the understanding of the indices values with respect to the catchment outlet and headwaters, the flow distances of the catchment are plotted below.

In [10]: flwdst = np.where(model.mesh.active_cell==0, np.nan, model.mesh.flwdst)

In [11]: plt.imshow(flwdst);

In [12]: plt.colorbar(label="Flow distance (m)");

In [13]: plt.title("Cance - Flow distance");

../../_images/user_guide.in_depth.prcp_indices.flwdst.png

Precipitation spatial standard deviation (std)#

First getting the indexes (i.e. the time step) where occured the minimum and maximum. We use the methods numpy.nanargmin and numpy.nanargmax to find the indexes ignoring NaN’s.

In [14]: ind_min = np.nanargmin(res.std[0,:])

In [15]: ind_max = np.nanargmax(res.std[0,:])

In [16]: ind_min, ind_max
Out[16]: (168, 666)

Then, we can visualize the precipitation grids at this time steps masking the non active cells.

In [17]: f, ax = plt.subplots(1, 2, tight_layout=True)

In [18]: ma = (model.mesh.active_cell == 0)

In [19]: prcp_min = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_min])

In [20]: prcp_max = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_max])

In [21]: map_min = ax[0].imshow(prcp_min);

In [22]: f.colorbar(map_min, ax=ax[0], fraction=0.05);

In [23]: ax[0].set_title("Minimum std");

In [24]: map_max = ax[1].imshow(prcp_max);

In [25]: f.colorbar(map_max, ax=ax[1], fraction=0.05, label="Precipitation (mm)");

In [26]: ax[1].set_title("Maximum std");

../../_images/user_guide.in_depth.prcp_indices.std.png

And the associated indices values

In [27]: std_min = res.std[0, ind_min]

In [28]: std_max = res.std[0, ind_max]

In [29]: std_min, std_max
Out[29]: (0.0051030866, 8.763975)

Scaled moments (d1 and d2)#

The same applies to scaled moments, except that we will also visualize the precipitation maps where the scaled moments are closed to 1.

In [30]: ind_min = np.nanargmin(res.d1[0,:])

In [31]: ind_max = np.nanargmax(res.d1[0,:])

In [32]: ind_one = np.nanargmin(np.abs(res.d1[0,:] - 1))

In [33]: ind_min, ind_max, ind_one
Out[33]: (529, 168, 1203)

Then, we can visualize the precipitation grids at this times step masking the non active cells.

In [34]: f, ax = plt.subplots(2, 2, tight_layout=True)

In [35]: ma = (model.mesh.active_cell == 0)

In [36]: prcp_min = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_min])

In [37]: prcp_max = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_max])

In [38]: prcp_one = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_one])

In [39]: map_min = ax[0,0].imshow(prcp_min);

In [40]: f.colorbar(map_min, ax=ax[0,0]);

In [41]: ax[0,0].set_title("Minimum d1");

In [42]: map_max = ax[0,1].imshow(prcp_max);

In [43]: f.colorbar(map_max, ax=ax[0,1]);

In [44]: ax[0,1].set_title("Maximum d1");

In [45]: map_one = ax[1,0].imshow(prcp_one);

In [46]: f.colorbar(map_one, ax=ax[1,0], label="Precipitation (mm)");

In [47]: ax[1,0].set_title("Close to one d1");

In [48]: ax[1,1].axis('off');

../../_images/user_guide.in_depth.prcp_indices.d1.png

And the associated indices values

In [49]: d1_min = res.d1[0, ind_min]

In [50]: d1_one = res.d1[0, ind_one]

In [51]: d1_max = res.d1[0, ind_max]

In [52]: d1_min, d1_one, d1_max
Out[52]: (0.02219003, 1.0007353, 1.581147)

Applying the same for d2

In [53]: ind_min = np.nanargmin(res.d2[0,:])

In [54]: ind_max = np.nanargmax(res.d2[0,:])

In [55]: ind_one = np.nanargmin(np.abs(res.d2[0,:] - 1))

In [56]: ind_min, ind_max, ind_one
Out[56]: (648, 633, 527)

Then, we can visualize the precipitation grids at this time steps masking the non active cells.

In [57]: f, ax = plt.subplots(2, 2, tight_layout=True)

In [58]: ma = (model.mesh.active_cell == 0)

In [59]: prcp_min = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_min])

In [60]: prcp_max = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_max])

In [61]: prcp_one = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_one])

In [62]: map_min = ax[0,0].imshow(prcp_min);

In [63]: f.colorbar(map_min, ax=ax[0,0]);

In [64]: ax[0,0].set_title("Minimum d2");

In [65]: map_max = ax[0,1].imshow(prcp_max);

In [66]: f.colorbar(map_max, ax=ax[0,1]);

In [67]: ax[0,1].set_title("Maximum d2");

In [68]: map_one = ax[1,0].imshow(prcp_one);

In [69]: f.colorbar(map_one, ax=ax[1,0], label="Precipitation (mm)");

In [70]: ax[1,0].set_title("Close to one d2");

In [71]: ax[1,1].axis('off');

../../_images/user_guide.in_depth.prcp_indices.d2.png

And the associated indices values

In [72]: d2_min = res.d2[0, ind_min]

In [73]: d2_one = res.d2[0, ind_one]

In [74]: d2_max = res.d2[0, ind_max]

In [75]: d2_min, d2_one, d2_max
Out[75]: (1.6027083e-07, 0.99878937, 2.9966152)

Vertical gap (vg)#

Finally, the same applies to the vertical gap.

In [76]: ind_min = np.nanargmin(res.vg[0,:])

In [77]: ind_max = np.nanargmax(res.vg[0,:])

In [78]: ind_min, ind_max
Out[78]: (105, 168)

Then, we can visualize the precipitation grids at this time steps masking the non active cells.

In [79]: f, ax = plt.subplots(1, 2, tight_layout=True)

In [80]: ma = (model.mesh.active_cell == 0)

In [81]: prcp_min = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_min])

In [82]: prcp_max = np.where(ma, np.nan, model.input_data.prcp[:,:,ind_max])

In [83]: map_min = ax[0].imshow(prcp_min);

In [84]: f.colorbar(map_min, ax=ax[0], fraction=0.05);

In [85]: ax[0].set_title("Minimum vg");

In [86]: map_max = ax[1].imshow(prcp_max);

In [87]: f.colorbar(map_max, ax=ax[1], fraction=0.05, label="Precipitation (mm)");

In [88]: ax[1].set_title("Maximum vg");

../../_images/user_guide.in_depth.prcp_indices.vg.png

And the associated indices values

In [89]: vg_min = res.vg[0, ind_min]

In [90]: vg_max = res.vg[0, ind_max]

In [91]: vg_min, vg_max
Out[91]: (0.050956417, 0.8981723)