PrecipitationDistribution: Generate random sequence of precipitation events#

class PrecipitationDistribution(*args, **kwds)[source]#

Bases: Component

Generate precipitation events.

This component can generate a random storm duration, interstorm duration, precipitation intensity or storm depth from a Poisson distribution when given a mean value.

Examples

>>> from landlab.components.uniform_precip import PrecipitationDistribution
>>> import numpy as np
>>> np.random.seed(np.arange(10))

To use hard-coded values for mean storm, mean interstorm, mean depth, model run time and delta t… Say we use 1.5 for mean storm, 15 for mean interstorm, 0.5 for mean depth, 100 for model run time and 1 for delta t…

>>> precip = PrecipitationDistribution(
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=100.0,
...     delta_t=1.0,
... )
>>> for dt, rate in precip.yield_storm_interstorm_duration_intensity():
...     pass  # and so on
...

Alternatively, we can pass a grid to the component, and call yield_storms() to generate storm-interstorm float pairs while the intensity data is stored in the grid scalar field ‘rainfall__flux’:

>>> from landlab import RasterModelGrid
>>> mg = RasterModelGrid((4, 5))
>>> precip = PrecipitationDistribution(
...     mg,
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=46.0,
... )
>>> storm_dts = []
>>> interstorm_dts = []
>>> intensities = []
>>> precip.seed_generator(seedval=1)
>>> for storm_dt, interstorm_dt in precip.yield_storms():
...     storm_dts.append(storm_dt)
...     interstorm_dts.append(interstorm_dt)
...     intensities.append(mg.at_grid["rainfall__flux"])
...
>>> len(storm_dts) == 4  # 4 storms in the simulation
True
>>> len(interstorm_dts) == len(storm_dts)
True
>>> rf_intensities_to_test = np.array(
...     [
...         0.8138257984406472,
...         0.15929112025199238,
...         0.17254519305000884,
...         0.09817611240558813,
...     ]
... )
>>> np.allclose(intensities, rf_intensities_to_test)
True
>>> np.isclose(sum(storm_dts) + sum(interstorm_dts), 46.0)  # test total_t
True
>>> np.isclose(interstorm_dts[-1], 0.0)  # sequence truncated as necessary
True

We can also turn the generator straight into a list, like this:

>>> precip.seed_generator(seedval=1)
>>> steps = [(dt + istorm_dt) for (dt, istorm_dt) in precip.yield_storms()]
>>> np.isclose(sum(steps), 46.0)
True

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Eagleson, P. (1978). Climate, soil, and vegetation: 2. The distribution of annual precipitation derived from observed storm sequences. Water Resources Research 14(5), 713-721. https://dx.doi.org/10.1029/wr014i005p00713

Create the storm generator.

Parameters:
  • grid (ModelGrid) – A Landlab grid (optional). If provided, storm intensities will be stored as a grid scalar field as the component simulates storms.

  • mean_storm_duration (float) – Average duration of a precipitation event.

  • mean_interstorm_duration (float) – Average duration between precipitation events.

  • mean_storm_depth (float) – Average depth of precipitation events.

  • total_t (float, optional) – If generating a time series, the total amount of time.

  • delta_t (float or None, optional) – If you want to break up storms into determined subsections using yield_storm_interstorm_duration_intensity, a delta_t is needed.

  • random_seed (int or float, optional) – Seed value for random-number generator.

__init__(grid=None, mean_storm_duration=1.0, mean_interstorm_duration=1.0, mean_storm_depth=1.0, total_t=0.0, delta_t=None, random_seed=0)[source]#

Create the storm generator.

Parameters:
  • grid (ModelGrid) – A Landlab grid (optional). If provided, storm intensities will be stored as a grid scalar field as the component simulates storms.

  • mean_storm_duration (float) – Average duration of a precipitation event.

  • mean_interstorm_duration (float) – Average duration between precipitation events.

  • mean_storm_depth (float) – Average depth of precipitation events.

  • total_t (float, optional) – If generating a time series, the total amount of time.

  • delta_t (float or None, optional) – If you want to break up storms into determined subsections using yield_storm_interstorm_duration_intensity, a delta_t is needed.

  • random_seed (int or float, optional) – Seed value for random-number generator.

property elapsed_time#

Get the elapsed time recorded by the module.

This will be particularly useful in the midst of a yield loop.

generate_from_stretched_exponential(scale, shape)[source]#

Generate and return a random variable from a stretched exponential distribution with given scale and shape.

Examples

>>> np.random.seed(0)
>>> np.round(np.random.rand(3), 6)  # these are our 3 rand #s to test
array([ 0.548814,  0.715189,  0.602763])
>>> from landlab.components import PrecipitationDistribution
>>> pd = PrecipitationDistribution(
...     mean_storm_duration=1.0,
...     mean_interstorm_duration=1.0,
...     mean_storm_depth=1.0,
... )
>>> np.random.seed(0)  # re-set seed so we get the same 3 #s
>>> np.round(1000 * pd.generate_from_stretched_exponential(2.0, 0.5))
720.0
>>> np.round(1000 * pd.generate_from_stretched_exponential(2.0, 0.5))
225.0
>>> np.round(1000 * pd.generate_from_stretched_exponential(2.0, 0.5))
513.0
get_interstorm_event_duration()[source]#

