2.6. Calculating descriptive statistics: the snow melt model

2.6.1. Over the spatial domain

For analysis or application of model outputs, it is often needed to calculate statistical values (e.g., mean, maximum value) over regions or time periods. We focus first on regions.

Display the map skiregio.map. It contains two skiing regions in the study area.

Open the script runoff.py and save as ski.py. For evaluation of the regions for possibilities for skiing, calculate the average snowpack (m water equivalent) in each of these regions, for each time step. Use the function areaaverage. Write the results to disk and display. Plot a time series (right-click on the legend) and click on one of the regions to get the time series for that particular region.

---

Question: Which of the two regions has in general the largest amount of snow?

1. Region 1.

2. Region 2.

Correct answers: b.

Feedback: None

---

Although average snow pack is important, it is more important that the snowpack thickness is above a certain threshold required for skiing. Calculate for each time step:

• A map of the study area that is True at cells with a snowpack that is thicker than 0.2 m (water equivalent thickness), and

• A map with the proportion (a value between zero and one) of each area that has a snowpack that is thicker than 0.2 m.

For the second step, convert the map created in the first step to scalar data type using scalar. Now you can use it in an areatotal function. To get the fraction, you need the number of cells for each region, which can be calculated by:

areaarea(skiregio.map)/cellarea()

Store under the same filename ski.py and run.

---

Question: Compare the two regions by plotting timeseries. For approximately how many days is each region fully covered with a snow pack that is sufficiently thick for skiing ?

1. Based on the results we can conclude that there is never sufficient snow in both areas.

2. 5 days.

3. 10 days.

4. 40 days.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow=0.0

    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')
    self.skiRegions = self.readmap('skiregio')
    self.numberCellsRegion = areaarea(self.skiRegions)/cellarea()

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'freez')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

    self.snow = self.snow+snowFall

    potentialMelt = ifthenelse(pcrnot(freezing),temperature*0.01,0)
    self.report(potentialMelt,'pmelt')
    actualMelt = min(self.snow, potentialMelt)
    self.report(actualMelt,'amelt')

    self.snow = self.snow - actualMelt

    runoffGenerated = actualMelt + rainFall
    self.report(runoffGenerated,'rg')

    discharge=accuflux(self.ldd,runoffGenerated*cellarea())
    self.report(discharge,'q')

    self.report(self.snow,'snow')

    snowInSkiRegion = areaaverage(self.snow,self.skiRegions)
    self.report(snowInSkiRegion,'avsn')
    thickEnough = ifthenelse(self.snow > 0.2, boolean(1), boolean(0))
    self.report(thickEnough,'th')
    numberCellsThickEnough = areatotal(scalar(thickEnough),self.skiRegions)
    self.report(numberCellsThickEnough,'nth')
    proportionThickEnough = numberCellsThickEnough / self.numberCellsRegion
    self.report(proportionThickEnough,'ski')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

2.6.2. Over time

It is also possible to calculate statistics over time.

Calculate a map with the maximum snowpack thickness (water equivalent) over time. Use the function max. Add the required statements to ski.py.

---

Question: What is the maximum snow thickness in the bottom right corner cell?

1. 0.30 m

2. 0.91 m

3. 0.27 m

4. 1.90 m

Correct answers: c.

Feedback:

In order to answer this question we need to add a new global variable to the initial:

self.maxSnowThickness = 0.0

At the end of the dynamic we add:

self.maxSnowThickness = max(self.snow,self.maxSnowThickness)
self.report(self.maxSnowThickness,'max')

If we visualize the last time step (aguila max00000.181) we can read the result: 0.27 m.

---