Import an annual maximum (AMAX) sample from NRFA peak flow .am files

Description

Imports the peak flows and dates from from NRFA peak flow .am files, excluding the rejected years

Usage

AMImport(x)

Arguments

Value

A data.frame with columns; Date and Flow

Author(s)

Anthony Hammond

Examples

# Import an AMAX sample and display the first six rows in the console
## Not run: 
am_4003 <- AMImport(r"{C:\Data\NRFAPeakFlow_v11\Suitable for QMED\4003.AM}")

## End(Not run)
## Not run: 
head(am_4003)

## End(Not run)

National River Flow Archive (NRFA) annual maximum peak flow data for sites suitable for QMED only, and those for pooling.

Description

A data.frame with three columns; Date, Flow, id. NRFA Peak Flow Dataset - Version (14)

Usage

AMPF

Format

A data frame with 43064 rows and 3 columns

Date

Date

Flow

Annual maximum peak flow, m3/s

id

Station identification number

Source

https://nrfa.ceh.ac.uk/data/peak-flow-dataset

Plot of the annual maximum sample

Description

Provides a barplot for an annual maximum sample

Usage

AMplot(x, ylab = "Discharge (m3/s)", xlab = "Hydrological year", main = NULL)

Arguments

Details

When used with a GetAM object or any data.frame with dates/POSIXct in the first column, the date-times are daily or sub-daily. Therefore, although it's an annual maximum (AM) sequence, some bars may be closer together depending on the number of days between them.

Value

A barplot of the annual maximum sample

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and plot
AMplot(GetAM(58002))

Areal reduction factor (ARF)

Description

The results of applying, to a rainfall depth, the ratio of the rainfall over an area to the rainfall depth of the same duration at a representative point in the area.

Usage

ARF(Depth, Area, D)

Arguments

Details

The ARF and it's use is detailed in the Flood Estimation Handbook (1999), volume 2. The DDF model is calibrated on point rainfall and the areal reduction factor converts it to a catchment rainfall for use with a rainfall runoff model such as ReFH (see details for ReFH function). The ReFH model includes a design rainfall profile for winter and summer but the depth duration frequency (DDF) model is calibrated on annual maximum peaks as opposed to seasonal peaks. A seasonal correction factor (SCF) is necessary to convert the DDF estimate to a seasonal one. The final depth, therefore is; Depth = DDFdepth x ARF x SCF.

Value

the rainfall depth or rainfall return period

Author(s)

Anthony Hammond

Examples

# Derive the ARF for a depth of 30, an area of 500km^2 and a duration of 12 hours
ARF(30, 500, 12)

Add an AMAX sample

Description

This function allows the user to add an AMAX sample and associated catchment descriptors for use with the FEH process.

Usage

AddGauge(CDs, AMAX, ID)

Arguments

Details

The function provides the necessary AMAX sample statistics and data.frame for adding catchment descriptors to the PeakFlowData data.frame. The user must then add these outputs using the rbind function (see example). The AMAX could be read in or pasted in by the user or imported using the AMImport function. Once they are added they can be used in the current R session. If a new session is started (rather than a saved workspace) the added AMAX would need to be added again.

Value

A list object. The first element is a data.frame which is a row of statistics and descriptors to be added to the PeakFlowData data.frame. The second element is the AMAX sample formatted to be added to the AMPF data.frame

Author(s)

Anthony Hammond

Examples

# Read in AMAX and catchment descriptors
## Not run: 
am_add <- AMImport(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.am}")
cds_add <- CDsXML(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.xml}")

## End(Not run)

# Apply the function and add the results to the necessary data frames
## Not run: 
gauge_27003 <- AddGauge(cds_add, am_add, ID = "27003")

## End(Not run)

# Append the descriptors and stats (element[[1]]) to PeakFlowData
## Not run: 
nrfa_data <- rbind(PeakFlowData, gauge_27003[[1]])

## End(Not run)

# Append the AMAX series (element[[2]]) to AMPF
## Not run: 
ampf <- rbind(AMPF, gauge_27003[[2]])

## End(Not run)

Aggregate a time series

Description

Aggregates time series data, creating hourly data from 15-minute data for example.

Usage

AggDayHour(x, func, Freq = "Day", hour = 9)

Arguments

Details

The function can be used with a data.frame with POSIXct in the first column and a variable in the second. You can choose the level of aggregation in hours, or you can choose daily. In the daily case you can choose which hour of the day to start the aggregation. For example, you might want mean flows from 09:00 rather than midnight. You can also choose the function used to aggregate the data. For example, you might want "sum" for rainfall, and "mean" for flow. When aggregating hourly the aggregation starts at whatever hour is in the first row of x and the associated time stamps will reflect this.

Value