Generate interstorm events.

This method takes one argument, the mean_interstorm_duration parameter. (In Eagleson (1978), this parameter was called Tb.)

This method is modeled identically to get_precipitation_event_duration()

This method finds a random value for interstorm_duration based on the poisson distribution about the mean. This is accomplished using the expovariate function from the “random” standard library.

Returns:

The interstorm duration.

Return type:

float

get_precipitation_event_duration()[source]#

This method is the storm generator.

This method has one argument: the mean_storm_duration parameter. (In Eagleson (1978), this parameter was called Tr.)

It finds a random storm_duration value based on the poisson distribution about the mean. This is accomplished using the expovariate function from the “random” standard library.

Returns:

The storm duration.

Return type:

float

get_storm_depth()[source]#

Generate storm depth.

Storm depth is used to generate a realistic intensity for different storm events.

(In Eagleson (1978) this parameter was called “h”)

This method requires storm_duration, mean_storm_duration and the mean_storm_depth. Storm_duration is generated through the initialize() or update() method.

Numpy has a random number generator to get values from a given Gamma distribution. It takes two arguments, alpha (or the shape parameter), which is the generated over the mean event and beta (or the scale parameter), which is the mean value These are all arguments in the function, which returns storm depth.

Returns:

The storm depth.

Return type:

float

get_storm_intensity()[source]#

Get the storm intensity.

This method draws storm intensity out of the storm depth generated by get_storm_depth.

This method requires the storm_depth and storm_duration and is the same as the parameter (“i”) in Eagleson (1978), but instead of being drawn from Poission, this is drawn from the Gamma distribution of (h), as \(h = i * Tr\).

Returns:

The storm intensity.

Return type:

float

get_storm_time_series()[source]#

Get a time series of storms.

This method creates a time series of storms based on storm_duration, and interstorm_duration. From these values it will calculate a complete time series.

The storm_time_series returned by this method is made up of sublists, each comprising of three sub-parts (e.g. [[x,y,z], [a,b,c]]) where x and a are the beginning times of a precipitation event, y and b are the ending times of the precipitation event and z and c represent the average intensity (mm/hr) of the storm lasting from x to y and a to be, respectively.

Even if a grid was passed to the component at instantiation, calling this method does not update the grid fields.

Returns:

containing several sub-arrays of events [start, finish, intensity]

Return type:

array

property intensity#

Get the intensity of the most recent storm simulated.

property interstorm_duration#

Interstorm duration.

[T]

seed_generator(seedval=0)[source]#

Seed the random-number generator.

The examples illustrate:

  1. That we can get the same sequence again by re-seeding with the same value (the default is zero)

  2. When we use a value other than the default, we get a different sequence

Examples

>>> precip = PrecipitationDistribution(
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=100.0,
...     delta_t=1.0,
... )
>>> round(precip.storm_duration, 2)
2.79
>>> round(precip.interstorm_duration, 2)
21.28
>>> round(precip.storm_depth, 2)
2.45
>>> round(precip.intensity, 2)
0.88
>>> precip.seed_generator()  # re-seed and get same sequence again
>>> round(precip.get_precipitation_event_duration(), 2)
2.79
>>> round(precip.get_interstorm_event_duration(), 2)
21.28
>>> round(precip.get_storm_depth(), 2)
2.45
>>> round(precip.get_storm_intensity(), 2)
0.88
>>> precip = PrecipitationDistribution(
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=100.0,
...     delta_t=1.0,
...     random_seed=1,
... )
>>> round(precip.storm_duration, 2)  # diff't vals with diff't seed
0.22
>>> round(precip.interstorm_duration, 2)
28.2
>>> round(precip.storm_depth, 4)
0.0012
>>> round(precip.intensity, 4)
0.0054
property storm_depth#

Depth of water in the storm.

[L]

property storm_duration#

Duration of storm.

[T]

update()[source]#

Update the storm values.

If new values for storm duration, interstorm duration, storm depth and intensity are needed, this method can be used to update those values one time.

Examples

>>> from landlab.components import PrecipitationDistribution
>>> precip = PrecipitationDistribution(
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=100.0,
...     delta_t=1,
... )

Additionally, if we wanted to update several times, a loop could be utilized to accomplish this. Say we want 5 storm_durations; this pseudo-code represents a way to accomplish this…

>>> storm_duration_list = []
>>> i = 0
>>> while i < 4:
...     storm_duration_list.append(precip.storm_duration)
...     precip.update()
...     i += 1
...

Note though that alternatively we could also do this, avoiding the method entirely…

>>> # ^^this lets you "manually" get the next item from the iterator
>>> precip = PrecipitationDistribution(
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=46.0,
... )
>>> storm_duration_list = []
>>> for i in range(5):
...     storm_duration_list.append(
...         next(precip.yield_storm_interstorm_duration_intensity())[0]
...     )
...

