3.1. Visualisation of stochastic data

3.1.1. Realizations

PCRaster includes the Aguila software for visualisation of spatio-temporal stochastic data. With Aguila, you can easily explore large multi-dimensional data sets. For the exercises below, you will need to look up some options in the Aguila manual, available at https://pcraster.geo.uu.nl/pcraster/latest/documentation/pcraster_aguila/index.html.

To learn how to visualize stochastic data you will use hydrological model outputs from a model of a catchment in Spain.

Aguila is started from a command shell. Open a command shell and go to the directory containing your data for the exercises. Type cd to change directories. In the directory, you will find the folder visualisation containing the data for the visualisation practical. In this directory, display the digital elevation model, the vegetation map, and the local drain direction map:

aguila dem.map veg.map ldd.map <Enter>

Use the menu items (or right click on the legend to get more options) to zoom in/out, to change the color palette, and to change the drawing mode to contours. Also try changing the number of classes shown (or number of contours).

To represent saturated conductivity as a stochastic variable, realizations were drawn from probability distributions. These were used in the stochastic hydrological model. In PCRaster Python, realizations are stored in subdirectories numbered 1, 2,..,N, with N the number of realizations. Use the command dir (or ls on linux) to check the number of realizations available in the data set.

Also, go to the folder 1 and type dir to have a look at its contents. You will see that it contains a ksat.map which is a realization of saturated conductivity (m/h).

Go back to the main folder visualisation.

You can display realizations of static data (without time component) by typing:

aguila --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' ksat

This opens the realizations in separate windows. It is more convenient to display in one window, by typing:

aguila --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4 ksat

You can combine the visualisation with the vegetation map:

aguila veg.map --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4 ksat

Question: What is the value of saturated conductivity at the 25th, 50th and 75th percentile, respectively, in the unit dense herbs? Make a rough estimate. Use the realizations at a single cell location.

  1. 25th percentile: 0.064; 50th percentile: 0.061; 75th percentile: 0.059.

  2. 25th percentile: 0.059; 50th percentile: 0.061; 75th percentile: 0.064.

  3. 25th percentile: 0.059; 50th percentile: 0.059; 75th percentile: 0.059.

  4. 25th percentile: 0.061; 50th percentile: 0.059; 75th percentile: 0.064.

Correct answers: b.

Feedback: None

Now, lets have a look at temporal data. Display the timeseries of precipitation:

aguila precip.tss

It contains precipitation (mm/h) for 10 hours. It was used as input to a stochastic model calculating actual infiltration for each time step, using the realizations of infiltration capacity.

The folders 1 to 12 contain realizations of actual infiltration (mm/h), numbered from 1 to 10 (i0000000.001, i0000000.002, etc), i.e. for each time step of an hour one map. Check this.

To animate these realizations, type in the visualisation folder:

aguila veg.map --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4  --timesteps=[1,10,1] i

Plot a timeseries by right clicking on the legend of the infiltration maps i.

Question: By clicking on cells in the ‘Dense herbs’ area, have a look at the timeseries for this region. When is the uncertainty largest, with large amounts of rain or with low amounts of rain, or does rain not affect the uncertainty? Can you explain this?

  1. The uncertainty is largest for low amounts of rainfall. This is because for these cases the variation in rainfall is largest.

  2. The uncertainty is largest for high amounts of rainfall. This is because for these cases the variation in rainfall is largest.

  3. The uncertainty does not vary with rainfall amount.

  4. The uncertainty is largest during peak rainfall. With large amounts of rain, the differences among realizations of infiltration become clear.

Correct answers: d.

Feedback: None

It is also possible to plot maps over each other, e.g. type:

aguila --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4  --timesteps=[1,10,1] q + ldd.map

The maps q contain discharge (m3/h). Zoom in to see the content of the map. Click on the main channel and plot a timeseries. Animate.

3.1.2. Probability distributions and exceedance probabilities

Probability distributions are stored as a set of percentile maps. Unlike realizations, these are stored in the main folder visualisation. Type:

dir ksat_*map

Or in linux you would use ls. Note that * is a wildcard, representing all possible characters. These are maps containing the percentile values. The percentiles are given as fractional values between zero and 1 in the file name. You can display these maps just like other maps, e.g.:

aguila ksat_0.5.map

However, it is also possible to display them at once, resulting in a plot of the probability density distribution:

aguila --quantiles=[0.025,0.975,0.005] ksat veg.map

This opens all ksat percentile maps (and the vegetation map), and displays the probability distribution between the 2.5% and 97.5% percentiles, with a step of 0.5%. Note that, as the information is only available at a limited number of percentile values, aguila interpolates. You can open the probability density view by right clicking on the legend. Read through the Aguila manual at https://pcraster.geo.uu.nl/pcraster/latest/documentation/pcraster_aguila/Views.html#probability-graph-view for an explanation.

Question: Give the boundaries of the 95% confidence interval (i.e., between 2.5 and 97.5% percentiles) of saturated conductivity in the dense herbs area.

  1. 2.5%: 0.0345 m/h; 97.5%: 0.0694 m/h

  2. 2.5%: 0.0485 m/h; 97.5%: 0.0485 m/h

  3. 2.5%: 0.0485 m/h; 97.5%: 0.0694 m/h

  4. 2.5%: 0.0694 m/h; 97.5%: 0.0485 m/h

Correct answers: c.

Feedback: None

You can retrieve cumulative probabilities by clicking on the + in the menu at the top of the probability density distribution view. After this, you can slide the vertical line in the probability density view, to get cumulative probabilities for different threshold values. If exceedance intervals are needed, right-click on the legend, select Edit draw properties, and select ‘Exceedance probabilities’.

Question: Give the probability that saturated conductivity exceeds 0.056 in the dense herbs area.

  1. The probability is 0.71

  2. The probability is 0.29

  3. The probability is 0.0071

  4. The probability is 0.0029

Correct answers: a.

Feedback: None

Probability density distributions can also be displayed for dynamic data.

dir i_*

These are the percentile maps of modelled actual infiltration, for each time step of an hour. Filenames refer to time steps and the percentile. To display and animate these, type:

aguila --quantiles=[0.025,0.975,0.005] --timesteps=[1,10,1] i veg.map

By right-clicking on the legend of i, open timeseries and the probability distribution plot. Animate. Click on different locations of the map to evaluate the uncertainty in infiltration.

In a similar way, create an animation of the probability distribution of the discharge, q.

Question: Change the map view such that it shows the cumulative probability for a threshold value of 20000 m3/h. Note that the probability that q exceeds the threshold is one minus the cumulative probability. Animate the cumulative probability maps. Describe the pattern in the exceedance probabilities (in space and time).

  1. The higher the discharge, the lower the exceedance probability.

  2. The higher the discharge, the higher the exceedance probability.

Correct answers: b.

Feedback: None

3.1.3. Confidence intervals

To plot confidence intervals, type:

aguila --quantiles=[0.025,0.975,0.005] ksat veg.map

Open the probability density plot. Click on the ‘+’ in the probability density plot window to switch to cumulative probabilities. Now, right click on the legend, select Edit draw properties, and in the colour assignment panel, select ‘Confidence interval’. If you want you can change the alpha level of the confidence interval in the field below. Click ‘OK’. Now, move the threshold value in the probability density plot (the vertical line), and Aguila interactively shows the confidence interval for that threshold value.

Question: Which vegetation units are certainly (i.e. with an alpha value of 0.05) above a threshold value of 0.01 ?

  1. All units.

  2. Pine forest, dense herbs, and open shrubs.

  3. Open shrubs.

  4. All units except open shrubs.

Correct answers: c.

Feedback: None