A data.frame with POSIXct in the first column (unless daily is chosen, then it's Date class), and the aggregated variable in the second column

Author(s)

Anthony Hammond

Examples

# Create a data frame with a normally distributed variable at
# a 15 minute sampling rate
ts_seq <- seq(
  as.POSIXct("2000-01-01 00:00:00", tz = "Europe/London"),
  as.POSIXct("2001-01-01 00:00:00", tz = "Europe/London"),
  by = 60 * 15
)
ts_df <- data.frame(DateTime = ts_seq, Var = rnorm(length(ts_seq), 10, 2))

# Aggregate to an hourly sampling rate, taking the maximum of each hour
hourly <- AggDayHour(ts_df, func = max, Freq = "Hour")

# Aggregate with the mean at a daily scale
daily <- AggDayHour(ts_df, func = mean, Freq = "Day")

Annual statistics extraction

Description

Extracts annual statistics (default maximums) from a data.frame which has dates (or POSIXct) in the first column and variable in the second.

Usage

AnnualStat(
  x,
  Stat = max,
  Truncate = TRUE,
  Mon = 10,
  Hr = 9,
  Sliding = FALSE,
  N = 24,
  ...
)

Arguments

Details

The statistics are extracted based on the UK hydrological year by default (start month = 10). Month can be changed using the Mon argument. A year is from Mon-Hr to Mon-(Hr-1). For example, the 2018 hydrological year with Hr = 9 would be from 2018-10-01 09:00:00 to 2019-10-01 08:00:00. If Hr = 0, then it would be from 2018-10-01 00:00:00 to 2019-09-30 23:00:00. Data before the first occurrence of the 'start month' and after and including the last occurrence of the 'start month' is not included in the calculation of the statistic.

Value

a data.frame with columns; DateTime and Result. By default Result is the annual maximum sample, but will be any statistic used as the Stat argument.

Author(s)

Anthony Hammond

Examples

# Extract the Thames AMAX daily mean flow and display the first six rows
thames_am <- AnnualStat(ThamesPQ[, c(1, 3)])
head(thames_am)

# Extract the annual rainfall totals
thames_annual_p <- AnnualStat(ThamesPQ[, 1:2], Stat = sum)

# Extract maximum five day rainfall totals from the Thames rainfall series
thames_5day_am <- AnnualStat(ThamesPQ[, 1:2], Stat = sum, Sliding = TRUE, N = 5)

Baseflow index (BFI)

Description

Calculates the baseflow index from a daily mean flow series

Usage

BFI(Q, PlotTitle = "Baseflow plot", Plot = TRUE, ReturnData = FALSE)

Arguments

Details

The baseflow index is calculated using the method outlined in Gustard, A. Bullock, A. Dixon, J. M.. (1992). Low flow estimation in the United Kingdom. Wallingford, Institute of Hydrology, 88pp. (IH Report No.108)

Value

the baseflow index and if Plot equals TRUE, a plot showing the flow time series (black) and the associated baseflow (red). If ReturnData is set to TRUE, the result is a list with two elements. The first is the BFI, the second is a data frame with the flow data and baseflow data.

Author(s)

Anthony Hammond

Examples

# Calculate the BFI from the daily discharge at Kingston upon Thames,
# which is in column three of the ThamesPQ data
BFI(ThamesPQ[, 3])

Bootstrap

Description

Resampling with replacement to approximate the sampling distribution of a statistic and quantify uncertainty.

Usage

Bootstrap(x, Stat, n = 500, Conf = 0.95, ReturnSD = FALSE, ...)

Arguments

Details

The bootstrapping procedure resamples from a sample length(x) * n times with replacement. After splitting into n samples of size length(x), the statistic of interest is calculated on each.

Value

If ReturnSD is FALSE a data.frame is returned with one row and three columns; central, lower, and upper. If ReturnSD is TRUE, the sampling distribution is returned.

Author(s)

Anthony Hammond

Examples

# Extract an AMAX sample and quantify uncertainty for the Gumbel estimated 50-year flow
am_203018 <- GetAM(203018)
Bootstrap(am_203018$Flow, Stat = GumbelAM, RP = 50)

# Quantify uncertainty for the sample standard deviation at the 90% confidence level
Bootstrap(am_203018$Flow, Stat = sd, Conf = 0.90)

# Return the sampling distribution of the mean and plot an associated histogram
samp_dist <- Bootstrap(am_203018$Flow, Stat = mean, ReturnSD = TRUE)
hist(samp_dist)

Import catchment descriptors from .xml files

Description

Imports catchment descriptors from xml files either from an FEH webservice download or from the Peakflows dataset downloaded from the national river flow archive (NRFA) website

Usage

CDsXML(x)

Arguments

Value

A data.frame with columns; Descriptor and Value.

Author(s)

Anthony Hammond

Examples

# Import catchment descriptors from a FEH webserver XML file and display XML in the console
## Not run: 
cds_my_site <- CDsXML(r"{C:\Data\FEH_Catchment_384200_458200.xml}")
cds_my_site

## End(Not run)

Import catchment descriptors from older .xml files

Description

Imports catchment descriptors from xml files (prior to FEH2025) either from an FEH webservice download or from the Peakflows dataset downloaded from the national river flow archive (NRFA) website

Usage

CDsXML_Legacy(x)

Arguments

Details

This function is to allow users to import catchment descriptors in the format prior to the 2025 update.

Value

A data.frame with columns; Descriptor and Value.

Author(s)

Anthony Hammond

Examples

# Import catchment descriptors from a FEH webserver XML file and display XML in the console
## Not run: 
cds_my_site <- CDsXML(r"{C:\Data\FEH_Catchment_384200_458200.xml}")
cds_my_site

## End(Not run)

Convert between British National Grid Reference (BNG) and Latitude and Longitude or Irish Grid references.

Description

Function to convert between BNG easting & northing and Latitude & Longitude (or vice versa). Or to convert between BNG and Irish national grid (or vice versa)

Usage

ConvertGridRef(x, fromBNG = TRUE, IGorLatLon = "LatLon")

Arguments

Details

To convert to Lat and Lon from BNG, ensure that the fromBNG argument is TRUE. To convert the other way, set fromBNG as FALSE. The same applies for converting between Irish grid and BNG. To convert Irish grid and BNG set the IGorLatLon argument to IG.

Value

A data.frame with the converted grid references. Either latitude and longitude if BNG = TRUE. Or BNG easting and northing if fromBNG = FALSE. Or, IG easting & northing if fromBNG = TRUE and IGorLatLon = "IG".

Author(s)

Anthony Hammond

Examples

# Convert a BNG numeric reference to latitude and longitude
ConvertGridRef(c(462899, 187850))

# Convert latitude and longitude to easting and northing
ConvertGridRef(c(51.6, -1), fromBNG = FALSE)

DDF results from a DDFImport object

Description

Extracts results from a data frame imported using the DDFImport function

Usage

DDF(x, duration, RP = 100)

Arguments

Details

The .xml files only provide a set number of durations and return periods for DDF13 and DDF22. This is an interpolator function to derive depths for intervening durations and return periods. The result is rounded to an integer.

Value

A DDF13 or DDF22 estimate of rainfall depth (mm)

Author(s)

Anthony Hammond

Examples

# Import DDF13 results from an NRFA Peak Flows XML file
## Not run: 
ddf13_4003 <- DDFImport("C:/Data/NRFAPeakFlow_v9/Suitable for QMED/04003.xml", DDFVersion = 13)

## End(Not run)

# Estimate the 20-year, 5-hour depth
## Not run: 
DDF(ddf13_4003, duration = 5, RP = 20)

## End(Not run)

FEH99 depth duration frequency precipitation model

Description

Estimation of design rainfall depths, and the rarity of observed rainfall

Usage

DDF99(Duration, RP, pars, Depth = NULL, disc = NULL)

Arguments

Details

The depth duration frequency rainfall model is detailed in the Flood Estimation Handbook (1999), volume 2. A note about the discretisation: The user can choose between "daily" or "hourly" for the sliding duration to fixed duration conversion. If the 'Depth' argument is used, it overrides the return period (RP) argument and provides RP as a function of depth. However, if both the 'Depth' and the 'disc' arguments are used, the sliding duration depth is provided as a function of the user input depth. This resulting depth can then be used without the 'disc' argument to determine the sliding duration RP.

Value

the rainfall depth or rainfall return period

Author(s)

Anthony Hammond

Examples

# Examples from FEH volume 2
# The parameters for these examples are from FEH v2

# What is the 2-day rainfall with return period 100 years for Norwich?
DDF99(Duration = 48, RP = 100, pars = c(-0.023, 0.273, 0.351, 0.236, 0.309, 2.488))

# What is the 4-hour rainfall with return period 20 years for a typical point in the Lyne catchment?
DDF99(Duration = 4, RP = 20, pars = c(-0.025, 0.344, 0.485, 0.402, 0.287, 2.374))

# How rare was the rainfall of 6th August 1978 at Broughshane, County Antrim?
DDF99(Duration = 5, Depth = 47.7, pars = c(-0.022, 0.412, 0.551, 0.276, 0.261, 2.252))

DDF99 parameters from .xml files

Description

Imports the FEH 1999 depth duration frequency parameters from xml files either from an FEH webservice download or from the Peakflows dataset downloaded from the national river flow archive (NRFA) website

Usage

DDF99Pars(x)

Arguments

Details

This function is coded to import DDF99 parameters from xml files from the NRFA or the FEH web-server. File paths for importing data require forward slashes. On some operating systems, such as windows, the copy and pasted file paths will have backward slashes and would need to be changed accordingly.

Value

A list with two elements, each a data frame with columns; parameters and associated values The first data frame is for the catchment average parameters (these still require an ARF adjustment where appropriate) and the second for the 1km2 grid point parameters.

Author(s)

Anthony Hammond

Examples

# Import DDF99 parameters from an NRFA Peak Flows XML file and display in console
## Not run: 
ddf99_4003 <- DDF99Pars("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/04003.xml")
ddf99_4003

## End(Not run)

# Import DDF99 parameters from a FEH webserver XML file and display in the console
## Not run: 
ddf99_my_site <- DDF99Pars("C:/Data/FEH_Catchment_384200_458200.xml")
ddf99_my_site

## End(Not run)

Derive and plot rainfall Depth Duration Frequency curves.

Description

Derive and plot rainfall Depth Duration Frequency curves from an input of rainfall data.

Usage

DDFExtract(x, Plot = TRUE, main = NULL, Truncate = TRUE)

Arguments

Details

The function works by extracting the annual maximum sample (by hydrological year - starting Oct 1st) of rainfall for a range of sliding durations (1 hour to 96 hours). It then calculates the median annual maximum rainfall depth (RMED) and a GEV growth curve for each duration. To ensure RMED increases with duration a power curve is fit as a function of duration to provide the final RMED estimates. Then the average growth factor for each return period (across the durations) is assumed.

Value

A dataframe with hours (1 to 96) in the first column then depths associated with a range of return periods (2 to 1000) in the remaining nine columns. If Plot = TRUE, a plot of the DDF curves is also returned.

Author(s)

Anthony Hammond

Examples

# Extract all available 15-minute rainfall from the St Ives (Cambridgeshire)
# rain gauge (WISKI ID = 179365).
## Not run: 
st_ives <- GetDataEA_Rain(WISKI_ID = "179365", Period = "15Mins")

## End(Not run)

# Apply the DDF function.
## Not run: 
DDFExtract(st_ives)

## End(Not run)

DDF13 or DDF22 results from .xml files

Description

Imports the depth duration frequency 2013 or 2022 results from xml files either from an FEH webservice download or from the Peakflows dataset downloaded from the national river flow archive (NRFA) website

Usage

DDFImport(x, ARF = TRUE, Plot = TRUE, DDFVersion = 22)

Arguments

Details

This function returns a data-frame of design rainfall estimates. For further durations and return periods, the separate DDF function can be applied with the data-frame as the argument/input.

Value

A data frame of DDF results (mm) with columns for duration and rows for return period. If Plot equals TRUE a DDF plot is also returned.

Author(s)

Anthony Hammond

Examples

# Import DDF22 results from an NRFA Peak Flows XML file and display them in console
## Not run: 
ddf22_4003 <- DDFImport(r"{C:\Data\NRFAPeakFlow_v11\Suitable for QMED\04003.xml}")
ddf22_4003

## End(Not run)

# Import DDF22 results from a FEH webserver XML file and display them in the console
## Not run: 
ddf22_my_site <- DDFImport(r"{C:\Data\FEH_Catchment_384200_458200.xml}")
ddf22_my_site

## End(Not run)

Linearly detrend a sample

Description

Applies a linear detrend to a sample

Usage

DeTrend(x)

Arguments

Details

Adjusts all the values in the sample, of size n, by the difference between the linearly modelled ith data point and the linearly modelled nth data point.

Value

A numeric vector which is a linearly detrended version of x.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum (AM) sample that looks to have a significant trend
am_21025 <- GetAM(21025)

# Plot the resulting AM as a bar plot. Then detrend and compare with another plot
plot(am_21025$Flow, type = "h", ylab = "Discharge (m^3/s)")
am_detrend <- DeTrend(am_21025$Flow)
plot(am_detrend, type = "h", ylab = "Discharge (m^3/s)")

Design hydrograph extraction

Description

Extracts a mean hydrograph from a flow series

Usage

DesHydro(
  x,
  Threshold = 0.975,
  EventSep,
  N = 10,
  Exclude = NULL,
  Plot = TRUE,
  main = "Design Hydrograph"
)

Arguments

Details

All the peaks over the threshold (default 0.975th) are identified and separated by a user defined value 'EventSep', which is a number of timesteps (peaks are separated by EventSep * 3). The top N peaks are selected and the hydrographs are then extracted. The hydrograph start is the time of peak minus EventSep. The End of the hydrograph is time of peak plus EventSep times 1.5. All events are scaled to have a peak flow of one, and the mean of these is taken as the scaled design hydrograph.

Value

a list of length three. The first element is a dataframe of the peaks of the hydrographs and the associated dates. The second element is a dataframe with all the scaled hydrographs, each column being a hydrograph. The third element is the averaged hydrograph

Author(s)

Anthony Hammond

Examples

# Extract a design hydrograph from the Thames daily mean flow and print the resulting hydrograph
thames_des_hydro <- DesHydro(ThamesPQ[, c(1, 3)], EventSep = 10, N = 10)

Diagnostic plots for pooling groups

Description

Provides 11 plots to compare the sites in the pooling group

Usage

DiagPlots(x, gauged = FALSE, UrbMax = 0.03)

Arguments

Value

Eleven diagnostic plots for pooling groups

Author(s)

Anthony Hammond

Examples

# Form a gauged pooling group and plot the diagnostics
pool_28015 <- Pool(GetCDs(28015))
DiagPlots(pool_28015, gauged = TRUE)

# Form an ungauged pooling group and plot the diagnostics
pool_28015 <- Pool(GetCDs(28015), exclude = 28015)
DiagPlots(pool_28015)

Donor adjustment candidates

Description

Provides donor adjustment candidates, with associated descriptors, in order of the proximity to the centroid of the subject catchment.

Usage

DonAdj(CDs, N = 10, UrbMax = 1)

Arguments

Details

This function provides the closest N catchments for consideration for QMED donor adjustment.

Value

A data.frame with rownames of site references and columns of catchment descriptors and distance from subject site.

Author(s)

Anthony Hammond

Examples

# Get some CDs and output candidate donor sites
cds_54022 <- GetCDs(54022)
DonAdj(cds_54022)

Extreme rank plot

Description

A plot to inspect the distribution of ordered data

Usage

ERPlot(x, dist = "GenLog", main = NULL, Pars = NULL, GF = NULL, ERType = 1)

Arguments

Details

By default this plot compares the percentage difference of simulated results with observed for each rank of the data. Another option (see ERType argument) compares the simulated flows for each rank of the sample with the observed of the same rank. For both plots 500 simulated samples are used. With the second option for each rank they are plotted and the mean of these is highlighted in red. There is a line of perfect fit so you can see how much this "cloud" of simulation differs from the observed. By default the parameters of the distribution for comparison with the sample are estimated from the sample. However, the pars argument can be used to compare the distribution with parameters estimated separately. Similarly the growth factor (GF) parameters, linear coefficient of variation (Lcv) & linear skewness (LSkew) with the median can be entered. In this way the pooling estimated distribution can be compared to the sample. This ERplot is an updated version of that described in Hammond, A. (2019). Proposal of the 'extreme rank plot' for extreme value analysis: with an emphasis on flood frequency studies. Hydrology Research, 50 (6), 1495-1507.

Value

The extreme rank plot as described in the details

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and plot
## Not run: 
am_27083 <- GetAM(27083)
ERPlot(am_27083$Flow)

## End(Not run)

# Assume pooled estimates of Lcv and LSkew of 0.23, 0.17, and a QMED of 12.
# Use these inputs for the GF argument and change the title
## Not run: 
ERPlot(am_27083$Flow, main = "Site 27083 pooled comparison", GF = c(0.23, 0.17, 12))

## End(Not run)

Extreme value plot (frequency and growth curves)

Description

Plots the extreme value frequency curve or growth curve with observed sample points.

Usage

EVPlot(
  x,
  dist = "GenLog",
  scaled = TRUE,
  Title = "Extreme value plot",
  ylabel = NULL,
  LineName = NULL,
  Unc = TRUE
)

Arguments

Details

The plotting has the option of generalised extreme value (GEV), generalised Pareto (GenPareto), Gumbel, or generalised logistic (GenLog) distributions. The uncertainty is quantified by bootstrapping.

Value

An extreme value plot (frequency or growth curve) with intervals to quantify uncertainty

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and plot the growth curve with the GEV distribution
am_203018 <- GetAM(203018)
EVPlot(am_203018$Flow, dist = "GEV")

Add lines and/or points to an extreme value plot

Description

Functionality to add extra lines or points to an extreme value plot (derived from the EVPlot function).

Usage

EVPlotAdd(
  Pars,
  dist = "GenLog",
  Name = "Adjusted",
  MED = NULL,
  xyleg = NULL,
  col = "red",
  lty = 1,
  pts = NULL,
  ptSym = NULL
)

Arguments

Details

A line can be added using the Lcv and Lskew based on one of four distributions (Generalised extreme value, Generalised logistic, Gumbel, Generalised Pareto). Points can be added as a numeric vector. If a single point is required, the base points() function can be used and the x axis will need to be log(RP-1).

Value

Additional, user specified line or points to an extreme value plot derived from the EVPlot function.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and plot the growth curve with the GEV distribution
am_203018 <- GetAM(203018)
EVPlot(am_203018$Flow, dist = "GEV")

# Now add a line (dotted and red) for the generalised logistic distribution
# First get the Lcv and Lskew using the L-moments function
pars <- as.numeric(LMoments(am_203018[, 2])[c(5, 6)])
EVPlotAdd(
  Pars = pars, dist = "GenLog", Name = "GenLog",
  xyleg = c(-5.2, 2.65), lty = 3
)

# Now add a line for the Gumbel distribution which is dark green and dashed
EVPlotAdd(
  Pars = pars[1], dist = "Gumbel", Name = "Gumbel",
  xyleg = c(-5.19, 2.5), lty = 3, col = "darkgreen"
)

# Now plot afresh and get another AMAX and add the points
EVPlot(am_203018$Flow, dist = "GEV")
am_27090 <- GetAM(27090)
EVPlotAdd(xyleg = c(-4.9, 2.65), pts = am_27090[, 2], Name = "27090")

Extreme value plot for pooling groups

Description

Plots the extreme value frequency curve or growth curve for gauged or ungauged pooled groups

Usage

EVPool(
  x,
  AMAX = NULL,
  gauged = FALSE,
  dist = "GenLog",
  QMED = NULL,
  Title = "Pooled growth curve",
  UrbAdj = FALSE,
  CDs
)

Arguments

Value

An extreme value plot for gauged or ungauged pooling groups

Author(s)

Anthony Hammond

Examples

# Get some CDs, form an ungauged pooling group and apply EVPool
cds_28015 <- GetCDs(28015)
pool_28015 <- Pool(cds_28015, exclude = 28015)
EVPool(pool_28015)

# Do the same for the gauged case, change the title, and convert with a QMED of 9.8
pool_g_28015 <- Pool(cds_28015)
EVPool(pool_g_28015, gauged = TRUE, Title = "Gauged frequency curve - Site 28015", QMED = 9.8)

Encounter probabilities

Description

Calculates the probability of experiencing at least n events with a given return period (RP), over a given number of years

Usage

EncProb(n, yrs, RP, dist = "Poisson")

Arguments

Details

The choice of binomial or Poisson distributions for calculating encounter probablities is akin to annual maximum (AM) versus peaks over threshold (POT) approaches to extreme value analysis. AM and binomial assume only one "event" can occur in the blocked time period. Whereas Poisson and POT don't make this assumption. In the case of most catchments in the UK, it is rare to have less than two independent "events" per year; in which case the Poisson and POT choices are more suitable. In large catchments, with seasonally distinctive baseflow, there may only be one independent peak in the year. However, the results from both methods converge with increasing magnitude, yielding insignificant difference beyond a 20-year return period.

Value

A probability

Author(s)

Anthony Hammond

Examples

# Calculate the probability of exceeding at least one 50-year RP event
# over a 10-year period, using the Poisson distribution
EncProb(n = 1, yrs = 10, RP = 50)

# Calculate the probability of exceeding at least two 100-year RP events
# over a 100-year period, using the Binomial distribution
EncProb(n = 2, yrs = 100, RP = 100, dist = "Binomial")

Flow duration curve

Description

A function to plot flow duration curves for a single flow series or flow duration curves from multiple flow series.

Usage

FlowDurationCurve(
  x = NULL,
  main = "Flow duration curve",
  CompareCurves = NULL,
  LegNames = NULL,
  Cols = NULL,
  AddQs = NULL,
  ReturnData = FALSE
)

Arguments

Details

The user can input a dataframe of dates (or POSIXct) and flow to return a plot of the flow duration curve for annual, winter and summer periods. Alternatively a list of flow series' (vectors) can be applied for a plot comparing the individual flow duration curves.

Value

If a dataframe of date in the first column and flow in the second is applied with the x argument a plot of the flow duration curves for the winter, summer and annual periods is returned. If a list of flow series is applied with the CompareCurves argument the associated flow duration curves are all plotted together. If ReturnData is TRUE, the plotted data is also returned.

Author(s)

Anthony Hammond

Examples

# Plot a flow duration curve for the Thames at Kingston October 2000 to September 2015
FlowDurationCurve(ThamesPQ[, c(1, 3)])

# Add two additional flow lines for the plot
FlowDurationCurve(ThamesPQ[, c(1, 3)], AddQs = c(25, 200))

# Compare flows from the rather wet 2013 water year (rows 4749 and 5114) with the rest of the flow
FlowDurationCurve(
  CompareCurves = list(
    ThamesPQ$Q[-seq(4749, 5114)],
    ThamesPQ$Q[4749:5114]
  ),
  LegNames = c("All but 2013", "Water year 2013")
)

Flow splitter

Description

A function to separate baseflow from runoff.

Usage

FlowSplit(
  x,
  BaseQUpper = NULL,
  AdjUp = NULL,
  ylab = "Value",
  xlab = "Time index"
)

Arguments

Details

The function is intended for the event scale as opposed to long term flow series. It works by linearly joining all the low points in the hydrograph - also the beginning and end points. Where a low point is any point with two higher points either side. Then any values above the hydrograph (xi) are set as xi. The baseflow point on the falling limb of the hydrograph/s can be raised using the AdjUp argument. The function works for any sampling frequency and arbitrary hydrograph length (although more suited for event scale and sub-annual events in general). This function is not design for deriving long term baseflow index. It could be used for such a purpose but careful consideration would be required for the BaseQUpper argument especially for comparison across river locations. If baseflow index is required the BFI function (with daily mean flow) is more suitable.

Value

A dataframe with the original flow (x) in the first column and the baseflow in the second. A plot of the original flow and the baseflow is also returned.

Author(s)

Anthony Hammond

Examples

# We'll extract a wet six month period on the Thames during the 2006-2007 hydrological year
thames_q <- subset(ThamesPQ[, c(1, 3)], Date >= "2006-11-04" & Date <= "2007-05-06")

# Then apply the flow split with default settings
q_split <- FlowSplit(thames_q$Q)

# Now do it with an upper baseflow level of 100 m^3/s
q_split <- FlowSplit(thames_q$Q, BaseQUpper = 100)


# Next we will get a single peaked "idealised" hydrograph using the ReFH function
q_refh <- ReFH(GetCDs(15006))
q_refh <- q_refh[[2]]$TotalFlow

# Now use the function with and without an upward adjustment of the baseflow on the falling limb
q_flow_split <- FlowSplit(q_refh)
q_flow_split <- FlowSplit(q_refh, AdjUp = 0.15)

Generalised extreme value distribution - estimates directly from sample

Description

Estimated quantiles as function of return period (RP) and vice versa, directly from the data

Usage

GEVAM(x, RP = 100, q = NULL)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The parameters are estimated by the method of L-moments, as detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the 50-year RP
am_27090 <- GetAM(27090)
GEVAM(am_27090$Flow, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GEVAM(am_27090$Flow, q = 600)

Generalised extreme value distribution estimates from parameters

Description

Estimated quantiles as function of return period (RP) and vice versa, from user input parameters

Usage

GEVEst(loc, scale, shape, q = NULL, RP = 100)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample, estimate the parameters, and estimate the 50-year RP
am_27090 <- GetAM(27090)
GEVPars(am_27090$Flow)

# Store the parameters in an object
pars <- as.numeric(GEVPars(am_27090$Flow))

# Get an estimate of 50-year flow
GEVEst(pars[1], pars[2], pars[3], RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GEVEst(pars[1], pars[2], pars[3], q = 600)

Generalised extreme value distribution growth factors

Description

Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)

Usage

GEVGF(lcv, lskew, RP)

Arguments

Details

Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.

Value

Generalised extreme value estimated growth factor

Author(s)

Anthony Hammond

Examples

# Estimate the 50-year growth factor from Lcv = 0.17 and Lskew = 0.04
GEVGF(0.17, 0.04, RP = 50)

Generalised extreme value distribution parameter estimates

Description

Estimated parameters from a sample (with Lmoments or maximum likelihood estimation) or from L1 (first L-moment), Lcv (linear coefficient of variation), and LSkew (linear skewness)

Usage

GEVPars(x = NULL, mle = FALSE, L1 = NULL, LCV = NULL, LSKEW = NULL)

Arguments

Details

The L-moment estimated parameters are by the method detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

Parameter estimates (location, scale, shape)

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the parameters using L-moments
am_27090 <- GetAM(27090)
GEVPars(am_27090$Flow)

# Estimate parameters using MLE
GEVPars(am_27090$Flow, mle = TRUE)

# Calculate L-moments and estimate the parameters with L1, Lcv, and Lskew
LMoments(am_27090$Flow)

# Store L-moments in an object
l_pars <- as.numeric(LMoments(am_27090$Flow))[c(1, 5, 6)]
GEVPars(L1 = l_pars[1], LCV = l_pars[2], LSKEW = l_pars[3])

Generalised logistic distribution - estimates directly from sample

Description

Estimated quantiles as a function of return period (RP) and vice versa, directly from the data

Usage

GenLogAM(x, RP = 100, q = NULL)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The parameters are estimated by the method of L-moments, as detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the 50-year RP
am_27090 <- GetAM(27090)
GenLogAM(am_27090$Flow, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GenLogAM(am_27090$Flow, q = 600)

Generalised logistic distribution estimates from parameters

Description

Estimated quantiles as function of return period (RP) and vice versa, from user input parameters

Usage

GenLogEst(loc, scale, shape, q = NULL, RP = 100)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample, estimate the parameters, and estimate 50-year RP
am_27090 <- GetAM(27090)
GenLogPars(am_27090$Flow)

# Store the parameters in an object
pars <- as.numeric(GenLogPars(am_27090$Flow))

# Get an estimate of 50-year flow
GenLogEst(pars[1], pars[2], pars[3], RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GenLogEst(pars[1], pars[2], pars[3], q = 600)

Generalised logistic distribution growth factors

Description

Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)

Usage

GenLogGF(lcv, lskew, RP)

Arguments

Details

Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.

Value

Generalised logistic estimated growth factor

Author(s)

Anthony Hammond

Examples

# Estimate the 50-year growth factor from an Lcv and Lskew of 0.17 and 0.04, respectively
GenLogGF(0.17, 0.04, RP = 50)

Generalised logistic distribution parameter estimates

Description

Estimated parameters from a sample (with Lmoments or maximum likelihood estimation) or from L1 (first L-moment), Lcv (linear coefficient of variation), and LSkew (linear skewness)

Usage

GenLogPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)

Arguments

Details

The L-moment estimated parameters are by the method detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

Parameter estimates (location, scale, shape)

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the parameters using L-moments
am_27090 <- GetAM(27090)
GenLogPars(am_27090$Flow)

# Estimate parameters using MLE
GenLogPars(am_27090$Flow, mle = TRUE)

# Calculate L-moments and estimate the parameters with L1, Lcv, and Lskew
LMoments(am_27090$Flow)

# Store L-moments in an object
l_pars <- as.numeric(LMoments(am_27090$Flow))[c(1, 5, 6)]
GenLogPars(L1 = l_pars[1], LCV = l_pars[2], LSKEW = l_pars[3])

Generalised Pareto distribution estimates from parameters

Description

Estimated quantiles as function of return period (RP) and vice versa, from user input parameters

Usage

GenParetoEst(loc, scale, shape, q = NULL, RP = 100, ppy = 1)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The average number of peaks per year argument (ppy) is necessary when ppy is not equal to 1.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get a POT sample, estimate the parameters, and estimate the 50-year RP
thames_pot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.90)
GenParetoPars(thames_pot$peak)

# Store the parameters in an object
pars <- as.numeric(GenParetoPars(thames_pot$peak))

# Get an estimate of 50-year flow
GenParetoEst(pars[1], pars[2], pars[3], ppy = 1.867, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GenParetoEst(pars[1], pars[2], pars[3], ppy = 1.867, q = 600)

Generalised Pareto distribution growth factors

Description

Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness). The Lcv and LSkew in this case should be calculated from peaks over threshold data and the ppy argument is necessary where the average number of peaks per year is not 1.

Usage

GenParetoGF(lcv, lskew, RP, ppy = 1)

Arguments

Details

Growth factors (GF) are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999. The average number of peaks per year argument (ppy) is for the function to convert from the peaks over threshold (POT) scale to the annual scale. For example, if there are 3 peaks per year, the probability associated with the 100-yr return period estimate would be 0.01/3 (i.e. an RP of 300 on the POT scale) rather than 0.01.

Value

Generalised Pareto estimated growth factor

Author(s)

Anthony Hammond

Examples

# Get POT flow data from the Thames at Kingston (noting the no. peaks per year).
# Then estimate the 100-year growth factor with lcv and lskew estimates
tpot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.90)
GenParetoGF(Lcv(tpot$peak), LSkew(tpot$peak), RP = 100, ppy = 1.867)

