3.3. Writing to disk, drawing realizations

3.3.1. Writing realizations to disk

Just like in deterministic modelling, report can be used to write data to disk.

In the framework folder, Open stochDynamicMod.py and save it as real.py. Add the following lines to the premcloop:

a = normal(1)
self.report(a,'a')

Add the same two statements to the initial and dynamic, using variable names b and c, respectively (instead of a). Save the script and run. Use dir to check the contents of the main folder framework and the folders containing the realizations, e.g. 1. Use aguila to display a, b, and c.

---

Question: What is stored by reporting in the dynamic?

1. It stores files in subdirectories referring to the Monte Carlo realizations, where each file has a number referring to the time step. Each file is a time series.

2. It stores files in subdirectories referring to the time steps, where each file has a number referring to the Monte Carlo realization. Each file is a time series.

3. It stores files in subdirectories referring to the Monte Carlo realizations, where each file has a number referring to the time step. Each file is a map.

4. It stores files in subdirectories referring to the time steps, where each file has a number referring to the Monte Carlo realization. Each file is a map.

Correct answers: c.

Feedback:

from pcraster import *
from pcraster.framework import *

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

  def premcloop(self):
    a = normal(1)
    self.report(a,'a')

  def initial(self):
    b = normal(1)
    self.report(b,'b')

  def dynamic(self):
    c = normal(1)
    self.report(c,'c')

  def postmcloop(self):
    print('running postmcloop')

nrOfSamples = 10
nrOfTimeSteps = 3
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel, nrOfTimeSteps)
mcModel = MonteCarloFramework(dynamicModel, nrOfSamples)
mcModel.run()

---

Opening files from disk can be done with self.readmap, just like in dynamic modelling.

3.3.2. Drawing realizations

The function normal draws realizations from a stochastic variable. PCRaster includes a number of other functions that draw realizations from stochastic variables, e.g. mapnormal, uniform. These were explained in the dynamic modelling exercises. Here, you will learn to use areanormal.

Open the map regions.map. Boolean True cells are ski regions. For the ski regions, vegetation height is a stochastic variable with a mean 0.2 and a standard deviation of 0.05 (m). For the remaining area, the mean is 0.3 with a standard deviation 0.1.

Open real.py and add statements (in the initial) to create realizations of vegetation height using areanormal. Use the variable name vegetationHeight and write to disk as h. Save the script and run. Use aguila to display 10 realizations of v. Hints:

• Read regions.map from disk in the premcloop. It is a deterministic map thus needs to be opened there.

• Use ifthenelse to create factors for mean and standard deviation.

---

Question: What is done by the areanormal function?

1. It draws a random value for each cell independently.

2. It draws a random value for each area independently, these are different between Monte Carlo realizations.

3. It normalizes input values within the boundaries of a normal distribution.

4. It draws a random value for each area independently, these are the same over the Monte Carlo realizations.

Correct answers: b.

Feedback:

from pcraster import *
from pcraster.framework import *

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

  def premcloop(self):
    a = normal(1)
    self.report(a,'a')
    self.regions=self.readmap('regions')

  def initial(self):
    b = normal(1)
    self.report(b,'b')

    normalAreas=areanormal(self.regions)
    sdFactor=ifthenelse(self.regions,0.05,scalar(0.1))
    meanFactor=ifthenelse(self.regions,0.2,scalar(0.3))
    vegetationHeight=normalAreas*sdFactor+meanFactor
    self.report(vegetationHeight,'h')

  def dynamic(self):
    c = normal(1)
    self.report(c,'c')

  def postmcloop(self):
    print('running postmcloop')

nrOfSamples = 10
nrOfTimeSteps = 3
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel, nrOfTimeSteps)
mcModel = MonteCarloFramework(dynamicModel, nrOfSamples)
mcModel.run()

---