3.6. Static stochastic modelling: neighbourhood operations¶

3.6.1. Drawing realizations of a stochastic digital elevation model¶

In the previous exercise, you used the following equation to calculate the velocity of forest fire spreading at a given cell, given in hours per metre. It was calculated as:

h = be^s/a

with, b, a parameter, s, the gradient of the topographical surface (m/m), a a parameter, and h velocity of forest fire spreading at the cell.

We assumed only a was unknown. However, as the digital elevation model will include error, we need to propagate the error in the digital elevation model to the output of the model h. Let’s assume the error in the elevation (m) is for each pixel modelled by a stochastic variable D:

D = e + Z

with, Z a stochastic variable with zero mean and a variance of 4.0.

In the folder fire, open stochStaticMod.py and save as slope.py. First, calculate realizations of D in the initial. In the premcloop, you still need to read the deterministic dem from disk. However, you need to add self. as the map is needed in the initial. So, the premcloop should look like this:

self.dem=self.readmap('dem')
self.gradient=slope(self.dem)
self.report(self.gradient,'gradient')

In the initial, you need to calculate the realizations by adding realizations of Z to the deterministic dem. As we need realizations on a cell-by-cell basis, you need to use normal. For D, use a variable name demStoch and write to disk as dS. Save the script and run.

Use Aguila to display 12 realizations of dS. Are the results approximately correct?

Now, calculate mean, variance, and percentiles of dS. Save the script and run. Display the variance map and the probability distribution plots.

Question: What is calculated by the script for the variance of dS? Is it similar to the value of 4.0 provided?

On average it is about 4, but it ranges between 3 and 5 depending on the cell location.
On average it is about 4, but it ranges between 1 and 7 depending on the cell location.
On average it is about 4, but it ranges between 3.9 and 4.1 depending on the cell location.
On average it is about 4, but it ranges between 3.99 and 4.01 depending on the cell location.

Correct answers: b.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(StaticModel, MonteCarloModel):
  def __init__(self):
    StaticModel.__init__(self)
    MonteCarloModel.__init__(self)
    setclone('clone.map')

  def premcloop(self):
    self.dem=self.readmap('dem')
    self.gradient=slope(self.dem)
    self.report(self.gradient,'gradient')

  def initial(self):
    demStoch=self.dem + normal(1) * 2.0
    self.report(demStoch,'dS')

    a=max(0.5 + mapnormal()/10,0.0001)
    self.report(a,'a')

    hoursPerMetre=0.002*exp(self.gradient/a)
    self.report(hoursPerMetre,'hpm')

  def postmcloop(self):
    names=['a','hpm','dS']
    sampleNumbers=self.sampleNumbers()
    print(sampleNumbers)
    timeSteps=[0]
    mcaveragevariance(names,sampleNumbers,timeSteps)
    percentiles=[0.025,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.975]
    mcpercentiles(names,percentiles,sampleNumbers,timeSteps)

nrOfSamples = 100
myModel = MyFirstModel()
staticModel = StaticFramework(myModel)
mcModel = MonteCarloFramework(staticModel, nrOfSamples)
mcModel.run()

3.6.2. Propagate DEM uncertainty through slope calculation¶

Using the stochastic digital elevation model defined in the previous section, you can now calculate the effect of this uncertainty on the uncertainty in the slope calculation. Thus far, slope was calculated in the premcloop. You need to move this calculation to the initial, using the realizations of the digital elevation model as input.

Open slope.py and make the required modifications. Be sure to calculate mean, variance, and percentiles of the slope. Save the script and run. Display the slope results to check your calculations.

Question: The uncertainty in fire spreading velocity of includes now the uncertainty of the parameter a and the uncertainy of the digitial elevation model. Display house.map and the variance map of the forest fire spreading velocity (h/m). What is the variance and standard deviation of the forest fire spreading velocity (m/h) at the location of the house? Compare the result with the run where we only included error in the a parameter (previous exercise). What is the contribution to uncertainty in forest fire spreading of the uncertainty of the digitial elevation model, approximately?

The uncertainty becomes zero.
The uncertainty has not changed (changes are due to numerical error).
The uncertainty is much lower now.
The uncertainty is much higher now.

Correct answers: d.

Feedback:

The uncertainty is much bigger. The variance is 1.696 x 10-6 and the standard deviation is 1.3 x 10-3.

from pcraster import *
from pcraster.framework import *