# Multiply by the median of the POT data for an estimate of the 100-year flood
GenParetoGF(Lcv(tpot$peak), LSkew(tpot$peak), RP = 100, ppy = 1.867) * median(tpot$peak)

Generalised Pareto distribution - estimates directly from sample

Description

Estimated quantiles as function of return period (RP) and vice versa, directly from the data

Usage

GenParetoPOT(x, ppy = 1, RP = 100, q = NULL)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The average number of peaks per year argument (ppy) is for the function to convert from the peaks over threshold (POT) scale to the annual scale. For example, if there are 3 peaks per year, the probability associated with the 100-yr return period estimate would be 0.01/3 (i.e. an RP of 300 on the POT scale) rather than 0.01. The parameters are estimated by the method of L-moments, as detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get a POT series and estimate the 50-year RP
thames_pot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.90)
GenParetoPOT(thames_pot$peak, ppy = 1.867, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GenParetoPOT(thames_pot$peak, ppy = 1.867, q = 600)

Generalised Pareto distribution parameter estimates

Description

Estimated parameters from a sample (with Lmoments or maximum likelihood estimation) or from L1 (first L-moment), Lcv (linear coefficient of variation), and LSkew (linear skewness)

Usage

GenParetoPars(x = NULL, L1, LCV, LSKEW)

Arguments

Details

The L-moment estimated parameters are by the method detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

Parameter estimates (location, scale, shape)

Author(s)

Anthony Hammond

Examples

# Get a peaks over threshold sample and estimate the parameters using L-moments
thames_pot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.90)
GenParetoPars(thames_pot$peak)

# Calculate L-moments and estimate the parameters with L1, Lcv, and Lskew
LMoments(thames_pot$peak)

# Store L-moments in an object
l_pars <- as.numeric(LMoments(thames_pot$peak))[c(1, 5, 6)]
GenParetoPars(L1 = l_pars[1], LCV = l_pars[2], LSKEW = l_pars[3])

