2.8. Probabilistic spatial models: forest fire and seed dispersal¶

2.8.1. Random variables in point models¶

2.8.1.1. Continuous random variables¶

PCRaster comes with a number of functions to draw realizations from continuous probability distributions.

Go to the folder /probab/fire in your exercises data set and open randomVar.py To the dynamic, add a line that draws a realization using the function normal Store the result. In a similar way, use uniform and store the result. Animate the outputs.

Question: What is calculated by the function uniform?

It draws realizations resulting in a uniform average value of 1.0.
It draws realizations from a uniform distribution between zero and one, independently for each cell and time step.
It draws a realization from a uniform distribution between zero and one, assigning it to all cells, independently for each time step.

Correct answers: b.

Feedback:

Look up the description of the function in the PCRaster manual pages if you do not understand the function.

Now, test the function mapuniform by adding it to the dynamic section and storing results. Look it up in the PCRaster manual pages if you need to find what inputs it requires.

Question: In what situations would you use uniform and in what situations would you use mapuniform?

The function uniform is used when error at a cell is independent from other cells. The function mapuniform represents uncertainty of an error source that is fully correlated over the area, i.e. over the whole area the error causes the value to be too high or too low; with uniform, the value would be too high or too low for each cell independently.
The function mapuniform is used when error at a cell is independent from other cells. The function uniform represents uncertainty of an error source that is fully correlated over the area, i.e. over the whole area the error causes the value to be too high or too low; with mapuniform, the value would be too high or too low for each cell independently.
The function uniform is used when the error is larger than indicated by the measurement instrument, as it adds uncertainty caused by spatial variation. The function mapuniform is used when the measurement instrument exactly represents the error over the whole area.

Correct answers: a.

Feedback:

None.

Now, add a statement that draws realizations from a Gaussian distributed variable with mean a mean of 10 and standard deviation of 5.

Question: What is the range of the values in the realizations that you find if you run the model for all time steps?

The values are approximately between 5.0 and 10.0
The values are approximately between 0.0 and 20.0
The values are approximately between 5.0 and 15.0
The values are approximately between -15.0 and 35.0

Correct answers: d.

Feedback:

The result should be realizations for each cell from a Gaussian distribution with a mean 10 and a standard deviation of 5. In theory, any value could be drawn (as the Gaussian distribution is not bounded), but with this number of cells and time steps, one find mostly values in the range between -15 and 35. If you did not find this, you most likely made an error in the script. In the dynamic, the following two lines are needed:
    # draw from a normal distribution
    # with mean zero and standard deviation one
    nor=normal(1)
    # adjust the realizations:
    # to get a standard deviation of five (instead of one),
    # we first multiply by five, which will change the sd,
    # then we add 10.0, to move the mean to ten.
    rv=10.0+nor*5.0
    self.report(rv,'rv')

2.8.1.2. Discrete random variables¶

Functions to draw realizations of discrete (classified) random variables are not provided. These can be derived from the functions for continuous random variables that you used in the previous section.

Open randomVar.py and add (a) statement(s) to the script that draws realizations (for each cell independently and for each timestep) of a discrete random variable with two possible outcomes, True and False:

P(True) = 0.8

P(False) = 0.2

Write results to disk, run the script, and display to check the result.

Question: How could you create the discrete random variable?

By drawing realizations from a uniform distribution between 0.2 and 0.8, converting the result to a Boolean.
By drawing realizations from a uniform distribution between 0.2 and 0.8, selecting the values above 0.5.
By drawing realizations from a normal distribution with a mean of 0.5 and a standard deviation of 1.0, selecting values above 0.2.
By drawing realizations from a uniform distribution between 0.0 and 1.0, selecting values below 0.8.

Correct answers: d.

Feedback:

    uni=uniform(1)
    aOne=uni < 0.8
    self.report(aOne,'aone')

Now, add statements that draw realizations of a discrete random variable with three possible outcomes, coded as 1, 2, and 3:

P(1) = 0.01
P(2) = 0.9
P(3) = 0.09

Question: How could you create the discrete random variable?

