import numpy as np
from corehydropy import (
Distribution,
TimeSeries,
autocorrelation,
fit_mle,
model_point_process,
model_univariate,
)30. Time series
The Numerics TimeSeries container, which had no R or Python binding until this release. It is the thing a flood-frequency analysis actually starts from: a gauge record with gaps in it, out of which you have to get an annual maximum series, or a set of independent peaks, before any of the fitting machinery is any use.
This page walks one thirty-year daily record from raw gauge output to two fitted flood-frequency models, and then generates synthetic records from it.
What you’ll learn
- How to build a
TimeSeries, and why the interval it carries is not just a label. - How to repair a record: filling dates that are missing entirely, and interpolating short runs of missing values – and what the library does when a run is too long.
- How to reduce a record to annual maxima over a water year, and what happens to a year that contains a missing value (it depends on the block function, and the reason is worth knowing).
- How to extract independent peaks over a threshold, and what the independence criterion is measured in.
- How to decompose a series into trend, seasonal and residual components.
- How to generate synthetic records with two resamplers, and why a seeded run gives identical numbers in Python and R.
Setup
A thirty-year daily record
The record below is illustrative rather than observed, but it is built to be exactly reproducible: a twelve-value monthly base flow, modulated by a small integer linear congruential generator written out in full, so Python and R compute byte-identical inputs without either of them calling a random number generator or a transcendental function.
def lcg(n, seed=7):
"""A deterministic uniform sequence in [0, 1).
The arithmetic stays well inside a double's exact integer range, so this is the same sequence
in every language.
"""
s = seed
out = np.empty(n)
for i in range(n):
s = (75 * s + 74) % 65537
out[i] = s / 65537
return out
dates = np.arange("1990-01-01", "2020-01-01", dtype="datetime64[D]")
n = dates.size
u = lcg(n)
# Monthly base flow (m^3/s), a snowmelt-dominated regime peaking in spring.
base = np.array([18, 22, 35, 48, 40, 26, 16, 12, 11, 14, 20, 19], dtype=float)
month = dates.astype("datetime64[M]").astype(int) % 12
flow = np.round(base[month] * (0.6 + 0.8 * u), 3)
# Storms: the top 1.5% of the modulating sequence become peaks, each scaled by a second
# deterministic sequence so the peaks have a spread rather than all landing on one multiple.
storm = u > 0.985
flow[storm] = np.round(flow[storm] * (2.5 + 7 * lcg(n, seed=913)[storm]), 3)
print(n, flow.min(), flow.max())10957 6.602 588.076Two kinds of damage, both of which a real record has:
# A short instrument outage: five consecutive missing days in 2002.
outage = (dates >= np.datetime64("2002-06-10")) & (dates <= np.datetime64("2002-06-14"))
flow[outage] = np.nan
# Three days that never made it into the file at all.
dropped = np.isin(dates, np.array(["2004-03-02", "2004-03-03", "2005-11-19"],
dtype="datetime64[D]"))
gauge = TimeSeries(dates[~dropped], flow[~dropped], interval="one_day")
gauge<TimeSeries 10954 ordinates, interval "one_day", 1990-01-01T00:00:00 to 2019-12-31T00:00:00, 5 missing>The interval is not decoration. It is what shift(by="start") walks, what convert_interval() converts from, and – the one that matters below – what peaks_over_threshold() measures its independence criterion in. A “minimum of 5 steps between events” means five days here and would mean five months on a monthly series.
Repairing the record
Missing DATES and missing VALUES are different problems and take different methods. fill_missing_dates() inserts the ordinates that are absent entirely; interpolate_missing() fills short runs of values that are present but unobserved.
filled = gauge.fill_missing_dates(gauge.dates.min(), gauge.dates.max())
print("before:", len(gauge), " after:", len(filled))before: 10954 after: 10957Now every day is present and the three dropped ones are missing values, alongside the five-day outage. Interpolating with a limit of 3 fills the dropped days and leaves the outage alone – the limit is the point, not an inconvenience: a five-day gap in a flashy record is not something linear interpolation should be inventing values for.
repaired = filled.interpolate_missing(max_missing=3)
def observed_in(ts, start, end):
inside = (ts.dates >= np.datetime64(start)) & (ts.dates <= np.datetime64(end))
return int(np.count_nonzero(~np.isnan(ts.values[inside])))
{
"dropped days (2)": observed_in(repaired, "2004-03-02", "2004-03-03"),
"dropped day (1)": observed_in(repaired, "2005-11-19", "2005-11-19"),
"outage (5)": observed_in(repaired, "2002-06-10", "2002-06-14"),
}{'dropped days (2)': 2, 'dropped day (1)': 1, 'outage (5)': 0}What is in the record
stats = repaired.statistics()
{k: round(stats[k], 3) for k in ["Record Length", "Missing Values", "Minimum", "Maximum",
"Mean", "Std Dev", "Skewness", "50%", "95%"]}{'Record Length': 10957.0,
'Missing Values': 5.0,
'Minimum': 6.602,
'Maximum': 588.076,
'Mean': 26.131,
'Std Dev': 29.84,
'Skewness': 8.649,
'50%': 19.73,
'95%': 53.273}Per-month statistics show the regime directly. Note that this method filters missing values out, while monthly_percentiles() does not – a difference inherited from the library, and the reason June’s percentiles below are missing while its summary row is not. The columns of monthly_statistics() are minimum, 5%, 25%, 50%, 75%, 95%, maximum and mean.
np.round(repaired.monthly_statistics()[:, [0, 3, 7, 6]], 2)array([[ 10.82, 17.77, 19.41, 228.25],
[ 13.21, 22.16, 24.95, 264.65],
[ 21.03, 34.49, 41.38, 455.42],
[ 28.82, 47.99, 52.57, 588.08],
[ 24. , 39.98, 44.02, 523.22],
[ 15.61, 25.77, 29.17, 337.67],
[ 9.67, 15.6 , 17.31, 186.83],
[ 7.23, 11.94, 13.74, 156.89],
[ 6.6 , 10.96, 11.46, 136.09],
[ 8.4 , 13.71, 15.76, 185.23],
[ 12. , 20.3 , 23.09, 256.51],
[ 11.4 , 19.02, 21. , 249.83]])np.round(repaired.monthly_percentiles([0.5, 0.95])[4:7], 2)array([[39.98, 54.34],
[25.82, 35.55],
[15.6 , 21.81]])The duration curve is the same information as a flow-exceedance plot:
dur = repaired.duration()
np.round(dur[[0, 29, dur.shape[0] // 2, dur.shape[0] - 1]], 3)array([[9.00000e-03, 5.88076e+02],
[2.74000e-01, 3.12738e+02],
[5.00050e+01, 1.97300e+01],
[9.99910e+01, 6.60200e+00]])Annual maxima
A flood-frequency analysis wants one value per year, and in most of the world it wants the year to start in autumn so a winter flood season is not split across two of them. That is water_year():
ams = repaired.water_year(block="maximum", start_month=10)
list(zip(ams.dates[:8].astype("datetime64[D]").astype(str), ams.values[:8]))[(np.str_('1990-06-27'), np.float64(288.599)),
(np.str_('1991-03-11'), np.float64(380.373)),
(np.str_('1992-02-16'), np.float64(252.031)),
(np.str_('1993-04-15'), np.float64(459.903)),
(np.str_('1994-03-07'), np.float64(453.897)),
(np.str_('1995-03-17'), np.float64(281.437)),
(np.str_('1996-03-03'), np.float64(372.686)),
(np.str_('1997-03-17'), np.float64(276.662))]Two things to notice. The result is an IRREGULAR series carrying each maximum’s own date, not the end of its year – so you can see when the annual flood happened, which is a seasonality signal in itself. And the record spans thirty-one water years, not thirty: it starts in January 1990, part-way through water year 1990, and ends in December 2019, part-way through water year 2020. Both partial years contribute their own maximum.
Compare it with a calendar-year reduction:
{
"water_year": np.round(repaired.water_year().values, 1),
"calendar_year": np.round(repaired.calendar_year().values, 1),
}{'water_year': array([288.6, 380.4, 252. , 459.9, 453.9, 281.4, 372.7, 276.7, 254.1,
228.3, 376.7, 238.6, 192.2, 252.4, 448.2, 186.3, 441.1, 313.1,
402.3, 525.6, 347.9, 469.6, 288.6, 103.6, 500.7, 202.8, 202. ,
588.1, 473. , 523.2, 215.3]),
'calendar_year': array([288.6, 380.4, 252. , 459.9, 453.9, 281.4, 372.7, 276.7, 254.1,
234.3, 376.7, 238.6, 237.7, 252.4, 448.2, 186.3, 441.1, 313.1,
402.3, 525.6, 347.9, 469.6, 288.6, 66.1, 500.7, 202.8, 202. ,
588.1, 473. , 523.2])}They differ in both length and content: a flood in October or later belongs to the NEXT water year but the same calendar year.
One behaviour is worth stating outright, because it decides whether a gappy record silently loses a year. A block containing a missing value is DROPPED by the sum and average block functions – the missing value propagates through the arithmetic – but KEPT by minimum and maximum, because the comparison is simply false against a missing value and it is never selected. This is the library’s behaviour, and it is why the maximum series above spans every year despite the outage:
{
"maximum_years": len(filled.water_year(block="maximum")),
"sum_years": len(filled.water_year(block="sum")),
}{'maximum_years': 31, 'sum_years': 28}Fitting is one call away:
fit = fit_mle(model_univariate("GeneralizedExtremeValue", ams.values))
gev_params = list(fit.parameters.values())
print(np.round(gev_params, 4))
np.round(Distribution("GeneralizedExtremeValue", gev_params).quantile([0.5, 0.9, 0.99]), 2)[2.953528e+02 1.167666e+02 2.544000e-01]array([336.21, 495.43, 611.95])Peaks over a threshold
Annual maxima throw away most of a record: thirty-one numbers out of nearly eleven thousand. A partial-duration series keeps every independent peak above a threshold instead, which is why the independence criterion matters – without it, one flood’s rising and falling limbs count as many events.
pot_1 = repaired.peaks_over_threshold(threshold=120, min_steps_between_events=1)
pot_10 = repaired.peaks_over_threshold(threshold=120, min_steps_between_events=10)
{"one_day_apart": len(pot_1), "ten_days_apart": len(pot_10)}{'one_day_apart': 124, 'ten_days_apart': 106}list(zip(pot_10.dates[:8].astype("datetime64[D]").astype(str),
np.round(pot_10.values[:8], 2)))[(np.str_('1990-06-27'), np.float64(288.6)),
(np.str_('1990-11-07'), np.float64(138.3)),
(np.str_('1990-11-23'), np.float64(184.79)),
(np.str_('1990-12-27'), np.float64(152.43)),
(np.str_('1991-03-11'), np.float64(380.37)),
(np.str_('1991-03-26'), np.float64(340.14)),
(np.str_('1991-04-21'), np.float64(371.64)),
(np.str_('1992-02-16'), np.float64(252.03))]min_steps_between_events is counted in the series’ own interval, so on this daily record ten steps is ten days. The extracted series feeds a point-process model directly:
pp = model_point_process(pot_10.values, threshold=120, total_years=30)
pp_params = list(fit_mle(pp).parameters.values())
np.round(pp_params, 4)array([2.957606e+02, 1.171909e+02, 2.657000e-01])Smoothing happens BEFORE the extraction, which is how an n-day event series is built: the peaks of a 3-day moving average are the largest 3-day volumes, not the largest instantaneous flows.
pot_3day = repaired.peaks_over_threshold(threshold=120, min_steps_between_events=10,
smoothing="moving_average", period=3)
{"instantaneous": len(pot_10), "three_day_mean": len(pot_3day)}{'instantaneous': 106, 'three_day_mean': 30}Far fewer events clear the same threshold, which is the expected result rather than a surprise: a three-day mean of a flashy record is a smaller number than its instantaneous peak. Choose the threshold for the series you are actually extracting from.
Decomposition
seasonal_decompose() splits the record into trend, seasonal and residual parts, keeping only the harmonics of the seasonal frequency for the seasonal component. On a monthly aggregation of this record the annual cycle is the whole story:
monthly_blocks = repaired.block_series(window="month", block="average")
monthly = TimeSeries(monthly_blocks.dates, monthly_blocks.values, interval="one_month")
d = monthly.seasonal_decompose(period=12)
defined = ~np.isnan(d["trend"])
rows = np.flatnonzero(defined)[:4]
{
"date": d["date"][rows].astype("datetime64[M]").astype(str).tolist(),
"trend": np.round(d["trend"][rows], 3).tolist(),
"seasonal": np.round(d["seasonal"][rows], 3).tolist(),
"residual": np.round(d["residual"][rows], 3).tolist(),
}{'date': ['1990-12', '1991-01', '1991-02', '1991-03'],
'trend': [26.149, 26.218, 26.25, 28.326],
'seasonal': [-11.293, -6.409, 3.929, 12.113],
'residual': [8.598, -1.459, -7.217, 20.96]}The trend is a moving average over period ordinates, so it – and the residual with it – is undefined for the first eleven months. Where all three components exist they add back to the original value exactly, which the reproduction check at the bottom asserts.
Synthetic records
Two resamplers generate records of arbitrary length from the observed one. resample_knn() is the conditional k-nearest-neighbour bootstrap of Lall and Sharma (1996): at each step it finds the k historical days closest to the current flow, picks one, and advances to whatever historically came NEXT. That last step is what preserves the recession behaviour – taking the neighbour’s own value instead would leave the trajectory hovering near where it started.
observed = repaired.replace_missing(0.0)
knn = observed.resample_knn(time_steps=365, k=20, seed=2024)
block = observed.resample_block_bootstrap(time_steps=365, block_size=30, seed=2024)
def lag1(ts):
return autocorrelation(ts, max_lag=1)["value"][1]
{
"observed": (round(observed.values.mean(), 3), round(observed.values.std(ddof=1), 3),
round(lag1(observed), 3)),
"knn": (round(knn.values.mean(), 3), round(knn.values.std(ddof=1), 3), round(lag1(knn), 3)),
"block": (round(block.values.mean(), 3), round(block.values.std(ddof=1), 3),
round(lag1(block), 3)),
}{'observed': (np.float64(26.119), np.float64(29.839), np.float64(0.194)),
'knn': (np.float64(23.174), np.float64(27.637), np.float64(0.133)),
'block': (np.float64(22.924), np.float64(28.45), np.float64(0.153))}Both preserve the marginal distribution to within sampling error on a 365-day draw, and both recover most of the lag-1 dependence. Neither recovers all of it, and it is worth being clear about why rather than reading the table as a ranking: this record’s day-to-day persistence comes almost entirely from its seasonal cycle rather than from a recession, since the modulating sequence is independent from day to day. A resampler that reorders time cannot keep a seasonal signal it was never conditioned on. On a real record with genuine recession behaviour the kNN resample is the one that holds lag-1 dependence, because it conditions each step on the current state; the block bootstrap keeps dependence only WITHIN a block and breaks it at every boundary.
Seeded runs are reproducible, and the reproduction is not only within one language: the whole computation, the Mersenne Twister included, happens in the shared C++ core, so these 365 numbers are bit-identical to the ones the R page prints.
np.round(knn.values[:8], 3)array([34.272, 53.309, 49.187, 16.272, 14.897, 17.273, 18.879, 9.222])Reproduction check
Every literal below was produced by this page. The R page asserts the same values.
def near(a, b, tol=1e-15):
return abs(a - b) <= tol * max(1.0, abs(b))
defined = ~np.isnan(d["trend"])
gev = gev_params
assert len(gauge) == 10954
assert len(filled) == 10957
assert int(np.count_nonzero(np.isnan(filled.values))) == 8
assert int(np.count_nonzero(np.isnan(repaired.values))) == 5
assert near(stats["Record Length"], 10957)
assert near(stats["Missing Values"], 5)
assert near(stats["Mean"], 26.130572315558798)
assert near(stats["Maximum"], 588.07600000000002)
assert len(ams) == 31
assert len(repaired.calendar_year()) == 30
assert near(ams.values[0], 288.59899999999999)
assert near(ams.values.max(), 588.07600000000002)
assert len(filled.water_year(block="maximum")) == 31
assert len(filled.water_year(block="sum")) == 28
assert near(gev[0], 295.35281082558026)
assert near(gev[1], 116.76657650129897)
assert near(gev[2], 0.25436343374399251)
assert len(pot_1) == 124
assert len(pot_10) == 106
assert len(pot_3day) == 30
assert near(pot_10.values[0], 288.59899999999999)
assert near(pp_params[0], 295.7606336759473)
assert len(monthly) == 359
assert int(np.count_nonzero(defined)) == 348
assert np.max(np.abs(d["trend"][defined] + d["seasonal"][defined] + d["residual"][defined]
- monthly.values[defined])) < 1e-9
assert len(knn) == 365
assert len(block) == 365
assert near(knn.values[0], 34.271999999999998)
assert near(knn.values[7], 9.2219999999999995)
assert near(block.values[0], 23.332999999999998)
print("All reproduction checks passed.")All reproduction checks passed.