Get an annual maximum sample from the National River Flow Archive sites suitable for pooling

Description

Extracts the annual maximum peak flow sample and associated dates for the site of interest.

Usage

GetAM(ref)

Arguments

Value

A data.frame with columns; Date, Flow, and id

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and display it in the console
GetAM(203018)

# Save an AMAX sample as an object
am_203018 <- GetAM(203018)

Get catchment descriptors from the National River Flow Archive sites considered suitable for median annual maximum flow estimation (QMED) and pooling.

Description

Extracts the catchment descriptors for a site of interest from the National River Flow Archive.

Usage

GetCDs(x)

Arguments

Details

If the site is considered suitable for QMED and pooling the CDs are extracted from the PeakFlowData data.frame. Otherwise they are extracted using the NRFA website. Note that if they are from the NRFA website then the 'easting' and 'northing' are not for the catchment centroid, they're for the gauge location. Also, where the gauge has NRFA peak flows available, but is not considered suitable for pooling or QMED, it will be derived from the NRFA webpage, and some descriptors differ a little between the data sets (NRFA website and NRFA peak flows), notably the Area.

Value

A data.frame with columns; Descriptor and Value.

Author(s)

Anthony Hammond

Examples

# Get CDs and display in the console
cds_203018 <- GetCDs(203018)
cds_203018

Get flow or level data from the Environment Agency's Hydrology Data Explorer

Description

Function to extract flow or level data from the Environment Agency's Hydrology Data Explorer.

Usage

GetDataEA_QH(
  Lat = 54,
  Lon = -2.25,
  Range = 20,
  RiverName = NULL,
  WISKI_ID = NULL,
  From = NULL,
  To = NULL,
  Type = "flow",
  Period = "DailyMean"
)

Arguments

Details

To find gauges you can input either a river name or a latitude and longitude. You can convert BNG to Lat and Lon using the ConvertGridRef function (you can also get lat and lon by left clicking on Google maps). The lat and lon option will provide all flow and level gauges within a specified range (default of 10km). This provides gauged details including the WISKI ID. You can get data from specific gauges using the WISKI_ID. Note that flow gauges also have level data available. You can get data from a date range using the From and To arguments or you can return all data by leaving the From and To as the default (NULL). Lastly, WISKI IDs are sometimes returned without a preceding 0 which might be necessary for the data extraction (oddly, most do have the necessary 0). If data extraction fails try adding a 0 to the beginning of the WISKI ID.

Value

If searching for gauge details with lat and lon or river name, then a list is returned. The first element is a dataframe with flow gauge details and the second is a dataframe of level gauge details. When extracting flow or level data with a WISKI ID then a dataframe with two columns is returned. The first being a Date or POSIXct column/vector and the second is the timeseries of interest.

Author(s)

Anthony Hammond

Examples

# Find gauges on the River Tame
## Not run: 
GetDataEA_QH(RiverName = "Tame")

## End(Not run)

# Find gauges within 10km of a latitude/longitude coordinate somewhere near the
# centre of Dartmoor
## Not run: 
GetDataEA_QH(Lat = 50.6, Lon = -3.9, Range = 10)

## End(Not run)

# Get all available daily maximum flow data from the Bellever gauge on the
# East Dart River
## Not run: 
bellever_max <- GetDataEA_QH(WISKI_ID = "SX67F051", Period = "DailyMax")

## End(Not run)

# Get 15-minute data from the Bellever for the November 2024 event
## Not run: 
bellever_nov_2024 <- GetDataEA_QH(
  WISKI_ID = "SX67F051",
  From = "2024-11-23", To = "2024-11-25", Period = "15Mins"
)

## End(Not run)

Get Environment Agency rainfall data (England).

Description

Extract rainfall data from the Environment Agency's Hydrology Data Explorer.

Usage

GetDataEA_Rain(
  Lat = 54,
  Lon = -2,
  Range = 10,
  WISKI_ID = NULL,
  Period = "Daily",
  From = NULL,
  To = NULL
)

Arguments

Details

The function provides one of two outputs. Either information about available local rain gauges, or the data from a specified gauge (specified by WISKI ID). The process is to find the local information (including WISKI ID) by using the latitude and longitude and range (You can convert BNG to Lat and Lon using the GridRefConvert function). Then use the WISKI ID to get the data. If data requested is not available, for example - outside the date range or not available at the requested sampling rate, an error message is returned stating "no lines available in input". To extract all the available data leave the From and To arguments as Null.

Value

If searching for rain gauge details with the Latitude and Longitude a dataframe of gauges is returned. If extracting rainfall using the WISKI_ID, a dataframe is returned with Date / POSIXct in the first columns and rainfall in the second.

Author(s)

Anthony Hammond

Examples

# Get information about available rain gauges
# within a 10km radius of latitude 54.5 and longitude -3.2
## Not run: 
GetDataEA_Rain(Lat = 54.5, Lon = -3.2)

## End(Not run)

# Use the WISKI reference ID for the Honister rain gauge
# to get some hourly rain data for December 2015
## Not run: 
honister_dec_2015 <- GetDataEA_Rain(
  WISKI_ID = "592463",
  Period = "Hourly", From = "2015-12-01", To = "2015-12-31"
)

## End(Not run)

# Inspect the first few rows and plot the data
## Not run: 
head(honister_dec_2015)
plot(honister_dec_2015, type = "h", ylab = "Rainfall (mm)")

## End(Not run)

Get regional Met Office average temperature or rainfall series (monthly, seasonal, and annual).

Description

Extracts regional mean temperature or rainfall from the Met Office UK & regional series. The total duration of bright sunshine is also available.

Usage

GetDataMetOffice(Variable, Region)

Arguments

Details

The function returns time series data from the 19th century through to the present month.

Value

A data.frame with 18 columns; year, months, seasons, and annual. Rows then represent each year of the timeseries.

Author(s)

Anthony Hammond

Examples

# Get the rainfall time series for the UK
## Not run: 
uk_rain <- GetDataMetOffice(Variable = "Rainfall", Region = "UK")

## End(Not run)

# Get the mean temperature data for East Anglia
## Not run: 
temp_east_anglia <- GetDataMetOffice(Variable = "Tmean", Region = "East_Anglia")

## End(Not run)

Get National River Flow Archive data using gauge ID.

Description

Extracts NRFA data using the API.

Usage

GetDataNRFA(ID = NULL, Type = "Q")

Arguments

Details

The function can be used to get daily catchment rainfall or mean flow, or both together (concurrent). It can also be used to get gaugings, AMAX, and POT data. Note that some sites have rejected peak flow years. In which case, if Type = AMAX or POT, the function returns a list, the first element of which is the rejected years, the second is the full AMAX or POT. Lastly if Type = "Catalogue" and ID is NULL, it will return a dataframe of all the NRFA gauges, associated details, comments, and descriptors. If Type equals "Catalogue" and a valid ID is used, then all these gauge details are provided for that gauge.

Value

A data.frame with date in the first columns and variable/s of interest in the remaining column/s. Except for the following circumstances: When Type = "Catalogue", then a large dataframe is returned with all the NRFA gauge metadata. When Type = "AMAX" or "POT" and there are rejected years a list is returned, where the first element is the dataframe of data and the second is rejected year/s (character string).

Author(s)

Anthony Hammond

Examples

# Get the concurrent rainfall (P) and mean flow (Q) series for the Tay at Ballathie (site 15006)
## Not run: 
ballathie_pq <- GetDataNRFA(15006, Type = "PQ")

## End(Not run)

# Get the gaugings
## Not run: 
ballathie_gaugings <- GetDataNRFA(15006, Type = "Gaugings")

## End(Not run)

Get Scottish Environment Protection Agency (SEPA) Flow, Level, or Rainfall data.

Description

Extract hydrometric data from SEPA's API.

Usage

GetDataSEPA(
  Lat = NULL,
  Lon = NULL,
  RiverName = NULL,
  Type = "Flow",
  StationID = NULL,
  Range = 20,
  From = NULL,
  To = NULL,
  Period = "Daily"
)

Arguments

Details

You can download data using the gauge station ID and you can find gauges within a given range using the latitude and longitude, or a river name. If the 'From' date is left as null, the earliest date of available data will be used. If the 'To' date is left as null, the most recent date of available data will be used.

Value

A data.frame with POSIXct in the first column, and rainfall in the second column. Unless the StationName provided is not in the available list, then the available list is returned.

Author(s)

Anthony Hammond

Examples

# search Rain gauges by Lat and Lon
## Not run: 
GetDataSEPA(Lat = 56, Lon = -4, Type = "Rain")

## End(Not run)

# Get daily rain from the Bannockburn station between two dates
## Not run: 
bannockburn <- GetDataSEPA(
  StationID = "36494",
  From = "1998-10-01", To = "1998-10-31"
)

## End(Not run)

# Inspect the first few rows and plot the data
## Not run: 
head(bannockburn)
plot(bannockburn, type = "h", ylab = "Rainfall (mm)")

## End(Not run)

QMED from a gauged site suitable for QMED

Description

Provides QMED (median annual maximum flow) from a site suitable for QMED, using the site reference. This provides the observed median of the annual maximum sample (excluding rejected observations) for the site of interest. It does not apply any adjustments or updates to account for non-stationarity.

Usage

GetQMED(x)

Arguments

Value

the median annual maximum

Author(s)

Anthony Hammond

Examples

# Get the observed QMED from site 55002
GetQMED(55002)

Goodness of fit comparison (single sample)

Description

compares the RMSE of four distribution fits for a single AMAX sample.

Usage

GoFCompare(x)

Arguments

Details

This function calculates an RMSE fit score for four distributions (GEV, GenLog, Gumbel, & Kappa3). The lowest RMSE is the best fit. It works as follows. For each distribution: Step1. Simulate 500 samples the same size as x. Step2. Calculate the mean across all 500 samples for each rank to create an ordered central estimate. Step3. Calculate the RMSE between the result of step 2 and the ordered x. Step4. Standardise the RMSE by dividing it by the mean of x and multiply it by 100 (RMSE as a percentage of mean). Note that this is not a hypothesis test. It is only for comparing the fit across the distributions.

Value

A list. The first element is a dataframe with four columns and one row of results. Each column has the standardised RMSE associated with one of the four distributions (GEV, GenLog, Gumbel, Kappa3). The second element is a character string stating the distribution with the best fit.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and then compare the fit
am_15006 <- GetAM(15006)
GoFCompare(am_15006$Flow)

Goodness of fit comparison (for a pooling group)

Description

compares the RMSE of four distribution fits for a pooling group.

Usage

GoFComparePool(x)

Arguments

Details

This function calculates an RMSE fit score for four distributions (GEV, GenLog, Gumbel, & Kappa3). The lowest RMSE is the best fit. It works for pooling groups created using the Pool or PoolSmall function. It uses the same method as GoFCompare (see the associated details of that function). It first standardises the pooled AMAX samples (by dividing them by median) and then treats them as a single large sample. Note that this is not a hypothesis test. It is only for comparing the fit across the distributions.

Value

A list. The first element is a dataframe with four columns and one row of results. Each column has the standardised RMSE associated with one of the four distributions (GEV, GenLog, Gumbel, Kappa3). The second element is a character string stating the distribution with the best fit.

Author(s)

Anthony Hammond

Examples

# Get a pooling group and then compare the fit
pool_60009 <- Pool(GetCDs(60009))
GoFComparePool(pool_60009)

Gumbel distribution - estimates directly from sample

Description

Estimated quantiles as a function of return period (RP) and vice versa, directly from the data

Usage

GumbelAM(x, RP = 100, q = NULL)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The parameters are estimated by the method of L-moments, as detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the 50-year RP
am_27090 <- GetAM(27090)
GumbelAM(am_27090$Flow, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GumbelAM(am_27090$Flow, q = 600)

Gumbel distribution estimates from parameters

Description

Estimated quantiles as function of return period (RP) and vice versa, from user input parameters

Usage

GumbelEst(loc, scale, q = NULL, RP = 100)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample, estimate the parameters, and estimate the 50-year RP
am_27090 <- GetAM(27090)
pars <- as.numeric(GumbelPars(am_27090$Flow))
GumbelEst(pars[1], pars[2], RP = 50)

# Estimate the RP for a 600 m^3/s discharge
GumbelEst(pars[1], pars[2], q = 600)

Gumbel distribution growth factors

Description

Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)

Usage

GumbelGF(lcv, RP)

Arguments

Details

Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.

Value

Gumbel estimated growth factor

Author(s)

Anthony Hammond

Examples