Using realizations from a uniform distribution, a table, and lookupnominal.
By multiplication of the outcome of the function uniform, three times.
Using realizations from a uniform distribution, and the maptotal function, to select values.
Using a window function.

Correct answers: a.

Feedback:

Just like in the previous exercise, we will derive the realizations from realizations of a uniform distribution between zero and one. First you need to create an ascii table to assign the correct values, save it as threeCl.tbl, for instance:
[0,0.01> 1
[0.01,0.91> 2
[0.91,1.0] 3
Then, use lookup to select the values from the uniform distribution:
    uni=uniform(1)
    classes=lookupnominal('threeCl.tbl', uni)
    self.report(classes,'classes')

2.8.2. Direct neighbourhood interaction: forest fire¶

2.8.2.1. Fire spreading¶

Direct neighbourhood functions (as used in cellular automata) are suitable for modelling forest fire. This is because the large scale spatio-temporal pattern of an evolving fire is driven by local interactions at the fire front. These local interactions will be represented here using a stochastic direct neighbourhood model. The model uses the following transition rules that are applied for each time step:

A cell that was burning at the previous time step stays burning.
Cells that are not burning at the previous time step potentially catch fire only when one or more neighbouring cells are burning at the previous time step. Neighbouring cells are the four direct neighbours in the x, y direction.
The cells under item 2 actually catch fire with a probability of 0.1.
A cell that was not burning at the previous time step and actually catched fire in the current timestep becomes a burning cell.

Display dem.map and start.map. The start.map contains True cells representing the starting area of a fire. Open fire.py. First add statements to represent rule 2.

Set the initial distribution of the fire (use the variable name fire) to start.map.
In the dynamic, add a window4total statement that calculates the number of neighbours that are burning.
In the dynamic, use the result of window4total to calculate a map neighbourBurns that is True when one or more neighbouring cells are burning.
Write the resulting maps to disk, save the script, and run.

Animate neighbourBurns. Do the cells with True also include cells that were already burning in the previous time step?

And finally for rule 2, add (a) statement(s) that calculate(s) potentialNewFire, a map that is True only for cells that potentially catch fire. To calculate this map, you need to combine fire and neighbourBurns using Boolean logic operations. Write to disk and check the result.

Question: How did you calculate potentialNewFire?

Using Boolean logic operations.
Using multiplication.
Using an additional window4total statement.
Using an additional window4total statement and areaaverage.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class Fire(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    self.fire=self.readmap('start')

  def dynamic(self):
    # nr of neighbours with fire
    nrOfBurningNeighbours=window4total(scalar(self.fire))
    # cells where at least one neighbour is burning
    neighbourBurns=nrOfBurningNeighbours > 0
    self.report(neighbourBurns,'nb')

    # cells that are not yet burning and where one neighbour burns
    potentialNewFire=neighbourBurns & ~ self.fire
    # if you want to avoid the use of operators the previous line
    # can also be written as follows:
    # potentialNewFire=pcrand(neighbourBurns,pcrnot(self.fire))
    
    self.report(potentialNewFire,'pnf')

nrOfTimeSteps=200
myModel = Fire()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

Now, let’s add rule 3.

Add the following statement to the dynamic calculating for each time step a map that is True with a probability of 0.1 and false elsewhere. Write the result to disk to check it.

realization = uniform(1) < 0.1

Add a statement that calculates newFire, a map that is True for cells that actually catch fire in the current time step. newFire represents rule 3. Write to disk, run, and check the result.

And finally, update fire by combining newFire and fire from the previous time step. Write to disk, store your model script, run it, and check the results.

Question: How does the fire evolve over time?

Due to the window shape of the window4total function it is covers an approximately rectangular block which increases in size over time, having a ragged edge due to the random process.
It covers a circle which increases in size, ragged at the edge due to the random process.
It covers a circle which increases in size, with a smooth edge.
Due to the window shape of the window4total function it is covers an approximately rectangular block which increases in size over time, having a smooth edge.

Correct answers: b.

Feedback:

from pcraster import *
from pcraster.framework import *

class Fire(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    self.fire=self.readmap('start')

  def dynamic(self):
    # nr of neighbours with fire
    nrOfBurningNeighbours=window4total(scalar(self.fire))
    # cells where at least one neighbour is burning
    neighbourBurns=nrOfBurningNeighbours > 0
    self.report(neighbourBurns,'nb')

    # cells that are not yet burning and where one neighbour burns
    potentialNewFire=neighbourBurns & ~ self.fire
    # if you want to avoid the use of operators the previous line
    # can also be written as follows:
    # potentialNewFire=pcrand(neighbourBurns,pcrnot(self.fire))
    
    self.report(potentialNewFire,'pnf')

    realization=uniform(1) < 0.1
    self.report(realization,'real')

    newFire=potentialNewFire & realization
    # or alternative notation:
    # newFire=pcrand(potentialNewFire,realization)
    self.report(newFire,'nf')

    self.fire=self.fire | newFire
    # or:
    # self.fire=pcror(self.fire,newFire)
    self.report(self.fire,'fire')

nrOfTimeSteps=200
myModel = Fire()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

2.8.2.2. Fire spreading determined by the topography¶

In the previous exercise, we assumed a probability of a cell catching fire to be everywhere the same, i.e. 0.1. In reality, however, a fire front mainly extends uphill. A very simple representation of this process is used here. The transition rules are kept the same, except for the probability used in rule 2. Now we will calculate this probability as:

All cells with a direct neighbour in downhill direction have a probability to catch fire of 0.8.
For all other cells, this probability is 0.1.

Open the script fire.py that was created in the previous exercise and save it as fireUphill.py. In the initial, create the map ldd using the digital elevation model dem.map. You need self.ldd as variable name as it will be used in the dynamic. Display the local drain direction map on top of the digital elevation model, zoom in if you get an almost black map!

In the dynamic, you need to add statements calculating for each map the probability for cells to catch fire. As a start, add two statements (insert them directly above the calculation of realization:

A statement calculating downhillBurning, a map that is True for cells with a cell in downhill direction that is burning. Use the function downstream.
A statement calculating probability, a map with the probabilities that a cell catches fire. Use downhillBurning as input.

Let the model write the results of these two calculations to disk. Save and run the model to see the results.

Now, you need to change the calculation of realization as it needs to take into account the spatial pattern (for each time step) of the probability. Do this by changing the statement. Save and run the model. Animate the resulting fire pattern together with the local drain direction map and the digital elevation model.

Question:: What is the pattern of the fire spreading over time? Do you think this is more realistic compared to the previous model? Give two suggestions to further improve the model.

Question: Have a look at the veg.map. At what time step does the fire reach the pine forest in the southern part of the area (approximately)?

At time step 50.
At time step 100.
At time step 150.
It does not reach the pine forest.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class Fire(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    self.fire=self.readmap('start')
    #
    dem=self.readmap('dem')
    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)

  def dynamic(self):
    nrOfBurningNeighbours=window4total(scalar(self.fire))
    neighbourBurns=nrOfBurningNeighbours > 0
    self.report(neighbourBurns,'nb')
    potentialNewFire=neighbourBurns & ~ self.fire
    self.report(potentialNewFire,'pnf')

    # 
    downhillBurning=downstream(self.ldd,self.fire)
    self.report(downhillBurning,'db')

    probability=ifthenelse(downhillBurning,scalar(0.8),0.1)
    self.report(probability,'prob')
    #

    #realization=uniform(1) < 0.1
    realization=uniform(1) < probability
    self.report(realization,'real')

    newFire=potentialNewFire & realization
    self.report(newFire,'nf')

    self.fire=self.fire | newFire
    self.report(self.fire,'fire')

nrOfTimeSteps=200
myModel = Fire()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

2.8.3. Neighbourhood interaction over a distance: plant seed dispersal¶

Plant dispersal occurs either as clonal growth or through seed dispersal. Clonal growth can be modelled using direct neighbourhood operations as new plants only occur directly neighbouring existing cells. Thus, it is actually a quite similar process, from the modelling point of view, as forest fire. Seed dispersal, however, allows plants to disperse over longer distances: new plants may establish also at cells that are not neighbours of cells containing plants in the previous time step. Here we will construct a simple seed dispersal model.

We use the following set of rules:

A single plant species disperses over an area.
For each time step, the probability P(E) that new plants establish at cells that are not yet occupied by the species is modelled as:

P(E) = 0.1^(d/r)

with, d, the distance away (m) from the nearest cell containing the plant, and r, a parameter (m, r = 40 m).

For each time step, the new distribution of the plants consists of the distribution of the plants in the previous time step combined with the distribution of the plants that established in the current time step.

Go to the directory probab/plants and display veg.map. We assume that the plant species to be modelled is artificially introduced (before dispersal starts) at all pixels containing dense herbs (class with code 5). Open plants.py. In the initial, create a Boolean map (use the variable name self.plants) that is True at dense herbs and false elsewhere. This is the map containing the plants that will disperse over time.

To represent the transition rules, add statements to the dynamic creating a map probability with P(E) for each time step. For now, you can assume that the plant does not yet spread (i.e., self.plants does not change), so probability will not change over time. Store the script, and run. Display the results.

Question: What is the probability in the left corner of the study area?

Below 0.01.
Between 0.01 and 0.8.
Above 0.8.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class Plants(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    vegetation=readmap('veg')
    self.plants=vegetation == 5
    self.range=40

  def dynamic(self):
    distance=spread(self.plants,0,1)
    self.report(distance,'dist')

    probability=0.1**(distance/self.range)
    self.report(probability,'prob')

    
nrOfTimeSteps=50
myModel = Plants()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

Using probability, create a new map newPlants that is True at cells where the plant establishes and false elsewhere. You need to use probability as input, and the uniform function. Write the result to disk and check the outcome.

Now, update self.plants by ‘adding’ the newly established plants to self.plants. Write self.plants to disk in the dynamic. Save the script, run, and animate the plants dispersal, together with the probability maps. Run the model two times and compare results.

Question: How long does it approximately take to reach the the bottom left side of the area?

5 time steps.
10 time steps.
25 time steps.
40 time steps.

Correct answers: c.

Feedback:

from pcraster import *
from pcraster.framework import *

class Plants(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    vegetation=readmap('veg')
    self.plants=vegetation == 5
    self.range=40

  def dynamic(self):
    distance=spread(self.plants,0,1)
    self.report(distance,'dist')

    probability=0.1**(distance/self.range)
    self.report(probability,'prob')

    newPlants=probability > uniform(1)  
    self.report(newPlants,'np')

    self.plants=pcror(self.plants, newPlants)

    self.report(self.plants,'plants')

    
nrOfTimeSteps=50
myModel = Plants()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

We can make all kinds of interesting extensions to the model. Let’s see what happens if we assume the plant can only grow at vegetation units Open Shrub (code 3) and Dense herbs (code 5). In the initial, create a map biotope that is True at cells containing Open shrub or Dense herbs. In the dynamic, modify the script such that the plant will only exist or establish at these open shrub or dens herbs cells. Save the script, run.

Question: How long does it approximately take to reach the Open shrub patch on the left side of the area? Explain.

It will not be reached for sure as the probabilities are too low.
Almost certainly within 10 years.
It might be reached although probabilities are extremely low.
Almost certainly within 10 to 50 years.

Correct answers: c.

Feedback:

from pcraster import *
from pcraster.framework import *

class Plants(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('clone.map')

  def initial(self):
    vegetation=readmap('veg')
    self.plants=vegetation == 5
    self.range=40

    self.biotope=pcror(vegetation == 5, vegetation == 3)
    self.report(self.biotope,'bio')

  def dynamic(self):
    distance=spread(self.plants,0,1)
    self.report(distance,'dist')

    probability=0.1**(distance/self.range)
    self.report(probability,'prob')

    newPlants=probability > uniform(1)  
    self.report(newPlants,'np')

    self.plants=pcror(self.plants, newPlants)
    self.plants=pcrand(self.plants,self.biotope)

    self.report(self.plants,'plants')

    
nrOfTimeSteps=50
myModel = Plants()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()