class MyFirstModel(StaticModel, MonteCarloModel):
  def __init__(self):
    StaticModel.__init__(self)
    MonteCarloModel.__init__(self)
    setclone('clone.map')

  def premcloop(self):
    self.dem=self.readmap('dem')

  def initial(self):
    demStoch=self.dem + normal(1) * 2.0
    self.report(demStoch,'dS')

    gradient=slope(self.dem)
    self.report(gradient,'gradient')

    a=max(0.5 + mapnormal()/10,0.0001)
    self.report(a,'a')

    hoursPerMetre=0.002*exp(gradient/a)
    self.report(hoursPerMetre,'hpm')

  def postmcloop(self):
    names=['a','hpm','dS', 'gradient']
    sampleNumbers=self.sampleNumbers()
    print(sampleNumbers)
    timeSteps=[0]
    mcaveragevariance(names,sampleNumbers,timeSteps)
    percentiles=[0.025,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.975]
    mcpercentiles(names,percentiles,sampleNumbers,timeSteps)

nrOfSamples = 100
myModel = MyFirstModel()
staticModel = StaticFramework(myModel)
mcModel = MonteCarloFramework(staticModel, nrOfSamples)
mcModel.run()

3.6.3. Neighbourhood defined by topology: relative distance calculation¶

With the fire spreading velocity (h/m) you can calculate a map with the time it takes until the fire reaches a certain location, given a specified starting location of the fire. We oversimplify here somewhat, by assuming that fire burns only in uphill direction. The ldddist function can be used to calculate a relative distance, i.e. the absolute distance (m) that is travelled through a cell multiplied by a relative ‘weighting’ for that cell, here the fire spreading velocity (h/m). The result is the time (h) for the fire to reach a location. It restricts the distances calculated to uphill local drain directions.

Have a look at the description of the function ldddist in the PCRaster manual. Have a look at start.map, it is the starting point of a fire.

Open slope.py and save as fire.py. Make the following additions or changes:

Read start.map from disk in the premcloop.
Calculate the local drain direction map in premcloop and write to disk.
In initial, use ldddist to calculate a map with the time the fire takes to reach a cell. As inputs, you need the local drain direction map, and the starting point of the fire, and the fire velocity map in hours per metre. Name the resulting map time and write to disk with report as time.
In the postmcloop, calculate sampling statistics (mean, variance, percentiles) of fire.

Save the script and run. Display the mean, variance, and the probability density function of time, together with the local drain direction map and house.map.

Question: What is the mean and standard deviation of the time it takes until the fire reaches the house?

Mean: 2.5 h. Standard deviation: 0.30 h.
Mean: 2.5 h. Standard deviation: 0.92 h.
Mean: 2.5 h. Standard deviation: 3.7 h.
Mean: 4.5 h. Standard deviation: 0.37 h.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(StaticModel, MonteCarloModel):
  def __init__(self):
    StaticModel.__init__(self)
    MonteCarloModel.__init__(self)
    setclone('clone.map')

  def premcloop(self):
    self.dem=self.readmap('dem')
    self.ldd=lddcreate(self.dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')
    self.startPoint=self.readmap('start')

  def initial(self):
    demStoch=self.dem + normal(1) * 2.0
    self.report(demStoch,'dS')

    gradient=slope(self.dem)
    self.report(gradient,'gradient')

    a=max(0.5 + mapnormal()/10,0.0001)
    self.report(a,'a')

    hoursPerMetre=0.002*exp(gradient/a)
    self.report(hoursPerMetre,'hpm')

    time=ldddist(self.ldd,self.startPoint,hoursPerMetre)
    self.report(time,'time')

  def postmcloop(self):
    names=['a','hpm','dS', 'gradient', 'time']
    sampleNumbers=self.sampleNumbers()
    print(sampleNumbers)
    timeSteps=[0]
    mcaveragevariance(names,sampleNumbers,timeSteps)
    percentiles=[0.025,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.975]
    mcpercentiles(names,percentiles,sampleNumbers,timeSteps)

nrOfSamples = 100
myModel = MyFirstModel()
staticModel = StaticFramework(myModel)
mcModel = MonteCarloFramework(staticModel, nrOfSamples)
mcModel.run()

Question: Assume that it takes two hours (after the initiation of the fire) until the owners of the house can be informed they should leave the house. Create the confidence interval plot for a threshold of two hours, using an alpha value of 0.1. Can we be certain with this model that the owners can be informed in time?

Yes. The alpha value is below 2 hours, so we do even not need to run the model.
Well, it is unclear. The fire might be there within two hours, but it can also take longer.
No, according to the model we can be certain that the fire takes less than two hours to reach the house. The owners cannot be informed in time.
Yes, according to the model we can be certain that the fire takes more than two hours to reach the house. The owners can be informed in time.

Correct answers: d.

Feedback: none