# Estimate the 50-year growth factor from an Lcv of 0.17
GumbelGF(0.17, RP = 50)

Gumbel distribution parameter estimates

Description

Estimated parameters from a sample (with Lmoments or maximum likelihood estimation) or from L1 (first L-moment), Lcv (linear coefficient of variation)

Usage

GumbelPars(x = NULL, mle = FALSE, L1, LCV)

Arguments

Details

The L-moment estimated parameters are by the method detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, New York'.

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

Parameter estimates (location, scale)

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the parameters using L-moments
am_27090 <- GetAM(27090)
GumbelPars(am_27090$Flow)

# Estimate parameters using MLE
GumbelPars(am_27090$Flow, mle = TRUE)

# Calculate L-moments and estimate the parameters with L1 and Lcv
pars <- as.numeric(LMoments(am_27090$Flow)[c(1, 5)])
GumbelPars(L1 = pars[1], LCV = pars[2])

Heterogeneity measure (H2) for pooling groups.

Description

Quantifies the heterogeneity of a pooled group

Usage

H2(x, H1 = FALSE)

Arguments

Details

The H2 measure was developed by Hosking & Wallis and can be found in their book 'Regional Frequency Analysis: An Approach Based on L-Moments' (1997). It was also adopted for use by the Flood Estimation Handbook (1999) and is described in volume 3. It works by recreating 500 pooling groups with the same sample sizes, assuming a four parameter Kappa distribution (parameters from the pooled L-moments). L-moment ratios are calculated for each of the 500 simulated pooling groups. The heterogeneity is determined by comparing the variance of L-moment ratios in the observed pooling group with the variance of the L-moment ratios across the simulated pooling groups. The simulations are homogeneous, therefore if the observed pooling group is homogeneous the expectation is that the variance will be similar to the average of the simulated variance.

Value

A vector of two characters; the first representing the H2 score and the second stating a qualitative measure of heterogeneity.

Author(s)

Anthony Hammond

Examples

# Get CDs, form a pooling group, and calculate H2
cds_203018 <- GetCDs(203018)
pool_203018 <- Pool(cds_203018)
H2(pool_203018)

Historic flood maximum likelihood estimation

Description

A function to estimate parameters from an annual maximum sample and a known number of historic floods.

Usage

HistoricMLE(x, k = NULL, h, threshold, dist = "GenLog", Uncertainty = FALSE)

Arguments

Details

This function applies the case where only the number of exceedances are known. Not the case where the discharge of the historic floods is known. This latter functionality will be added at a later date. Note that if Uncertainty is set to TRUE, a range of return periods and associated estimates are returned along with uncertainty - quantified as the FSE. In some cases the uncertainty can increase. This happens when the additional information (number of exceedances and time period) does not outweigh an increase to the scale of skew parameter. The uncertainty calculated is a function of sample size and variance.

Value

A data.frame with one column of flow estimates. The row names denote the name of each estimate. If Uncertainty is TRUE, a list is returned and the second element is a dataframe with estimates and factorial standard errors.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and assume 3 exceedances over 100 years
# with a threhsold of 140m 3/s
AM71011 <- GetAM(71011)
HistoricMLE(AM71011$Flow, k = 3, h = 100, threshold = 140)


# Now estimate again but set Uncertainty = TRUE
HistoricMLE(AM71011$Flow, k = 3, h = 100, threshold = 140, Uncertainty = TRUE)

Hydrological plot of concurrent discharge and precipitation

Description

Plots concurrent precipitation and discharge, with precipitation along the top and discharge along the bottom

Usage

HydroPlot(
  x,
  main = "Concurrent Rainfall & Discharge",
  ylab = "Discharge (m3/s)",
  From = NULL,
  To = NULL,
  adj.y = 1.5,
  RainAxisMax = NULL,
  plw = 1,
  qlw = 1.8,
  Return = FALSE
)

Arguments

Details

The input of x is a dataframe with the first column being time. If the data is sub daily this should be class POSIXct with time as well as date.

Value

A plot of concurrent precipitation and discharge, with the former at the top and the latter at the bottom. If the Return argument equals true the associated data-frame is also returned.

Author(s)

Anthony Hammond

Examples

# Plot the Thames precipitation and discharge for the 2013 hydrological year,
# adjusting the y axis to 1.8
HydroPlot(ThamesPQ, From = "2013-10-01", To = "2014-09-30", adj.y = 1.8)

Kappa3 distribution - estimates directly from sample

Description

Estimated quantiles as a function of return period (RP) and vice versa, directly from the data

Usage

Kappa3AM(x, RP = 100, q = NULL)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. The parameters are estimated by the method of L-moments, as detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'. This is the Kappa3 distribution as defined in Kjeldsen, T. (2019), 'The 3-parameter Kappa distribution as an alternative for use with FEH pooling groups.' (Circulation - The Newsletter of the British Hydrological Society, no. 142).

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the 50-year RP
am_27090 <- GetAM(27090)
Kappa3AM(am_27090$Flow, RP = 50)

# Estimate the RP for a 600 m^3/s discharge
Kappa3AM(am_27090$Flow, q = 600)

Kappa3 distribution estimates from parameters

Description

Estimated quantiles as function of return period (RP) and vice versa, from user input parameters

Usage

Kappa3Est(loc, scale, shape, q = NULL, RP = 100)

Arguments

Details

If the argument q is used, it overrides RP and provides RP as a function of q (magnitude of variable) as opposed to q as a function of RP. This is the Kappa3 distribution as defined in Kjeldsen, T (2019), 'The 3-parameter Kappa distribution as an alternative for use with FEH pooling groups.' (Circulation - The Newsletter of the British Hydrological Society, no. 142).

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

quantile as a function of RP or vice versa

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample, estimate the parameters, and the 50-year RP flow
am_27090 <- GetAM(27090)

# Get the parameters and store in an object
pars <- as.numeric(Kappa3Pars(am_27090$Flow))

# Get an estimate of the 50-year flow
Kappa3Est(pars[1], pars[2], pars[3], RP = 50)

# Estimate the RP for a 600 m^3/s discharge
Kappa3Est(pars[1], pars[2], pars[3], q = 600)

Kappa3 distribution growth factors

Description

Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)

Usage

Kappa3GF(lcv, lskew, RP)

Arguments

Details

Growth factors are calculated by the method outlined in Kjeldsen, T (2019), 'The 3-parameter Kappa distribution as an alternative for use with FEH pooling groups.'Circulation - The Newsletter of the British Hydrological Society, no. 142

Value

Kappa3 distribution estimated growth factor

Author(s)

Anthony Hammond

Examples

# Calculate Kappa growth factor for the 100-year flood
#assuming LCV and LSKEW of 0.165 and 0.17
Kappa3GF(0.165, 0.17, RP = 100)

Kappa3 distribution parameter estimates

Description

Estimated parameters from a sample (using Lmoments) or from user supplied L1 (first L-moment), Lcv (linear coefficient of variation), and LSkew (linear skewness)

Usage

Kappa3Pars(x = NULL, L1 = NULL, LCV = NULL, LSKEW = NULL)

Arguments

Details

The L-moment estimated parameters are by the method detailed in 'Hosking J. and Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'. This is the Kappa3 distribution as defined in Kjeldsen, T (2019), 'The 3-parameter Kappa distribution as an alternative for use with FEH pooling groups.' (Circulation - The Newsletter of the British Hydrological Society, no. 142).

This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

Parameter estimates (location, scale, shape)

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample and estimate the parameters
am_27090 <- GetAM(27090)
Kappa3Pars(am_27090$Flow)

# Calculate L-moments and estimate the parameters with L1, LCV, and LSKEW
l_pars <- as.numeric(LMoments(am_27090$Flow))[c(1, 5, 6)]
Kappa3Pars(L1 = l_pars[1], LCV = l_pars[2], LSKEW = l_pars[3])

Linear Kurtosis (LKurt)

Description

Calculates the LKurtosis from a sample of data

Usage

LKurt(x)

Arguments

Details

LKurtosis calculated according to methods outlined by Hosking & Wallis (1997): Regional Frequency Analysis and approach based on LMoments. Also in the Flood Estimation Handbook (1999), volume 3.

Value

Numeric. The LSkew of a sample.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and calculate the L-moments
am_27051 <- GetAM(27051)
LKurt(am_27051$Flow)

Lmoments & Lmoment ratios

Description

Calculates the Lmoments and Lmoment ratios from a sample of data

Usage

LMoments(x)

Arguments

Details

Lmoments calculated according to methods outlined by Hosking & Wallis (1997): Regional Frequency Analysis and approach based on LMoments. Also in the Flood Estimation Handbook (1999), volume 3.

Value

A data.frame with one row and column headings; L1, L2, L3, L4, Lcv, LSkew, and LKurt. The first four are the Lmoments and the next three are the Lmoment ratios.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and calculate the L-moments
am_27051 <- GetAM(27051)
LMoments(am_27051$Flow)

Adjust L-Ratios in a pooling group

Description

Adjusts the linear coefficient of variation (Lcv) and the linear skewness (LSkew) for a chosen site in a pooling group

Usage

LRatioChange(x, SiteID, lcv, lskew)

Arguments

Details

Pooling groups are formed from the PeakFlowData data.frame and all the Lcv and LSkew values are precalculated using the National River Flow Archive Peak flow dataset noted in the description file. The resulting pooled growth curve is calculated using the Lcv and Lskew in the pooled group. The user may have further data and be able to add further peak flows to the annual maximum samples within a pooling group. If that is the case a new Lcv and Lskew can be determined using the LMoments function. These new values can be added to the pooling group with this LRatioChange function. Also the non-flood years adjustment function may have been applied to a site, which provides a new Lcv and LSkew. In which case, the LRatioChange function can be applied. The function creates a new pooling group object and x will still exist in it's original state after the function is applied.

Value

A new pooling group, the same as x except for the user adjusted Lcv and Lskew for the user selected site.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors and create a pooling group
cds_39001 <- GetCDs(39001)
pool_39001 <- Pool(cds_39001, include = 39001)

# Apply the function to create a new adjusted pooling group,
# changing the subject site's lcv and lskew to 0.187 and 0.164, respectively
pool_39001_adj <- LRatioChange(pool_39001, SiteID = 39001, lcv = 0.187, lskew = 0.164)

Linear Skewness (LSkew)

Description

Calculates the LSkew from a sample of data

Usage

LSkew(x)

Arguments

Details

LSkew calculated according to methods outlined by Hosking & Wallis (1997): Regional Frequency Analysis and approach based on LMoments. Also in the Flood Estimation Handbook (1999), volume 3.

Value

Numeric. The LSkew of a sample.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and calculate the L-moments
am_27051 <- GetAM(27051)
LSkew(am_27051$Flow)

Linear coefficient of variation (Lcv)

Description

Calculates the Lcv from a sample of data

Usage

Lcv(x)

Arguments

Details

Lcv calculated according to methods outlined by Hosking & Wallis (1997): Regional Frequency Analysis and approach based on LMoments. Also in the Flood Estimation Handbook (1999), volume 3.

Value

Numeric. The Lcv of a sample.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and calculate the L-moments
am_27051 <- GetAM(27051)
Lcv(am_27051$Flow)

Urban adjustment for the linear coefficient of variation (Lcv)

Description

Urbanises or de-urbanises the Lcv using the FEH2025 methods

Usage

LcvUrb(LCV, URBEXT, DeUrb = FALSE)

Arguments

Value

The urban adjust Lcv or the de-urbanised Lcv

Author(s)

Anthony Hammond

Examples

# Apply a de-urbanisation with an LCV of 0.21 and an URBEXT of 0.1138
LcvUrb(0.21, 0.1138, DeUrb = TRUE)

# Apply and urban adjustment using LCV 0.196 and URBEXT of 0.1138
LcvUrb(0.196, 0.1138)

Low Flows

Description

A function to estimate lower flow quantiles in ungauged catchments.

Usage

LowFlows(CDs = NULL, AREA = NULL, SAAR = NULL, BFIHOST = NULL, Exclude = NULL)

Arguments

Details

This function provides estimates of the mean flow, Q95, Q70, Q50, Q10, and Q5. The function works by finding the 30 catchments in the NRFA data set with the most similar SAAR9120 to the subject site (via the API). The observed flows for those catchments are scaled by the catchment area. Then a weighted average is taken and multiplied by the subject site catchment area for the final estimate. The weighting is done by Eucidean distance based on SAAR9120 and BFIHOST19scaled. These are weighted based on the correlation of these descriptors to the scaled flows.

Value

A data.frame with one column of flow estimates. The row names denote the name of each estimate.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors, then estimate the flows
## Not run: 
CDs_27083 <- GetCDs(27083)
LowFlows(CDs_27083)