Notice that doing this will not automatically stop the iteration, however - it will continue ad infinitum.

yield_storm_interstorm_duration_intensity(subdivide_interstorms=False)[source]#

Iterator for a time series of storms interspersed with interstorms.

This method is intended to be equivalent to get_storm_time_series, but instead offers a generator functionality. This will be useful in cases where the whole sequence of storms and interstorms doesn’t need to be stored, where we can save memory this way.

The method keeps track of the delta_t such that if a storm needs to be generated longer than this supplied model timestep, the generator will return the storm in “chunks”, until there is no more storm duration. e.g., storm of intensity 1. is 4.5 long, the delta_t is 2., the generator yields (2.,1.) -> (2.,1.) -> (0.5,1.) -> …

If delta_t is None or not supplied, no subdivision occurs.

Once a storm has been generated, this method will follow it with the next interstorm, yielded as (interstorm_duration, 0.). Note that the interstorm will NOT be subdivided according to delta_t unless you set the flag subdivide_interstorms to True.

The method will keep yielding until it reaches the RUN_TIME, where it will terminate.

Yields:

tuple of float – (interval_duration, rainfall_rate_in_interval)

Notes

One recommended procedure is to instantiate the generator, then call instance.next() (in Python 2) or next(instance) (in Python 3) repeatedly to get the sequence.

yield_storms()[source]#

Iterator for a time series of storm-interstorm pairs.

This method is very similar to this component’s other generator, yield_storm_interstorm_duration_intensity(), but the way it yields is slightly different. Instead of yielding (interval_duration, rf_rate), with interstorms represented as intervals with rf_rate = 0, it yields:

(storm_duration, interstorm_duration)

When each tuple pair is yielded, the grid scalar field ‘rainfall__flux’ is updated with the rainfall rate occuring during storm_duration.

This generator method is designed for direct equivalence with the spatially resolved generators found elsewhere in Landlab.

This generator will be useful in cases where the whole sequence of storms and interstorms doesn’t need to be stored, where we can save memory this way.

This method does not attempt to subdivide timesteps. If you want that, provide a delta_t at instantiation, and use yield_storm_interstorm_duration_intensity(subdivide_interstorms=True).

The method will keep yielding until it reaches the RUN_TIME, where it will terminate. If it terminates during a storm, the final tuple will be (truncated_storm_duration, 0.). Otherwise it will be (storm_duration, truncated_interstorm_duration.)

Yields:

tuple of float – (storm_duration, interstorm_duration)

Notes

One recommended procedure is to instantiate the generator, then call next(instance) (in Python 3) repeatedly to get the sequence (See Examples, below).

Examples

>>> from landlab import RasterModelGrid
>>> mg = RasterModelGrid((4, 5))
>>> precip = PrecipitationDistribution(
...     mg,
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=46.0,
... )
>>> storm_dts = []
>>> interstorm_dts = []
>>> intensities = []
>>> precip.seed_generator(seedval=1)
>>> for storm_dt, interstorm_dt in precip.yield_storms():
...     storm_dts.append(storm_dt)
...     interstorm_dts.append(interstorm_dt)
...     intensities.append(mg.at_grid["rainfall__flux"])
...
>>> len(storm_dts) == 4  # 4 storms in the simulation
True
>>> len(interstorm_dts) == len(storm_dts)
True
>>> rf_intensities_to_test = np.array(
...     [
...         0.8138257984406472,
...         0.15929112025199238,
...         0.17254519305000884,
...         0.09817611240558813,
...     ]
... )
>>> np.allclose(intensities, rf_intensities_to_test)
True
>>> np.isclose(sum(storm_dts) + sum(interstorm_dts), 46.0)  # total_t
True
>>> np.isclose(interstorm_dts[-1], 0.0)  # sequence truncated
True

An alternative way to use the generator might be:

>>> # ^^this lets you "manually" get the next item from the iterator
>>> flux = mg.at_grid.pop("rainfall__flux")  # remove the existing field
>>> precip = PrecipitationDistribution(
...     mg,
...     mean_storm_duration=1.5,
...     mean_interstorm_duration=15.0,
...     mean_storm_depth=0.5,
...     total_t=46.0,
... )
>>> precip.seed_generator(seedval=1)
>>> mystorm_generator = precip.yield_storms()
>>> my_list_of_storms = []
>>> for i in range(4):
...     my_list_of_storms.append(next(mystorm_generator))
...

Note that an exception is thrown if you go too far:

>>> my_list_of_storms.append(next(mystorm_generator))  # the 5th iter
Traceback (most recent call last):
  ...
StopIteration

Also note that the generator won’t terminate if you try this without first instantiating the generator:

>>> allmytimes = []
>>> for i in range(20):  # this will run just fine
...     allmytimes.append(next(precip.yield_storms()))
...
>>> total_t = sum([sum(storm) for storm in allmytimes])
>>> total_t > 46.0
True