# -*- coding: utf-8 -*-
'''
Copyright 2015 by Tobias Houska
This file is part of Statistical Parameter Estimation Tool (SPOTPY).
:author: Philipp Kraft
A one dimensional cmf model analysing data from the Schwingbach hydrological observatory.
You need to have cmf and pandas installed on your system: svn checkout svn://fb09-pasig.umwelt.uni-giessen.de/cmf/trunk cmf
'''
import cmf
from datetime import datetime, timedelta
import numpy as np
import spotpy
import os, sys
class _CmfProject:
"""
This class describes the cmf setup
"""
def __init__(self, par):
cmf.set_parallel_threads(1)
# run the model
self.project = cmf.project()
self.cell = self.project.NewCell(x=0, y=0, z=238.628, area=1000, with_surfacewater=True)
c = self.cell
r_curve = cmf.VanGenuchtenMualem(Ksat=10**par.pKsat, phi=par.porosity, alpha=par.alpha, n=par.n)
# Make the layer boundaries and create soil layers
lthickness = [.01] * 5 + [0.025] * 6 + [0.05] * 6 + [0.1] * 5
ldepth = np.cumsum(lthickness)
# Create soil layers and pick the new soil layers if they are inside of the evaluation depth's
for d in ldepth:
c.add_layer(d, r_curve)
# Use a Richards model
c.install_connection(cmf.Richards)
# Use shuttleworth wallace
self.ET = c.install_connection(cmf.ShuttleworthWallace)
c.saturated_depth = 0.5
self.gw = self.project.NewOutlet('groundwater', x=0, y=0, z=.9)
cmf.Richards(c.layers[-1], self.gw)
self.gw.potential = c.z - 0.5 #IMPORTANT
self.gw.is_source = True
self.gwhead = cmf.timeseries.from_scalar(c.z - 0.5)
self.solver = cmf.CVodeIntegrator(self.project, 1e-9)
def set_parameters(self, par):
"""
Sets the parameters of the model
:param par: The parameter value object (named tuple)
:return:
"""
try:
for l in self.cell.layers:
r_curve = cmf.VanGenuchtenMualem(Ksat=10**par.pKsat, phi=par.porosity, alpha=par.alpha, n=par.n)
r_curve.w0 = r_curve.fit_w0()
l.soil = r_curve
self.cell.saturated_depth = 0.5
self.gw.potential = self.cell.z - 0.5
except RuntimeError as e:
sys.stderr.write('Set parameters failed with:\n' + str(par) + '\n' + str(e))
raise
def load_meteo(self, driver_data):
datastart = driver_data.time[0]
# Create meteo station for project
meteo = self.project.meteo_stations.add_station('Schwingbach', position=(0, 0, 0), tz=1, timestep=cmf.h)
rain = cmf.timeseries.from_array(datastart, cmf.h, driver_data.rain_mmday)
meteo.rHmean = cmf.timeseries.from_array(datastart, cmf.h, driver_data.relhum_perc)
meteo.Windspeed = cmf.timeseries.from_array(datastart, cmf.h, driver_data.windspeed_ms)
meteo.Rs = cmf.timeseries.from_array(datastart, cmf.h, driver_data.solarrad_wm2 * 86400e-6)
meteo.T = cmf.timeseries.from_array(datastart, cmf.h, driver_data.airtemp_degc)
meteo.Tmax = meteo.T.floating_max(cmf.day)
meteo.Tmin = meteo.T.floating_min(cmf.day)
self.project.rainfall_stations.add('Schwingbach', rain, (0, 0, 0))
self.project.use_nearest_rainfall()
# Use the meteorological station for each cell of the project
self.project.use_nearest_meteo()
self.meteo = meteo
self.gwhead = cmf.timeseries.from_array(datastart, cmf.h, driver_data.gwhead_m)
def set_boundary(self, t):
"""
Sets the boundary conditions for time t
:param t: a cmf.Time
:return: None
"""
gw_level = self.gwhead[t]
if np.isfinite(gw_level):
self.gw.potential = gw_level
def _load_data(filename):
"""
Loads data from csv, where the first column has the measurement date and all the other columns additional data
:return: Rec-Array with date as a first column holding datetimes
"""
def str2date(s):
"""Converts a string to a datetime"""
return datetime.strptime(s.decode(), '%Y-%m-%d %H:%M:%S')
# Load the data
return np.recfromcsv(filename, converters={0: str2date}, comments='#')
class _ProgressReporter:
"""
Simple helper class to report progress and check for too long runtime
"""
def __init__(self, start, end, frequency=cmf.week, max_runtime=15 * 60, verbose=False):
self.verbose = verbose
self.frequency = frequency
self.max_runtime = max_runtime
self.stopwatch = cmf.StopWatch(cmf.AsCMFtime(start), cmf.AsCMFtime(end))
def __call__(self, t):
elapsed, total, remaining = self.stopwatch(t)
if self.verbose and not t % cmf.week:
print('{modeltime:%Y-%m-%d}: {elapsed}/{total}'
.format(modeltime=t.AsPython(),
elapsed=elapsed * cmf.sec,
total=total * cmf.sec))
if elapsed > self.max_runtime:
raise RuntimeError('{:%Y-%m-%d %H:%S} model took {:0.0f}min until now, stopping'
.format(t.AsPython(), elapsed / 60))
class Cmf1d_Model(object):
"""
A 1d Richards based soilmoisture model for Schwingbach site #24
"""
alpha = spotpy.parameter.Uniform(0.0001, 0.2, optguess=0.1156, doc=u'α in 1/cm for van Genuchten Mualem model')
pKsat = spotpy.parameter.Uniform(-2, 2, optguess=0, doc=u'log10 of saturated conductivity of the soil in m/day')
n = spotpy.parameter.Uniform(1.08, 1.8, optguess=1.1787, doc=u'van Genuchten-Mualem n')
porosity = spotpy.parameter.Uniform(0.3, 0.65, optguess=0.43359, doc=u'φ in m³/m³')
def __init__(self, days=None):
self.driver_data = _load_data('cmf_data/driver_data_site24.csv')
self.evaluation_data = np.loadtxt('cmf_data/soilmoisture_site24.csv', delimiter=',',
comments='#', usecols=[1,2,3])
self.datastart = self.driver_data.time[0]
# Set the end point for the model runtime, either by parameter days or
# run to the end of available driver data
if days is None:
self.dataend = self.driver_data.time[-1]
else:
self.dataend = self.datastart + timedelta(days=days)
# The depth below ground in m where the evaluation data belongs to
self.eval_depth = [0.1, 0.25, 0.4]
# Make the model
self.model = _CmfProject(self.make_parameters())
# Load meteo data
self.model.load_meteo(driver_data=self.driver_data)
self.__doc__ += '\n\n' + cmf.describe(self.model.project)
def make_parameters(self, random=False, **kwargs):
return spotpy.parameter.create_set(self)
def evaluation(self):
"""
:return: The evaluation soilmoisture as a 2d array
"""
return self.evaluation_data
def objectivefunction(self, simulation, evaluation):
"""
Returns the negative RMSE for all values
:param simulation:
:param evaluation:
:return:
"""
# Find all data positions where simulation and evaluation data is present
take = np.isfinite(simulation) & np.isfinite(evaluation)
rmse = -spotpy.objectivefunctions.rmse(evaluation=evaluation[take], simulation=simulation[take])
return rmse
def get_eval_layers(self, eval_depth=None):
"""
Returns the list of layers, that should be evaluated
:param eval_depth: List of depth below ground where evaluation should be done. If None, all layers are evaluated
:return:
"""
c = self.model.cell
if eval_depth:
edi = 0 # current index of the evaluation depth
eval_layers = [] # the layers to do the evaluation are stored here
for l in c.layers:
if edi < len(self.eval_depth) and l.upper_boundary <= self.eval_depth[edi] < l.lower_boundary:
eval_layers.append(l)
edi += 1
else:
eval_layers = list(c.layers)
return eval_layers
def simulation(self, par=None, verbose=False):
'''
Runs the model instance
:param par: A namedtuple of model parameters
:param verbose: Write the progress for the single run
:return: A 2d array of simulated soil moisture values per depth
'''
# Set the parameters
par = par or self.make_parameters()
self.model.set_parameters(par)
if verbose:
print('Parameters:')
for k, v in par._asdict().items():
print(' {} = {:0.4g}'.format(k,v))
eval_layers = self.get_eval_layers(self.eval_depth)
# Prepare result array with nan's
result = np.nan * np.ones(shape=(len(self.evaluation_data), len(eval_layers)))
# Save the starting conditions
result[0] = [l.theta for l in eval_layers]
reporter = _ProgressReporter(self.datastart, self.dataend, max_runtime=15 * 60, verbose=verbose)
try:
for i, t in enumerate(self.model.solver.run(self.datastart, self.dataend, timedelta(hours=1))):
self.model.set_boundary(t)
# Get result
result[i+1] = [l.theta for l in eval_layers]
# Raises if the runtime is too long
reporter(t)
except RuntimeError as e:
sys.stderr.write(str(par))
sys.stderr.write(str(e))
return result
if __name__ == '__main__':
print(spotpy.__file__, spotpy.__version__)
if len(sys.argv) > 1:
runs = int(sys.argv[1])
else:
runs = 1
model = Cmf1d_Model()
print(spotpy.describe.setup(model))
if runs > 1:
parallel = 'mpi' if 'OMPI_COMM_WORLD_SIZE' in os.environ else 'seq'
sampler = spotpy.algorithms.mc(model,
dbformat='csv',
dbname='cmf1d_mc',
parallel=parallel,
save_sim=False)
sampler.sample(runs)
else:
from spotpy.gui.mpl import GUI
gui = GUI(model)
gui.show()