## End(Not run)
# Now estimate again but remove gauge 27083 from the analysis
## Not run: 
LowFlows(CDs_27083, Exclude = 27083)

## End(Not run)

Monthly Statistics

Description

Derives monthly statistics from a data.frame with Dates or POSIXct in the first column and variable of interest in the second

Usage

MonthlyStats(
  x,
  Stat,
  AggStat = NULL,
  TS = FALSE,
  Plot = FALSE,
  ylab = "Magnitude",
  main = "Monthly Statistics",
  col = "grey"
)

Arguments

Details

The statistic of interest for each month is calculated for each calendar year in the data.frame. An aggregated result is also calculated for each month using an aggregating statistic (the mean by default). The data.frame is first truncated at the first occurrence of January 1st and last occurrence of December 31st.

Value

A list with two elements. The first element is a data.frame with year in the first column and months in the next 12 (i.e. each row has the monthly stats for the year). The second element is a dataframe with month in the first column and the associated aggregated statistic in the second. i.e. the aggregated statistic (default is the mean) for each month is provided. However, of TS = TRUE, a monthly time series is returned - as a dataframe with date in the first column and monthly value in the second.

Author(s)

Anthony Hammond

Examples

# Get the mean flows for each month for the Thames at Kingston
qm_on_thames <- MonthlyStats(ThamesPQ[, c(1, 3)],
  Stat = mean,
  ylab = "Discharge (m^3/s)", main = "Thames at Kingston monthly mean flow", Plot = TRUE
)

# Get the monthly sums of rainfall for the Thames at Kingston
pm_on_thames <- MonthlyStats(ThamesPQ[, c(1, 2)],
  Stat = sum,
  ylab = "Rainfall (mm)", main = "Thames as Kingston monthly rainfall", Plot = TRUE
)

British national grid reference (NGR) distances

Description

Calculates the Euclidean distance between two British national grid reference points using the Pythagorean/Euclidean method.

Usage

NGRDist(i, j)

Arguments

Details

Note, that the result is converted to km from m.

Value

A distance in kilometres (if British national grid easting and northing are applied)

Author(s)

Anthony Hammond

Examples

# Calculate the distance between the catchment centroid for the
# Kingston upon Thames river gauge and the catchment centroid for the
# gauge at Ardlethen on the River Ythan.
# First retrieve the catchment descriptors (CDs) to obtain eastings and northings
GetCDs(10001)
GetCDs(39001)

# Calculate the distance between two centroids (eastings and northings)
NGRDist(i = c(381355, 839183), j = c(462899, 187850))

Non-flood adjustment

Description

Adjusts the linear coefficient of variation (Lcv) and the linear skewness (LSkew) to account for non-flood years

Usage

NonFloodAdj(x)

Arguments

Details

The method is the "permeable adjustment method" detailed in chapter 19, volume three of the Flood Estimation Handbook, 1999. The method makes no difference for sites where there are no annual maximums (AM) in the sample that are < median(AM)/2. Once applied the results can be used with the LRatioChange function to update the associated member of a pooling group. There is also the NonFloodAdjPool() function which can be used for multiple sites in a pooling group. The non flood adjustment procedure makes the assumption that annual maxima below QMED/2 are not from the same distribution and will result in a biased estimate. In turn it assumes that the AMAX are from a stationary process. The process adds uncertainty to the usual fitting process for three main reasons. Firstly, the definition of non-flood year (QMED/2). Secondly, the reduced sample size. Thirdly, the calculation process is based, in part, on the proportion of non-flood years to flood years. This proportion has uncertainty as a function of the sample size and the proportion because the standard error of a proportion (p) = sqrt((p * (1 - p)) / n).

Value

A list is returned. The first element of the list is a dataframe with one row and two columns - the adjusted Lcv in the first column and Lskew in the second. The second element of the list is another dataframe with one row and three columns. Number of non-flood years in the first column, sample size in the second and the percent of non-flood year in the third.

Author(s)

Anthony Hammond

Examples

# Get an annual maximum sample with a BFIHOST above 0.65 and with some
# annual maxima lower than half the median of the AMAX series, then apply the function
NonFloodAdj(GetAM(44013)[, 2])

Non-flood adjustment for pooling groups

Description

Applies the NonFloodAdj function to adjust the LCV and LSKEW of one or more sites in a pooling group.

Usage

NonFloodAdjPool(x, Index = NULL, AutoP = NULL, ReturnStats = FALSE)

Arguments

Details

For more details of the method for individual sites see the details section of the NonFloodAdj function. As a default this function applies NonFloodAdj to every member of the pooling group. Index can be supplied which is the row name/s of the members you wish to adjust. Or AutoP can be applied and is a percentage. Any member with a greater percentage of non-flood years than AutoP is then adjusted. The non-flood adjustment procedure makes the assumption that annual maxima below QMED/2 are not from the same distribution and will result in a biased estimate. In turn it assumes that the AMAX are from a stationary process. The process adds uncertainty to the usual fitting process for three main reasons. Firstly, the definition of non-flood year (QMED/2). Secondly, the reduced sample size. Thirdly, the calculation process is based, in part, on the proportion of non-flood years to flood years. This proportion has uncertainty as a function of the sample size and the proportion because the standard error of a proportion (p) = sqrt((p * (1 - p)) / n).

Value

By default the pooling group is returned with adjusted LCVs and LSKEWs for all sites indexed (or all sites when Index = NULL), or all sites with percentage of non-flood years above AutoP. No difference will be seen for sites with no AMAX < 0.5QMED. If ReturnStats is set to TRUE, a dataframe with Non-flood year stats is returned. The dataframe has a row for each site in the pooling group and three columns. The first is the number of non-flood years, the second is the number of years, and the third is the associated percentage.

Author(s)

Anthony Hammond

Examples

# Set up a pooling group for site 44013, then apply the function
pool_44013 <- Pool(GetCDs(44013), N = 500)
pool_nf <- NonFloodAdjPool(pool_44013)

# Return the non-flood stats for the pooling group
NonFloodAdjPool(pool_44013, ReturnStats = TRUE)

Optimise distribution parameters

Description

Estimates the parameters of the Generalised extreme value, generalised logistic, Kappa3, or Gumbel distribution from known return period estimates

Usage

OptimPars(x, dist = "GenLog")

Arguments

Details

Given a dataframe with return periods (RPs) in the first column and associated estimates in the second column, this function provides an estimate of the distribution parameters. Ideally the first RP should be 2. Extrapolation outside the RPs used for calibration comes with greater uncertainty.

Value

The parameters of one of four user chosen distributions; Generalised logistic, generalised extreme value, Gumbel, and Kappa3.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors and some quick results
# Then estimate the GenLog parameters
results <- QuickResults(GetCDs(27051))[[1]]
OptimPars(results[, 1:2])

Peaks over threshold (POT) data extraction

Description

Extracts independent peaks over a threshold from a sample

Usage

POTextract(
  x,
  div = NULL,
  TimeDiv = NULL,
  thresh = 0.975,
  Plot = TRUE,
  ylab = "Magnitude",
  xlab = "Time",
  main = "Peaks over threshold"
)

Arguments

Details

If the x argument is a numeric vector, the peaks will be extracted with no time information. x can instead be a data.frame with dates in the first column and the numeric vector in the second. In this latter case, the peaks will be time-stamped and a hydrograph, including POT, will be plotted by default. The method of extracting independent peaks assumes that there is a value either side of which events can be considered independent. For example, if two peaks above the chosen threshold are separated by the mean flow, they could be considered independent, but not if flow hasn't returned to the mean at any time between the peaks. Mean flow may not always be appropriate, in which case the 'div' argument can be applied (and is a percentile). The TimeDiv argument can also be applied to ensure the peaks are separated by a number of time-steps either side of the peaks. For extracting POT rainfall a div of zero could be used and TimeDiv can be used for further separation - which would be necessary for sub-daily time-series. Alternatively the POTt function can be applied (see associated details). When plotted, the blue line is the threshold, and the green line is the independence line (div).

Value

Prints the number of peaks per year and returns a data.frame with columns; Date and peak, with the option of a plot. Or a numeric vector of peaks is returned if only a numeric vector of the hydrological variable is input.

Author(s)

Anthony Hammond

Examples

# Extract POT data from Thames mean daily flow 2000-10-01 to 2015-09-30 with
# div = mean (default) and threshold = 0.95, and display the first six rows
thames_q_pot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.95)
head(thames_q_pot)

# Extract Thames POT from only the numeric vector of flows and display the
# first six rows
thames_q_pot <- POTextract(ThamesPQ[, 3], thresh = 0.9)
head(thames_q_pot)

# Extract the Thames POT precipitation with a div of 0, the default
# threshold, and five timesteps (days) either side of the peak, and display the first six rows
thames_p_pot <- POTextract(ThamesPQ[, c(1, 2)], div = 0, TimeDiv = 5)
head(thames_p_pot)

Peaks over threshold (POT) data extraction (quick)

Description

Extracts independent peaks over a threshold from a sample, using time as the independence criteria.

Usage

POTt(
  x,
  threshold = 0.975,
  div,
  Plot = TRUE,
  PlotType = "l",
  main = "Peaks over threhsold",
  ylab = "Magnitude",
  xlab = "Time"
)

Arguments

Details

This provides a quicker option than the POTextract function - useful for very long time series'. It only has the option of time division to ensure independence between peaks. If the x argument is a numeric vector, the peaks will be extracted with no time information. x can instead be a data.frame with dates in the first column and the numeric vector in the second. In this latter case, the peaks will be time-stamped and a hydrograph, including POT, will be plotted by default.

Value

A data.frame with columns; Date and peak, with the option of a plot. Or a numeric vector of peaks is returned if only a numeric vector of the variable is input as x.

Author(s)

Anthony Hammond

Examples

# Extract POT data from Thames catchment daily rainfall 2000-10-01 to 2015-09-30 with
# div = 14 (14 days) and threshold = 0.975, and display the first six rows
thames_p_pot <- POTt(ThamesPQ[, c(1, 2)], div = 14, threshold = 0.975)
head(thames_p_pot)

# Extract Thames rainfall POT from the numeric vector of rainfall, with threshold
# set to 0.95 and div set to 14, and display the first six rows
thames_p_pot <- POTt(ThamesPQ[, 2], threshold = 0.95, div = 14)
head(thames_p_pot)

National River Flow Archive descriptors and calculated statistics for sites suitable for QMED & pooling

Description

A data.frame of catchment & data descriptors. NRFA Peak Flow Dataset - Version 14

Usage

PeakFlowData

Format

A data frame with 902 rows and 36 variables

Details

Most of the FEH statistical method functions rely on this data. However, the data frame is open for manipulation in case the user wishes to add sites that aren't included, or change parts where local knowledge has improved on the data. If changes are made, they will only remain within the workspace. If a new workspace is opened and the UKFE package is loaded, the data frame will have returned to it's original state.

Source

https://nrfa.ceh.ac.uk/data/peak-flow-dataset

Create pooling group

Description

Function to develop a pooling group based on catchment descriptors

Usage

Pool(CDs, N = 800, UrbMax = 0.03, DeUrb = TRUE, exclude = NULL, include = NULL)

Arguments

Details

A pooling group is created from a CDs object, derived from GetCDs or CDsXML, or specifically with the catchment descriptors (see arguments). To change the default pooling group, one or more sites can be excluded using the 'exclude' option, which requires either a site reference or multiple site references in a vector. If this is done, the site with the next lowest similarity distance measure is added to the group (until the total number of years is at least N). Similarly a site can be included specifically by using the include argument. Sites with URBEXT2015 (urban extent) > 0.03 are excluded from the pooling group by default. This threshold can be adjusted with UrbMax. If DeUrb is set as TRUE (the default), the LCV values for sites in the pooling group are de-urbanised. If the user has more data available for a particular site within the pooling group, the LCV and LSKEW for the site can be updated after the group has been finalised using the LRatioChange function.

The pooling method is as specified by FEH2025. The de-urbanisation functionality assumes that the growth curve associated with an annual maximum flow sample is impacted by urbanisation and that this impact can be modelled as a function of the catchment URBEXT. The method for pooling the catchments together is based on the similarity of AREA, SAAR, FARL, FPEXT, and BFIHOST. These were seen to have the most significant impact on the LCV and LSKEW - and ultimately to provide the lowest 'Pooling Uncertainty Measure' (a statistic for assessing the similarity between pooled and single site gauged estimates).

Value

A data.frame of the pooling group with site reference row names and 24 columns, each providing catchment & gauge details for the sites in the pooling group.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors
cds_73005 <- GetCDs(73005)

