# coding: utf8
# !/usr/env/python
"""Base class for common functions of terrainbento stochastic erosion
models."""
import os
import textwrap
import numpy as np
import scipy.stats as stats
from landlab.components import PrecipitationDistribution
from terrainbento.base_class import ErosionModel
_STRING_LENGTH = 80
[docs]class StochasticErosionModel(ErosionModel):
"""Base class for stochastic-precipitation terrainbento models.
A **StochasticErosionModel** inherits from
:py:class:`ErosionModel` and provides functionality needed by all
stochastic-precipitation models.
This is a base class that handles processes related to the generation of
preciptiation events.
Two primary options are avaliable for the stochastic erosion models. When
``opt_stochastic_duration=True`` the model will use the
`PrecipitationDistribution <https://landlab.readthedocs.io/en/master/reference/components/uniform_precip.html>`_
Landlab component to generate a random storm duration, interstorm duration,
and precipitation intensity or storm depth from an exponential
distribution. When this option is selected, the following parameters are
used:
- mean_storm_duration
- mean_interstorm_duration
- mean_storm_depth
When ``opt_stochastic_duration==False`` the model will have uniform
timesteps but generate rainfall from a stretched exponential distribution.
The duration indicated by the parameter ``step`` will first be split into a
series of sub-timesteps based on the parameter
``number_of_sub_time_steps``, and then each of these sub-timesteps will
experience a duration of rain and no-rain based on the value of
``rainfall_intermittency_factor``. The duration of rain and no-rain will
not change, but the intensity of rain will vary based on a stretched
exponential distribution described by the shape factor
``rainfall__shape_factor`` and with a scale factor calculated so that the
mean of the distribution has the value given by ``rainfall__mean_rate``.
The following parameters are used:
- rainfall__shape_factor
- number_of_sub_time_steps
- rainfall_intermittency_factor
- rainfall__mean_rate
The hydrology uses calculation of drainage area using the user-specified
routing method. It then performs one of two options, depending on the
user's choice of ``opt_stochastic_duration`` (True or False).
If the user requests stochastic duration, the model iterates through a sequence
of storm and interstorm periods. Storm depth is drawn at random from a gamma
distribution, and storm duration from an exponential distribution; storm
intensity is then depth divided by duration. This sequencing is implemented by
overriding the run_for method.
If the user does not request stochastic duration (indicated by setting
``opt_stochastic_duration`` to ``False``), then the default
(**erosion_model** base class) **run_for** method is used. Whenever
**run_one_step** is called, storm intensity is generated at random from an
exponential distribution with mean given by the parameter
``rainfall__mean_rate``. The stream power component is run for only a
fraction of the time step duration step, as specified by the parameter
``rainfall_intermittency_factor``. For example, if ``step`` is 10 years and
the intermittency factor is 0.25, then the stream power component is run
for only 2.5 years.
In either case, given a storm precipitation intensity :math:`P`, the runoff
production rate :math:`R` [L/T] is calculated using:
.. math::
R = P - I (1 - \exp ( -P / I ))
where :math:`I` is the soil infiltration capacity. At the sub-grid scale, soil
infiltration capacity is assumed to have an exponential distribution of which
:math:`I` is the mean. Hence, there are always some spots within any given grid cell
that will generate runoff. This approach yields a smooth transition from
near-zero runoff (when :math:`I>>P`) to :math:`R \\approx P`
(when :math:`P>>I`), without a "hard threshold."
The following at-node fields must be specified in the grid:
- ``topographic__elevation``
"""
_required_fields = ["topographic__elevation"]
[docs] def __init__(
self,
clock,
grid,
random_seed=0,
record_rain=False,
opt_stochastic_duration=False,
mean_storm_duration=1,
mean_interstorm_duration=1,
mean_storm_depth=1,
rainfall__shape_factor=1,
number_of_sub_time_steps=1,
rainfall_intermittency_factor=1,
rainfall__mean_rate=1,
storm_sequence_filename="storm_sequence.txt",
frequency_filename="exceedance_summary.txt",
**kwargs
):
"""
Parameters
----------
clock : terrainbento Clock instance
grid : landlab model grid instance
The grid must have all required fields.
random_seed, int, optional
Random seed. Default is 0.
opt_stochastic_duration : bool, optional
Flag indicating if timestep is stochastic or constant. Default is
False.
mean_storm_duration : float, optional
Average duration of a precipitation event. Default is 1.
mean_interstorm_duration : float, optional
Average duration between precipitation events. Default is 1.
mean_storm_depth : float, optional
Average depth of precipitation events. Default is 1.
number_of_sub_time_steps : int, optional
Number of sub-timesteps. Default is 1.
rainfall_intermittency_factor : float, optional
Value between zero and one that indicates the proportion of time
rain occurs. A value of 0 means it never rains and a value of 1
means that rain never ceases. Default is 1.
rainfall__mean_rate : float, optional
Mean of the precipitation distribution. Default is 1.
rainfall__shape_factor : float, optional
Shape factor of the precipitation distribution. Default is 1.
record_rain : boolean
Flag to indicate if a sequence of storms should be saved. Default
is False.
storm_sequence_filename : str
Storm sequence filename. Default is "storm_sequence.txt"
frequency_filename : str
Filename for precipitation exceedance frequency summary. Default
value is "exceedance_summary.txt"
**kwargs :
Keyword arguments to pass to
:py:class:`ErosionModel`
Returns
-------
StochasticErosionModel : object
Examples
--------
This model is a base class and is not designed to be run on its own. We
recommend that you look at the terrainbento tutorials for examples of
usage.
"""
# Call StochasticErosionModel init
super().__init__(clock, grid, **kwargs)
# ensure Precipitator and RunoffGenerator are vanilla
self._ensure_precip_runoff_are_vanilla()
self.opt_stochastic_duration = opt_stochastic_duration
# verify that opt_stochastic_duration and PrecipChanger are consistent
if self.opt_stochastic_duration and (
"PrecipChanger" in self.boundary_handlers
):
msg = (
"terrainbento StochasticErosionModel: setting "
"opt_stochastic_duration=True and using the PrecipChanger "
"boundary condition handler are not compatible."
)
raise ValueError(msg)
self.seed = int(random_seed)
self.random_seed = random_seed
self.frequency_filename = frequency_filename
self.storm_sequence_filename = storm_sequence_filename
self.mean_storm_duration = mean_storm_duration
self.mean_interstorm_duration = mean_interstorm_duration
self.mean_storm_depth = mean_storm_depth
self.shape_factor = rainfall__shape_factor
self.number_of_sub_time_steps = number_of_sub_time_steps
self.rainfall_intermittency_factor = rainfall_intermittency_factor
self.rainfall__mean_rate = rainfall__mean_rate
# initialize record for storms. Depending on how this model is run
# (stochastic time, number_time_steps>1, more manually) the step may
# change. Thus, rather than writing routines to reconstruct the time
# series of precipitation from the step could change based on users use,
# we"ll record this with the model run instead of re-running.
# make this the non-default option.
# Second, test that
if record_rain:
self.record_rain = True
self.rain_record = {
"event_start_time": [],
"event_duration": [],
"rainfall_rate": [],
"runoff_rate": [],
}
else:
self.record_rain = False
self.rain_record = None
[docs] def calc_runoff_and_discharge(self):
"""Calculate runoff rate and discharge; return runoff."""
if self.rain_rate > 0.0 and self.infilt > 0.0:
runoff = self.rain_rate - (
self.infilt * (1.0 - np.exp(-self.rain_rate / self.infilt))
)
if runoff <= 0:
runoff = 0 # pragma: no cover
else:
runoff = self.rain_rate
self.grid.at_node["surface_water__discharge"][:] = (
runoff * self.grid.at_node["drainage_area"]
)
return runoff
[docs] def run_for_stochastic(self, step, runtime):
"""**Run_for** with stochastic duration.
Run model without interruption for a specified time period, using
random storm/interstorm sequence.
**run_for_stochastic** runs the model for the duration ``runtime`` with
model time steps given by the PrecipitationDistribution component.
Model run steps will not exceed the duration given by ``step``.
Parameters
----------
step : float
Model run timestep,
runtime : float
Total duration for which to run model.
"""
self.rain_generator._delta_t = step
self.rain_generator._run_time = runtime
for (
tr,
p,
) in self.rain_generator.yield_storm_interstorm_duration_intensity():
self.rain_rate = p
self.run_one_step(tr)
[docs] def instantiate_rain_generator(self):
"""Instantiate component used to generate storm sequence."""
# Handle option for duration.
if self.opt_stochastic_duration:
self.rain_generator = PrecipitationDistribution(
mean_storm_duration=self.mean_storm_duration,
mean_interstorm_duration=self.mean_interstorm_duration,
mean_storm_depth=self.mean_storm_depth,
total_t=self.clock.stop,
delta_t=self.clock.step,
random_seed=self.seed,
)
self.run_for = self.run_for_stochastic # override base method
else:
from scipy.special import gamma
self.rain_generator = PrecipitationDistribution(
mean_storm_duration=1.0,
mean_interstorm_duration=1.0,
mean_storm_depth=1.0,
random_seed=self.seed,
)
self.scale_factor = self.rainfall__mean_rate / gamma(
1.0 + (1.0 / self.shape_factor)
)
if (
isinstance(self.number_of_sub_time_steps, (int, np.integer))
is False
):
raise ValueError(
("number_of_sub_time_steps must be of type integer.")
)
self.n_sub_steps = self.number_of_sub_time_steps
[docs] def reset_random_seed(self):
"""Reset the random number generation sequence."""
self.rain_generator.seed_generator(seedval=self.seed)
def _pre_water_erosion_steps(self):
"""Convenience function for pre-water erosion steps.
If a model needs to do anything before each erosion step is run,
e.g. recalculate a threshold value, that model should overwrite
this function.
"""
pass
[docs] def handle_water_erosion(self, step):
"""Handle water erosion for stochastic models.
If we are running stochastic duration, then self.rain_rate will
have been calculated already. It might be zero, in which case we
are between storms, so we don't do water erosion.
If we're NOT doing stochastic duration, then we'll run water
erosion for one or more sub-time steps, each with its own
randomly drawn precipitation intensity.
This routine assumes that a model-specific method:
**calc_runoff_and_discharge()**
will have been defined. Additionally a model eroder must also have been
defined.
For example, BasicStVs calculated runoff and discharge in a different
way than the other models.
If the model has a function **update_threshold_field**, this
function will test for it and run it. This is presently done in
BasicDdSt.
Parameters
----------
step : float
Model run timestep.
"""
# (if we're varying precipitation parameters through time, update them)
if "PrecipChanger" in self.boundary_handlers:
(
self.rainfall_intermittency_factor,
self.rainfall__mean_rate,
) = self.boundary_handlers[
"PrecipChanger"
].get_current_precip_params()
if self.opt_stochastic_duration and self.rain_rate > 0.0:
self._pre_water_erosion_steps()
runoff = self.calc_runoff_and_discharge()
self.eroder.run_one_step(step)
if self.record_rain:
# save record into the rain record
self.record_rain_event(
self.model_time, step, self.rain_rate, runoff
)
elif self.opt_stochastic_duration and self.rain_rate <= 0.0:
# calculate and record the time with no rain:
if self.record_rain:
self.record_rain_event(self.model_time, step, 0, 0)
elif not self.opt_stochastic_duration:
dt_water = (step * self.rainfall_intermittency_factor) / float(
self.n_sub_steps
)
for i in range(self.n_sub_steps):
self.rain_rate = self.rain_generator.generate_from_stretched_exponential(
self.scale_factor, self.shape_factor
)
self._pre_water_erosion_steps()
runoff = self.calc_runoff_and_discharge()
self.eroder.run_one_step(dt_water)
# save record into the rain record
if self.record_rain:
event_start_time = self.model_time + (i * dt_water)
self.record_rain_event(
event_start_time, dt_water, self.rain_rate, runoff
)
# once all the rain time_steps are complete,
# calculate and record the time with no rain:
if self.record_rain:
# calculate dry time
dt_dry = step * (1 - self.rainfall_intermittency_factor)
# if dry time is greater than zero, record.
if dt_dry > 0:
event_start_time = self.model_time + (
self.n_sub_steps * dt_water
)
self.record_rain_event(event_start_time, dt_dry, 0.0, 0.0)
[docs] def finalize(self):
"""Finalize stochastic erosion models.
The finalization step of stochastic erosion models in
terrainbento results in writing out the storm sequence file and
the precipitation exceedence statistics summary if
``record_rain`` was set to ``True``.
"""
# if rain was recorded, write it out.
if self.record_rain:
self.write_storm_sequence_to_file(
filename=self.storm_sequence_filename
)
if self.opt_stochastic_duration is False:
# if opt_stochastic_duration is False, calculate exceedance
# frequencies and write out.
try:
self.write_exceedance_frequency_file(
filename=self.frequency_filename
)
except IndexError:
msg = (
"terrainbento stochastic model: the rain record was "
"too short to calculate exceedance frequency statistics."
)
os.remove(self.frequency_filename)
raise RuntimeError(msg)
[docs] def record_rain_event(
self, event_start_time, event_duration, rainfall_rate, runoff_rate
):
"""Record rain events.
Create a record of event start time, event duration, rainfall rate, and
runoff rate.
Parameters
----------
event_start_time : float
event_duration : float
rainfall_rate : float
runoff_rate : float
"""
self.rain_record["event_start_time"].append(event_start_time)
self.rain_record["event_duration"].append(event_duration)
self.rain_record["rainfall_rate"].append(rainfall_rate)
self.rain_record["runoff_rate"].append(runoff_rate)
[docs] def write_storm_sequence_to_file(self, filename="storm_sequence.txt"):
"""Write event duration and intensity to a formatted text file.
Parameters
----------
filename : str
Default value is "storm_sequence.txt"
"""
# Open a file for writing
if self.record_rain is False:
raise ValueError(
"Rain was not recorded when the model run. To "
"record rain, set the parameter 'record_rain'"
"to True."
)
with open(filename, "w") as stormfile:
stormfile.write(
"event_start_time"
+ ","
+ "event_duration"
+ ","
+ "rainfall_rate"
+ ","
+ "runoff_rate"
+ "\n"
)
n_events = len(self.rain_record["event_start_time"])
for i in range(n_events):
stormfile.write(
str(
np.around(
self.rain_record["event_start_time"][i], decimals=5
)
)
+ ","
+ str(
np.around(
self.rain_record["event_duration"][i], decimals=5
)
)
+ ","
+ str(
np.around(
self.rain_record["rainfall_rate"][i], decimals=5
)
)
+ ","
+ str(
np.around(
self.rain_record["runoff_rate"][i], decimals=5
)
)
+ "\n"
)
[docs] def write_exceedance_frequency_file(
self, filename="exceedance_summary.txt"
):
"""Write summary of rainfall exceedance statistics to file.
Parameters
----------
filename : str
Default value is "exceedance_summary.txt"
"""
if self.record_rain is False:
raise ValueError(
"Rain was not recorded when the model run. To "
"record rain, set the parameter 'record_rain'"
"to True."
)
# calculate the number of wet days per year.
number_of_days_per_year = 365
nwet = int(
np.ceil(
self.rainfall_intermittency_factor * number_of_days_per_year
)
)
if nwet == 0:
raise ValueError(
"No rain fell, which makes calculating exceedance "
"frequencies problematic. We recommend that you "
"check the valude of rainfall_intermittency_factor."
)
with open(filename, "w") as exceedance_file:
# ndry = int(number_of_days_per_year - nwet)
# Write some basic information about the distribution to the file.
exceedance_file.write("Section 1: Distribution Description\n")
exceedance_file.write(
"Scale Factor: " + str(self.scale_factor) + "\n"
)
exceedance_file.write(
"Shape Factor: " + str(self.shape_factor) + "\n"
)
exceedance_file.write(
(
"Intermittency Factor: "
+ str(self.rainfall_intermittency_factor)
+ "\n"
)
)
exceedance_file.write(
("Number of wet days per year: " + str(nwet) + "\n\n")
)
message_text = (
"The scale factor that describes this distribution is "
+ "calculated based on a provided value for the mean wet day rainfall."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
exceedance_file.write(
(
"This provided value was:\n"
+ str(self.rainfall__mean_rate)
+ "\n"
)
)
# calculate the predictions for 10, 25, and 100 year event based on
# the analytical form of the exceedance function.
event_intervals = np.array([10.0, 25, 100.0])
# calculate the probability of each event based on the number of years
# and the number of wet days per year.
daily_distribution_exceedance_probabilities = 1.0 / (
nwet * event_intervals
)
# exceedance probability is given as
# Probability of daily rainfall of p exceeding a value of po is given as:
#
# P(p>po) = e^(-(po/P)^c)
# P = scale
# c = shape
#
# this can be re-arranged to
#
# po = P * (- ln (P(p>po))) ^ (1 / c)
expected_rainfall = self.scale_factor * (
-1.0 * np.log(daily_distribution_exceedance_probabilities)
) ** (1.0 / self.shape_factor)
exceedance_file.write("\n\nSection 2: Theoretical Predictions\n")
message_text = (
"Based on the analytical form of the wet day rainfall "
+ "distribution, we can calculate theoretical predictions "
+ "of the daily rainfall amounts associated with N-year events."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
for i in range(len(daily_distribution_exceedance_probabilities)):
exceedance_file.write(
(
"Expected value for the wet day total of the "
+ str(event_intervals[i])
+ " year event is: "
+ str(np.round(expected_rainfall[i], decimals=3))
+ "\n"
)
)
# get rainfall record and filter out time without any rain
all_precipitation = np.array(self.rain_record["rainfall_rate"])
rainy_day_inds = np.where(all_precipitation > 0)
wet_day_totals = all_precipitation[rainy_day_inds]
num_days = len(wet_day_totals)
# construct the distribution of yearly maxima.
# here an effective year is represented by the number of draws implied
# by the intermittency factor
# first calculate the number of effective years.
num_effective_years = int(np.floor(wet_day_totals.size / nwet))
# write out the calculated event only if the duration
exceedance_file.write("\n\n")
message_text = (
"Section 3: Predicted 95% confidence bounds on the "
+ "exceedance values based on number of samples drawn."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
message_text = (
"The ability to empirically estimate the rainfall "
+ "associated with an N-year event depends on the "
+ "probability of that event occurring and the number of "
+ "draws from the probability distribution. The ability "
+ "to estimate increases with an increasing number of samples "
+ "and decreases with decreasing probability of event "
+ "occurrence."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
message_text = (
"Exceedance values calculated from "
+ str(len(wet_day_totals))
+ " draws from the daily-rainfall probability distribution. "
+ "This corresponds to "
+ str(num_effective_years)
+ " effective years."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
# For a general probability distribution, f, with a continuous not zero
# quantile function at F-1(p), the order statistic associated with the
# p percentile given n draws from the distribution is given as:
# X[np] ~ AN ( F-1(p), (p * (p - 1 ))/ (n * [f (F-1 (p)) ]**2))
# where AN is the asymptotic normal. The value for the variance is more
# intuitive once you consider that [f (F-1 (p)) ] is the probability
# that an event of percentile p will occur. Thus the variance increases
# non-linearly with decreasing event probability and decreases linearly
# with increaseing observations.
# we"ve already calculated F-1(p) for our events, and it is represented
# by the variable expected_rainfall
daily_distribution_event_percentile = (
1.0 - daily_distribution_exceedance_probabilities
)
event_probability = (
(self.shape_factor / self.scale_factor)
* (
(expected_rainfall / self.scale_factor)
** (self.shape_factor - 1.0)
)
* (
np.exp(
-1.0
* (expected_rainfall / self.scale_factor)
** self.shape_factor
)
)
)
event_variance = (
daily_distribution_event_percentile
* (1.0 - daily_distribution_event_percentile)
) / (num_days * (event_probability ** 2))
event_std = event_variance ** 0.5
t_statistic = stats.t.ppf(
0.975, num_effective_years, loc=0, scale=1
)
exceedance_file.write("\n")
message_text = (
"For the given number of samples, the 95% "
+ "confidence bounds for the following event "
+ "return intervals are as follows: "
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
for i in range(len(event_intervals)):
min_expected_val = (
expected_rainfall[i] - t_statistic * event_std[i]
)
max_expected_val = (
expected_rainfall[i] + t_statistic * event_std[i]
)
exceedance_file.write(
(
"Expected range for the wet day total of the "
+ str(event_intervals[i])
+ " year event is: ("
+ str(np.round(min_expected_val, decimals=3))
+ ", "
+ str(np.round(max_expected_val, decimals=3))
+ ")\n"
)
)
# next, calculate the emperical exceedance values, if a sufficient record
# exists.
# inititialize a container for the maximum yearly precipitation.
maximum_yearly_precipitation = np.nan * np.zeros(
(num_effective_years)
)
for yi in range(num_effective_years):
# identify the starting and ending index coorisponding to the
# year
starting_index = yi * nwet
ending_index = starting_index + nwet
# select the years portion of the wet_day_totals
selected_wet_day_totals = wet_day_totals[
starting_index:ending_index
]
# record the yearly maximum precipitation
maximum_yearly_precipitation[
yi
] = selected_wet_day_totals.max()
# calculate the distribution percentiles associated with each interval
event_percentiles = (1.0 - (1.0 / event_intervals)) * 100.0
# calculated the event magnitudes associated with the percentiles.
event_magnitudes = np.percentile(
maximum_yearly_precipitation, event_percentiles
)
# write out the calculated event only if the duration
exceedance_file.write("\n\nSection 4: Empirical Values\n")
message_text = (
"These empirical values should be interpreted in the "
+ "context of the expected ranges printed in Section 3. "
+ "If the expected range is large, consider using a longer "
+ "record of rainfall. The empirical values should fall "
+ "within the expected range at a 95% confidence level."
)
exceedance_file.write(
"\n".join(textwrap.wrap(message_text, _STRING_LENGTH))
)
exceedance_file.write("\n")
for i in range(len(event_percentiles)):
exceedance_file.write(
(
"Estimated value for the wet day total of the "
+ str(np.round(event_intervals[i], decimals=3))
+ " year event is: "
+ str(np.round(event_magnitudes[i], decimals=3))
+ "\n"
)
)
[docs]def main(): # pragma: no cover
"""Executes model."""
import sys
try:
infile = sys.argv[1]
except IndexError:
print("Must include input file name on command line")
sys.exit(1)
sm = StochasticErosionModel.from_file(infile)
sm.run()
if __name__ == "__main__":
main()
