The model development cycle consists typically of two steps: the conversion of the conceptual model structure into a numerical implementation, and the assimilation or calibration of the model with observed data (Karssenberg, in prep). Most of the available tools are either dedicated to model development (e.g. ModelMaker, Stella, Matlab), or to calibration or assimilation (e.g. PEST, ReBEL). The desired combination of those tools usually requires an additional software layer.
This framework offers a combined interface for the task of model construction and optimisation. Spatial operations can be used from the PCRaster library – temporal, deterministic and stochastic modelling capabilities are added with this framework. Modelling options are Monte Carlo simulations (e.g. Doucet et al, 2001), Particle filter (e.g. Xiong, 2006) and Ensemble Kalman filter (e.g Simon, 2006).
The example script shows a spatio-temporal model executed in a Monte Carlo simulation. The process descriptions simulating the snow pack of a catchment are given in the initial and dynamic sections of the SnowModel class. The Monte Carlo framework executes the premcloop to generate initial values valid for all samples, and afterwards 100 realisations each running 180 time steps. Finally the postmcloop calculates statistics of the ensemble results.
For a more detailed description of the model see Karssenberg (submitted) and the manual.
from pcraster import * class SnowModel(object): def __init__(self): setclone("clone.map") def premcloop(self): dem = self.readmap("dem") self.ldd = lddcreate(dem, 1e31, 1e31, 1e31, 1e31) elevationMeteoStation = scalar(2058.1) self.elevationAboveMeteoStation = dem - elevationMeteoStation self.degreeDayFactor = 0.01 def initial(self): self.snow = scalar(0) self.temperatureLapseRate = 0.005 + (mapnormal() * 0.001) self.report(self.temperatureLapseRate, "lapse") self.temperatureCorrection = self.elevationAboveMeteoStation * self.temperatureLapseRate def dynamic(self): temperatureObserved = self.readDeterministic("tavgo") precipitationObserved = self.readDeterministic("pr") precipitation = max(0, precipitationObserved * (mapnormal() * 0.2 + 1.0)) temperature = temperatureObserved - self.temperatureCorrection snowFall = ifthenelse(temperature < 0, precipitation, 0) self.snow = self.snow + snowFall potentialMelt = ifthenelse(temperature > 0, temperature * self.degreeDayFactor, 0) actualMelt = min(self.snow, potentialMelt) self.snow = max(0, self.snow - actualMelt) rain = ifthenelse(temperature >= 0, precipitation, 0) discharge = accuflux(self.ldd, actualMelt + rain) self.report(self.snow, "s") self.report(discharge, "q") def postmcloop(self): names = ["s", "q"] mcaveragevariance(names, self.sampleNumbers(), self.timeSteps()) percentiles = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] mcpercentiles(names, percentiles, self.sampleNumbers(), self.timeSteps()) myModel = snowModel() dynamicModel = DynamicFramework(myModel, 180) mcModel = MonteCarloFramework(dynamicModel, 100) mcModel.run()
The image below shows a comparison of model results simulated in Monte Carlo mode (upper row) and with the particle filter (lower row, filter moments at time steps 70, 100 and 150). The left images show the snow pack (m) results of 10 realisations, the centre image the map view with currently selected cell and the right images show the cumulative probability density function.
Distribution of software
The software is part of the current release of PCRaster.
This research was supported by the Dutch research programme Space for Geo-Information (Ruimte voor Geo-Informatie, RGI) as RGI project 313 (“Online coupling of spatial optimisation tools and spatially distributed simulation models”) with the cooperation partners PCRaster Environmental Software, Department of Physical Geography, Utrecht University and the Joint Research Centre.
Doucet, A., de Freitas, N. and Gordon, N. (2001):
Sequential Monte Carlo Methods in Practice. Springer
Xiong, X., Navon, I. M. and Uzunoglu, B. (2006):
A note on the particle filter with posterior Gaussian resampling. Tellus Series A: Dynamic Meteorology and Oceanography
Simon, D. (2006):
Optimal State Estimation: Kalman, H and Nonlinear Approaches. Wiley Interscience
Karssenberg, D. et al, submitted:
A software framework for construction of stochastic spatio-temporal models assimilated or calibrated with observational data.