# Set up a pooling group object called pool_73005 excluding sites 79005 & 46003
# Then print the group to the console
pool_73005 <- Pool(cds_73005, exclude = c(79005, 46003))
pool_73005

Pooled flood estimates

Description

Provides pooled results from a pooling group.

Usage

PoolEst(
  x,
  dist = "GenLog",
  CDs,
  QMEDEstimate,
  UrbAdj = TRUE,
  URBEXT = NULL,
  Gauged = FALSE,
  fseQMED = 1.5,
  Uncertainty = TRUE
)

Arguments

Details

PoolEst is a function to provide results from a pooling group derived using the Pool function. QMED (median annual maximum flow) needs to be supplied and can be derived from the QMED function for ungauged estimates or the annual maximum sample for gauged estimates. The method applied is based on FEH2025. The methods for estimating the L-moments and growth factors are outlined in the Flood Estimation Handbook (1999), volume 3. The estimation procedure assumes that the pooled AMAX samples are from the same underlying distribution (aside from the QMED scaling factor), that the distribution is correctly specified, that the individual samples are all independent and identically distributed, and that the samples are independent of each other. The urban adjustment (which is applied as default) assumes that the growth curve associated with an annual maximum flow sample is impacted by urbanisation and that this impact can be modelled as a function of the catchment URBEXT. A quantification of uncertainty is provided if Uncertainty is set to TRUE (the default). For more information see the Uncertainty function.

Value

A list of length 3. Element one is a data frame with columns; return period (a range from 2 - 1000), peak flow estimates (Q), growth factor estimates (GF), and factorial standard error (if Uncertainty = TRUE). The second element is the estimated Lcv and Lskew. The third provides distribution parameters for the frequency curve.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors and form a pooling group.
cds_27083 <- GetCDs(27083)
pool_27083 <- Pool(cds_27083)

#Get results assuming a GEV distribution
PoolEst(pool_27083, CDs = cds_27083, dist = "GEV", QMEDEstimate = 12, Uncertainty = FALSE)

QMED (median annual maximum flow) estimate from catchment descriptors

Description

Estimated median annual maximum flow from catchment descriptors and donor sites

Usage

QMED(
  CDs = NULL,
  DonorIDs = NULL,
  no.Donors = NULL,
  alpha = TRUE,
  UrbAdj = TRUE,
  UrbanExpansion = TRUE,
  UrbMax = 1,
  ReturnDetails = FALSE,
  AREA,
  SAAR,
  FARL,
  BFIHOST,
  URBEXT,
  Exclude = NULL
)

Arguments

Details

QMED is estimated from catchment descriptors: QMED = 6.8247*AREA^0.8499*0.1780^(1000/SAAR)*FARL^3.0450*0.0321^(BFIHOST^2) as specified by FEH2025. This QMED model is a multiple linear regression with transformed predictor variables and is trained on log transformed observed QMED values. The following assumptions are therefore applied: the relationship between the transformed independent variables (predictors) and the logarithmically transformed dependent variable (QMED) is linear, observations (observed QMEDs used in calibration) are independent of each other, model and sampling errors are independent of each other, model and sampling errors are normally distributed and have a mean of zero, predictor variables are independent, the cross-correlation of model errors can be described by distance between catchment centroids and the form of the associated correlation matrix is known prior to and for the calibration process. When UrbAdj = TRUE, urban adjustment is applied to the QMED estimate according to the FEH2025 method. The observed donor QMEDs are de-urbanised for the donor process and the donor adjusted rural estimate can then be urban adjusted accordingly (this is done by default). The use of the urban adjustment factor (UrbAdj) assumes that QMED is impacted by urbanisation and this impact can be determined by the URBEXT catchment descriptor. Use of the UrbanExpansion option applies a nationally averaged urban expansion factor to the URBEXT value, tending to overall underestimated urbanisation in more urban catchments and overestimated urbanisation in more rural catchments. Note that the distance-dependent moderation term (alpha) in the one-donor adjustment is not always appropriate, for example in some situations where the subject site is on the same watercourse as the donor. Similarly the two-donor distance-weighting method can give unsuitable results in some situations, for example where a subject site is in between the two donors on the same watercourse. Finally, for flexibility there is the option to input the relevant catchment descriptors directly rather than using a CDs object.

To derive an appropriate estimate when the donor catchment is urban ensure that DonUrbAdj is TRUE.

Value

An estimate of QMED from catchment descriptors. If donors are applied and ReturnDetails is TRUE, then an additional data frame is provided with associated details.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors and calculate QMED as if it was ungauged, with
# no donors, with one donor, two donors, and 8 donors.
cds_55004 <- GetCDs(55004)
QMED(cds_55004)
QMED(cds_55004, DonorIDs = 55012)
QMED(cds_55004, DonorIDs = c(55012, 60007))
QMED(cds_55004, no.Donors = 8)

QMED donor adjustment

Description

Applies a donor adjustment to the median annual maximum flow (QMED) estimate

Usage

QMEDDonEq(
  QMED.scd,
  QMEDgObs,
  QMEDgCds,
  Distance = NULL,
  xSI,
  ySI,
  xDon,
  yDon,
  alpha = TRUE
)

Arguments

Details

Although a single donor adjustment can be applied with the QMED function, this additional function is provided for flexibility. The method is that of FEH2025.

Author(s)

Anthony Hammond

Examples

# Get observed QMED for site 15006
q_ob <- GetQMED(15006)

# Get QMED equation estimated QMED for the donor site
q_cd <- QMED(CDs = GetCDs(15006))



# Apply the QMEDDonEq function with the information gained, assuming
# a distance of 30km and subject site QMED estimate of 3.9
QMEDDonEq(
  QMED.scd = 3.9, QMEDgObs = q_ob, QMEDgCds = q_cd,
  Distance = 30
)

QMED Linking equation

Description

Estimates the median annual maximum flow (QMED) from non-flood flows

Usage

QMEDLink(Q5dmf, Q10dmf, DPSBAR, BFI)

Arguments

Details

The QMED Linking equation estimates QMED as a function of the flow that is exceeded five percent of the time, the flow that is exceeded 10 percent of the time, the baseflow index, and the catchment descriptor drainage path slope (DPSBAR). All of these can be found for sites on the National River Flow Archive (NRFA) website. The method is provided in the guidance note 'WINFAP 4 QMED Linking equation' (2016) by Wallingford HydroSolutions.

Author(s)

Anthony Hammond

Examples

# Calculate the QMED for site 1001 (Wick at Tarroul)
QMEDLink(10.14, 7.352, 29.90, 0.39)

Empirical estimate of QMED from peaks over threshold (POT) data

Description

Estimates the median annual maximum flow (QMED) from peaks over threshold data

Usage

QMEDPOT(x, ppy)

Arguments

Details

If there are multiple peaks per year, the peaks per year (ppy) argument is used to convert to the annual scale to derive QMED. If ppy is one, then the median of the POT sample is returned (the median of x).

Author(s)

Anthony Hammond

Examples

# Extract some POT data and estimate QMED
thames_pot <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.90)
QMEDPOT(thames_pot$peak, ppy = 1.867263)

QMED factorial standard error for gauged sites

Description

Estimates the median annual maximum flow (QMED) factorial standard error (FSE) by bootstrapping the sample

Usage

QMEDfseSS(x)

Arguments

Details

The bootstrapping procedure resamples from the sample N*500 times with replacement. After splitting into 500 samples of size N, the median is calculated for each. Then the exponent of the standard deviation of the log transformed residuals is taken as the FSE. i.e. exp(sd(log(x)-mean(log(x)))), where x is the bootstrapped medians.

Value

The factorial standard error for the median of a sample.

Author(s)

Anthony Hammond

Examples

# Extract an AMAX sample and estimate the QMED factorial standard error
am_203018 <- GetAM(203018)
QMEDfseSS(am_203018$Flow)

Quick pooled results

Description

Provides pooled results, directly from the catchment descriptors

Usage

QuickResults(
  CDs,
  no.Donors = 8,
  dist = "GenLog",
  Qmed = NULL,
  UrbMax = 0.03,
  Include = NULL
)

Arguments

Details

The quick results function provides results with a default pooling group. Sites are chosen from those with URBEXT2015 below or equal to UrbMax (the default is 0.03).The LCVs in the pooling group are 'de-urbanised'. The final LCV estimate is then urban adjusted. QMED is estimated using the QMED function with eight donors, all of which have a de-urbanised observed QMED for the donor process. Then the QMED estimate has an urban adjustment applied. If the CDs are for a site suitable for pooling/QMED, this QMED estimate converges to the observed.

Value

A list of length three. Element one is a data frame with columns; return period (RP), peak flow estimates (Q) and growth factor estimates (GF). The second element is the estimated Lcv and Lskew (linear coefficient of variation and skewness). The third element is a dataframe with the distribution parameters.

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors
cds_73005 <- GetCDs(73005)

# Get results
QuickResults(cds_73005)

# Get results with a GEV distribution
QuickResults(cds_73005, dist = "GEV")

Stage-Discharge equation optimisation

Description

Optimises a power law rating equation from observed discharge and stage

Usage

Rating(x, a = NULL)

Arguments

Details

The power law rating equation optimised here has the form q = c(h+a)^n; where 'q' is flow, 'h' is the stage, c' and 'n' are constants, and 'a' is the stage when flow is zero. The optimisation uses all the data provided in the dataframe (x). If separate rating limbs are necessary, x can be subset per limb. i.e. the rating function would be used multiple times, once for each subset of x. There is the option, with the 'a' argument, to hold the stage correction parameter (a), at a user defined level. If 'a' is NULL it will be calibrated with 'c' & 'n' as part of the optimisation procedure. Note that this is a purely statistical procedure and hydraulic considerations may prove useful for improving results (particularly where extrapolation is required).

Value

A list with three elements. The first is a vector of the three calibrated rating parameters. The second is the rating equation; discharge as a function of stage. The third is the rating equation; stage as a function of discharge. A rating plot is also returned.

Author(s)

Anthony Hammond

Examples

# Create some dummy data
flow <- c(177.685, 240.898, 221.954, 205.55, 383.051, 154.061, 216.582)
stage <- c(1.855, 2.109, 2.037, 1.972, 2.574, 1.748, 2.016)
observations <- data.frame(flow, stage)

# Apply the rating function
Rating(observations)

# Apply the rating function with the stage correction at zero
Rating(observations, a = 0)

Revitalised Flood Hydrograph Model (ReFH)

Description

Provides outputs for the ReFH model from catchment descriptors or user defined inputs.

Usage

ReFH(
  CDs = NULL,
  Depth = NULL,
  Duration = NULL,
  Timestep = NULL,
  RainProfile = "FSR",
  PlotTitle = NULL,
  RPa = NULL,
  alpha = FALSE,
  WaterBalance = FALSE,
  Season = NULL,
  AREA = NULL,
  TP = NULL,
  BR = NULL,
  BL = NULL,
  Cmax = NULL,
  Cini = NULL,
  BFini = NULL,
  Rain = NULL,
  UHShape = "KT",
  UrbanLoss = FALSE,
  Loss = NULL,
  LossCini = NULL
)

Arguments

Details

As default this function is the ReFH model as described in the Flood Estimation Handbook Supplementary Report No.1 (2007). However, optional extras have been added such as an urban loss model and an option to ensure a water balance if it has been violated. The urban loss model is applied with default urban parameters (IF = 0.7, DS = 0.5, IRF = 0.4) and is described in the Wallingford Hydrosolutions report "ReFH2 Science Report Closing a Water Balance" (2019). The method to derive design rainfall profiles is described in the Flood Estimation Handbook (1999), volume 2. This has been slightly adjusted so that even profiles are possible. Users can also input their own rainfall with the 'Rain' argument. As a default, when catchment descriptors (CDs) are provided the ReFH function uses catchment descriptors to estimate the parameters of the ReFH model and an approximate two-year rainfall for the critical duration. If a parameter argument is used for one or more of the parameters, then these overwrite the CD derived parameters. This ReFH function is recommended for analysing the plausible catchment response to an input of rainfall. For this reason, and as noted, multiple additional features are available. The user can change components such as the unit hydrograph and the loss. The WaterBalance option can be applied to ensure it is not violated. Also, the baseflow component can be set as a constant (as opposed to a function of the runoff) by setting BR to zero. The rainfall can also be changed by choosing a range of different randomised profiles.

Value

A list with two elements, and a plot. First element of the list is a data.frame of parameters, initial conditions and the catchment area. The second is a data.frame with columns Rain, NetRain, Runoff, Baseflow, and TotalFlow. If the scale argument is used a numeric vector containing the scaled hydrograph is returned instead of the results dataframe. The plot is of the ReFH output, with rainfall, net-rainfall, baseflow, runoff and total flow. If the scaled argument is used, a scaled hydrograph is plotted.

Author(s)

Anthony Hammond

Examples

# Get CDs and apply the ReFH function
cds_203018 <- GetCDs(203018)
ReFH(cds_203018)

# Apply the ReFH function with a user defined initial baseflow
ReFH(cds_203018, BFini = 6)

Seasonal correction factor (SCF)

Description

The results of applying the ratio of the seasonal annual maximum rainfall for a given duration to the annual maximum rainfall for the same duration

Usage

SCF(SAAR, duration)

Arguments

Details

The SCF and its use is detailed in R&D Technical Report FD1913/TR - Revitalisation of the FSR/FEH rainfall runoff method (2005). The ReFH model has a design rainfall profile included for winter and summer but the depth duration frequency (DDF) model is calibrated on annual maximum peaks as opposed to seasonal peaks. The SCF is necessary to convert the DDF estimate to a seasonal one. Similarly, the DDF model is calibrated on point rainfall and the area reduction factor converts it to a catchment rainfall for use with a rainfall runoff model such as ReFH (see details of the ReFH function). The final depth, therefore, is; Depth = DDFdepth x ARF x SCF. Note that the SCF function (as detailed in FEH volume 2) was derived for durations of up to one day.

Value

A data.frame of one row and two columns: SCFSummer and SCFWinter.

Author(s)

Anthony Hammond

Examples

# Derive the SCF for a SAAR of 1981 and a duration of 6.5 hours
SCF(1981, 6.5)

Seasonality plot

Description

A plot to inspect the seasonality of peak flows

Usage

Seasonality(x, Lines = FALSE)

Arguments

Details

The dots (or dark lines if Lines = TRUE) show the season of individual peaks. The red line shows the average seasonality. The longer it is the more clustered in time the peaks are.

Value

A seasonality plot

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and plot the seasonality
am_27083 <- GetAM(27083)
Seasonality(am_27083)

# Now do the same with lines instead of dots
Seasonality(am_27083, Lines = TRUE)

Data simulator

Description

Simulation of a random sample from the generalised extreme value, generalised logistic, Gumbel, Kappa3, or generalised Pareto distributions

Usage

SimData(n, pars = NULL, dist = "GenLog", GF = NULL)

Arguments

Details

The simulated sample can be generated using the distribution parameters (pars) location, scale and shape, or the growth factor (GF) inputs linear coefficient of variation (Lcv), linear skewness (LSkew) & median annual maximum (QMED). This function applies a probability distribution model which assumes that the sample data is independent and identical, i.e. the assumption is that all observations in the sample would not impact or depend on any other. Furthermore, all observations are from the same underlying process which has not changed over the period of record (stationarity).

Value

A random sample of size n for the chosen distribution.

Author(s)

Anthony Hammond

Examples

# Simulate a sample of size 30 from a GenLog distribution with parameters 299, 51, -0.042
SimData(30, pars = c(299, 51, -0.042), dist = "GenLog")

# Now simulate using the Lcv, Lskew, and median (0.17, 0.04, 310)
SimData(30, GF = c(0.17, 0.04, 310), dist = "GenLog")

Kingston upon Thames daily flow and catchment precipitation 2000-10-01 to 2015-09-30

Description

A data.frame of four columns; Date, Precipitation (P), & daily mean flow (Q)

Usage

ThamesPQ

Format

A data frame with 5478 rows and 4 columns:

Date

Date

P

Precipitation, in mm

Q

Daily mean discharge, in m3/s

Source

https://nrfa.ceh.ac.uk/data/station/meanflow/39001

Trend hypothesis test

Description

A hypothesis test for trend

Usage

TrendTest(x, Variance = FALSE, method = "mk", alternative = "two.sided")

Arguments

Details

The test can be performed on a numeric vector, or a data.frame with dates in the first column and the associated variable of interest in the second. A choice can be made between Mann Kendall, Pearson, or Spearman tests. The Spearman and Mann Kendall are based on ranks and will therefore have the same results whether dates are included or not. The default is Mann Kendall. The default is to test for any trend (alternative = "two.sided"). For positive trend set alternative to "greater", and to test for negative trend set alternative to "less".

Interpretation: When testing for positive trend (alternative = "greater") the P_value is the probability of exceeding the observed statistic under the null hypothesis (that it is less than zero). The vice versa is true when testing for negative trend (alternative = "less"). For alternative = "two.sided" the P_value is the probability of exceeding the absolute value of the observed statistic under the null hypothesis (that it is different from zero). Low P values indicate that the null hypothesis is less likely.

Value

A data.frame with columns and associated values: P_value, statistic (Kendall's tau, Spearman's rho, or Pearson's correlation coefficient), and a standardised distribution value. The latter is either the z score (for MK test) or students 't' of the observed statistic under the null hypothesis.

Author(s)

Anthony Hammond

Examples

# Get an AMAX sample and apply a trend test with the default Mann-Kendall test
am_27083 <- GetAM(27083)
TrendTest(am_27083)

# Apply the test with the Pearson correlation method with dates
# included (full object) and not (flow values only)
TrendTest(am_27083, method = "pearson")
TrendTest(am_27083$Flow, method = "pearson")

# Apply the default Mann-Kendall test for positive trend
TrendTest(am_27083$Flow, alternative = "greater")

Urban adjustment factor (UAF)

Description

UAF from catchment descriptors for QMED estimation in ungauged urban catchments

Usage

UAF(CDs = NULL, URBEXT, BFIHOST, IF = 0.3, PRimp = 70)

Arguments

Details

The urban adjustment factor is to adjust the rural QMED estimates (as estimated using the QMED function) to urban estimates. This is necessary because the QMED equation is calibrated on rural catchments. The assumption is that the magnitude of QMED is impacted by urbanisation and that this impact can be modelled as a function of the catchment descriptors URBEXT and BFIHOST. This UAF function is based on URBEXT2015 and BFIHOST19scaled.

Value

The urban adjustment factor

Author(s)

Anthony Hammond

Examples

# Get some catchment descriptors for an urban catchment and calculate the UAF
cds_53006 <- GetCDs(53006)
UAF(cds_53006)

# Calculate UAF using a user input URBEXT2015 and BFIHOST19scaled
UAF(URBEXT = 0.1138, BFIHOST = 0.3620)

Urban expansion factor

Description

This function provides a coefficient to multiply by URBEXT2015 to adjust it to a given year

Usage

UEF(Year)

Arguments

Details

The urban expansion factor is that of the FEH2025 method. The urban expansion model assumes a national average expansion as a function of year. This means that on some catchments the value will be overestimated (primarily on rural ones) and on others the value will be underestimated (primarily on urban ones).

Value

A numeric urban expansion factor.

Author(s)

Anthony Hammond

Examples

# Get an expansion factor for the year 2025
UEF(2025)

UK outline

Description

Easting and northing national grid reference points around the coast of the UK

Usage

UKOutline

Format

A data frame with 3867 rows and 2 variables

X_BNG

Easting, British national grid reference

Y_BNG

Northing, British national grid reference

Source

https://environment.data.gov.uk/

Uncertainty quantification for gauged and ungauged pooled estimates

Description

Uncertainty for both the gauged and ungauged case are quantified specifically (bespoke) for the pooling group according to methods detailed in Hammond, A. (2021). Sampling uncertainty of UK design flood estimation. Hydrology Research. 1357-1371. 52 (6). Note that this function only quantifies sampling (aleatoric) uncertainty. It does not quantify uncertainty associated with models, model choices applied, or hydrometric data. The method assumes that AMAX samples within the pooling group are independent of each other and serially independent and identically distributed. The default ungauged QMED fse is 1.5. This is the FSE for the QMED 2025 regression model for all catchments suitable for QMED estimation. If Gauged = TRUE, the assumption is that the top site in the pooling group is the gauged subject site. In the Gauged case the pooling group is bootstrapped (parametrically) 200 times to create 200 pooling groups. Then the PoolEst function is applied to each of the 200 resampled pooling groups along with 200 resampled QMED from the gauged subject site. The FSE is calculated from the range of 200 estimates for each return period. In the ungauged case, each gauge is bootstrapped (parametrically) w x 200 times, where w is the weighting of the gauge. The growth factors for the associated 200 LMoment ratios are calculated and multiplied by a sampled QMED based on the fseQMED value. i.e SampledQMED = exp(rnorm(1, log(QMEDEstimate), log(fseQMED)). Then the FSEs are calculated across the 200 results for each return period.

Usage

Uncertainty(
  x,
  dist = "GenLog",
  Gauged = FALSE,
  QMEDEstimate = NULL,
  fseQMED = 1.5
)

Arguments

Value

A dataframe with 11 rows and two columns. Return period in the first column and factorial standard error in the second.

Author(s)

Anthony Hammond

Examples

# Derive a pooling group
pool_203018 <- Pool(GetCDs(203018))

# Calculate the factorial standard errors as if it it were ungauged.
Uncertainty(pool_203018, QMEDEstimate = QMED(GetCDs(203018)))

Weighted Lmoment ratios (LCV and LSKEW) from a pooling group

Description

Provides the weighted LCV and LSKEW for a pooling group

Usage

WeightedMoments(x)

Arguments

Details

Weighting method as according to: FEH2025

Value

A data.frame with site references in the first column and associated weights in the second

Author(s)

Anthony Hammond

Examples

# Get some CDs, form a gauged pooling group, and estimate gauged Lcv
cds_27051 <- GetCDs(27051)
pool_27051 <- Pool(cds_27051)
WeightedMoments(pool_27051)

Linear coefficient of variation (Lcv) weightings for a pooling group

Description

Provides the LCV weights for each site in a pooling group

Usage

WeightsLCV(x)

Arguments

Details

Weighting method for FEH2025

Value

A data.frame with site references in the first column and associated weights in the second

Author(s)

Anthony Hammond

Examples

# Get some CDs, form an ungauged pooling group, and estimate ungauged Lcv
cds_27051 <- GetCDs(27051)
pool_27051 <- Pool(cds_27051, exclude = 27051)
WeightsLCV(pool_27051)

Linear Skewness (LSKEW) weightings for a pooling group

Description

Provides the LSKEW weights for each site in a pooling group

Usage

WeightsLSKEW(x)

Arguments

Details

Weighting method for FEH2025.

Value

A data.frame with site references in the first column and associated weights in the second

Author(s)

Anthony Hammond

Examples

# Get some CDs, form an ungauged pooling group, and estimate ungauged LSkew
cds_27051 <- GetCDs(27051)
pool_27051 <- Pool(cds_27051, exclude = 27051)
WeightsLSKEW(pool_27051)

Zdist Goodness of fit measure for pooling groups

Description

Calculates the goodness of fit score for pooling groups.

Usage

Zdists(x)

Arguments

Details

The goodness of fit measure provides a Z-Score which quantifies the number of standard deviations from the mean of a normal distribution. To determine goodness of fit for a given distribution (assume GEV for this example), 500 pooling groups are formed which match the number of sites and samples sizes of the pooling group of interest. These are formed by simulation with the GEV distribution having LCV and LSKEW which are the weighted mean LCV and LSKEW of the pooling group (weighted by sample size) and a median of 1. The weighted mean L-Kurtosis of the observed pooling group (tR4) is compared to the mean and standard deviation (sd) of L-Kurtosis from the simulated pooling groups (tR4_Dist) by calculating the associated Z-score: (tR4 - mean(tR4_Dist)) / sd(tR4_Dist). The fit of the distribution can be considered acceptable if the absolute Z-Score is less than 1.645 (essentially a hypothesis test with alpha level equal to 0.1). This is done for all candidate distributions and the lowest absolute score is considered the best fit.

NOTE: This is slightly different from the zdist function described in the science report 'Improving the FEH statistical procedures for flood frequency estimation, Environment Agency (2008)'. That function assumes a theoretical LKurtosis as a function of the pooled LSKEW to compare with a distribution of LKurtosis from simulated pooling groups. This means that the Gumbel distribution cannot be compared (hence the change which is a recommendation in 'Regional Frequency Analysis' by Hosking & Wallis (1997)), i.e. the Gumbel distribution is now included whereas it previously could not be.

Value

A list with the first element a data.frame of four Z-Scores related to the columns; "GEV", "GenLog", "Gumbel", and "Kappa3". The second element is a character stating which has the best fit.

Author(s)

Anthony Hammond

Examples

# Get CDs, form a pooling group, and calculate the Z-dists
cds_39001 <- GetCDs(39001)
pool_39001 <- Pool(cds_39001, N = 500)
Zdists(pool_39001)