| Title: | UK Flood Estimation |
|---|---|
| Description: | Functions to implement the methods of the Flood Estimation Handbook (FEH), associated updates and the revitalised flood hydrograph model (ReFH). Currently the package uses NRFA peak flow dataset version 13. Aside from FEH functionality, further hydrological functions are available. Most of the methods implemented in this package are described in one or more of the following: "Flood Estimation Handbook", Centre for Ecology & Hydrology (1999, ISBN:0 948540 94 X). "Flood Estimation Handbook Supplementary Report No. 1", Kjeldsen (2007, ISBN:0 903741 15 7). "Regional Frequency Analysis - an approach based on L-moments", Hosking & Wallis (1997, ISBN: 978 0 521 01940 8). "Proposal of the extreme rank plot for extreme value analysis: with an emphasis on flood frequency studies", Hammond (2019, <doi:10.2166/nh.2019.157>). "Making better use of local data in flood frequency estimation", Environment Agency (2017, ISBN: 978 1 84911 387 8). "Sampling uncertainty of UK design flood estimation" , Hammond (2021, <doi:10.2166/nh.2021.059>). "Improving the FEH statistical procedures for flood frequency estimation", Environment Agency (2008, ISBN: 978 1 84432 920 5). "Low flow estimation in the United Kingdom", Institute of Hydrology (1992, ISBN 0 948540 45 1). Wallingford HydroSolutions, (2016, <http://software.hydrosolutions.co.uk/winfap4/Urban-Adjustment-Procedure-Technical-Note.pdf>). Data from the UK National River Flow Archive (<https://nrfa.ceh.ac.uk/>, terms and conditions: <https://nrfa.ceh.ac.uk/costs-terms-and-conditions>). |
| Authors: | Anthony Hammond [aut, cre] |
| Maintainer: | Anthony Hammond <[email protected]> |
| License: | GPL-3 |
| Version: | 0.3.9 |
| Built: | 2025-01-09 20:32:24 UTC |
| Source: | CRAN |
This function allows the user to add an AMAX sample and associated catchment descriptors for use with the FEH process.
AddGauge(CDs, AMAX, ID)AddGauge(CDs, AMAX, ID)
CDs |
catchment descriptor object imported using the CDsXML function. |
AMAX |
Either a data.frame with date (or POSIXct) in the first column and a numeric vector in the second (the AMAX). Or an AMAX sample (a numeric vector). |
ID |
This is a user supplied identification number for the AMAX. |
The function provides the necessary AMAX sample statistics and data.frame for adding catchment descriptors to the NRFAData 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.
A list object. The first element is a data.frame which is a row of statistics and descriptors to be added to the NRFAData data.frame. The second element is the AMAX sample formatted to be added to the AMSP data.frame
Anthony Hammond
#Read in AMAX and catchment descriptors ## Not run: AMAdd <- AMImport(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.am}") ## Not run: CDsAdd <- CDsXML(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.xml}") #Apply the function and then add the results to the necessary data.frames. ## Not run: Gauge27003 <- AddGauge(CDsAdd, AMAdd, ID = "27003") #Add the descriptors and stats (the first element of the output) to the NRFAData data.frame ## Not run: NRFAData <- rbind(NRFAData, Gauge27003[[1]]) #Add the AMAX to the AMSP data.frame. ## Not run: AMSP <- rbind(AMSP, Gauge27003[[2]])#Read in AMAX and catchment descriptors ## Not run: AMAdd <- AMImport(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.am}") ## Not run: CDsAdd <- CDsXML(r"{D:\NRFAPeakFlow_v12_1_0\suitable-for-neither\027003.xml}") #Apply the function and then add the results to the necessary data.frames. ## Not run: Gauge27003 <- AddGauge(CDsAdd, AMAdd, ID = "27003") #Add the descriptors and stats (the first element of the output) to the NRFAData data.frame ## Not run: NRFAData <- rbind(NRFAData, Gauge27003[[1]]) #Add the AMAX to the AMSP data.frame. ## Not run: AMSP <- rbind(AMSP, Gauge27003[[2]])
Aggregates time series data, creating hourly data from 15 minute data for example.
AggDayHour(x, func, Freq = "Day", hour = 9)AggDayHour(x, func, Freq = "Day", hour = 9)
x |
a data.frame with POSIXct in the first column and numeric vector in the second. |
func |
the function used for aggregation; mean, max, or sum, for example. |
Freq |
Choices are "Day", or "Hour". |
hour |
An integer between 0 and 23. This is used if "Day" is chosen in the Freq argument to determine when the day starts. |
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.
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
Anthony Hammond
#Create a data.frame with a normally distributed variable at #a 15 minute sampling rate. TS <- 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 <- data.frame(DateTime = TS, Var = rnorm(length(TS), 10, 2)) #use the function to aggregate to an hourly sampling rate, taking the maximum of each hour Hourly <- AggDayHour(TS, func = max, Freq = "Hour") #now aggregate with the mean at a daily scale Daily <- AggDayHour(TS, func = mean, Freq = "Day")#Create a data.frame with a normally distributed variable at #a 15 minute sampling rate. TS <- 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 <- data.frame(DateTime = TS, Var = rnorm(length(TS), 10, 2)) #use the function to aggregate to an hourly sampling rate, taking the maximum of each hour Hourly <- AggDayHour(TS, func = max, Freq = "Hour") #now aggregate with the mean at a daily scale Daily <- AggDayHour(TS, func = mean, Freq = "Day")
Imports the peak flows and dates from from NRFA peak flow .AM files, exluding the rejected years
AMImport(x)AMImport(x)
x |
the file path for the .AM file |
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.
A data.frame with columns; Date and Flow
Anthony Hammond
#Import an AMAX sample and display the first six rows in the console ## Not run: AM.4003 <- AMImport("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/4003.AM") ## Not run: head(AM.4003)]#Import an AMAX sample and display the first six rows in the console ## Not run: AM.4003 <- AMImport("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/4003.AM") ## Not run: head(AM.4003)]
Provides a barplot for an annual maximum sample
AMplot(x, ylab = "Discharge (m3/s)", xlab = "Hydrological year", main = NULL)AMplot(x, ylab = "Discharge (m3/s)", xlab = "Hydrological year", main = NULL)
x |
A data.frame with at least two columns. The first a date column and the second the annual maximum (AM) sequence. A third column with the station id can be applied which is then used for the default plot title. |
ylab |
Label for the y axis (character string). |
xlab |
Label for the x axis (character string). |
main |
Title for the plot (character string). The default is 'Annual maximum sample:', where : is followed by an ID number if this is included in a third column of the dataframe x. |
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.
A barplot of the annual maximum sample
Anthony Hammond
#Get an AMAX sample and plot AMplot(GetAM(58002))#Get an AMAX sample and plot AMplot(GetAM(58002))
A data.frame with three columns; Date, Flow, id. NRFA Peak Flow Dataset - Version 13.0.2
AMSPAMSP
A data frame with 26539 rows and 3 columns
Date
Annual maximum peak flow, m3/s
Station identification number
https://nrfa.ceh.ac.uk/peak-flow-dataset
Extracts annual statistics (default maximums) from a data.frame which has dates (or POSIXct) in the first column and variable in the second.
AnnualStat( x, Stat = max, Truncate = TRUE, Mon = 10, Hr = 9, Sliding = FALSE, N = 24, ... )AnnualStat( x, Stat = max, Truncate = TRUE, Mon = 10, Hr = 9, Sliding = FALSE, N = 24, ... )
x |
a data.frame with dates (or POSIXct) in the first column and variable in the second |
Stat |
A user chosen function to extract statistics, for example mean. The default it max. User supplied functions could also be used. |
Truncate |
Logical argument with a default of TRUE. If TRUE, then x is truncated to be within the first and last occurrence of the chosen month and time. If FALSE tuncation is not done and results from partial years will be included. |
Mon |
Choice of month as a numeric, from 1 to 12. The default is 10 which means the year starts October 1st. |
Hr |
Choice of hour to start the year (numeric from 0 to 23). The default is 9. |
Sliding |
Logical argument with a default of FALSE. This can be applied if you want the statistic over a sliding period. For example, deriving maximum annual rainfall totals over a 24 hour period, rather than the maximum daily totals. The number of periods (timesteps) is chosen with the N argument. If for example you want the annual maximum sum of rainfall over a 24 hour period, and you have 15minute data, the Stat input would be sum, and N would be 96 (because there are 96 15 minute periods in 24 hours). |
N |
Number of timesteps to slide over - used in conjunction with Sliding. The default is 24, make sure to adjust this depending on the duration of interest and the sampling rate of the input data. |
... |
further arguments for the stat function. Such as na.rm = TRUE. |
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.
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.
Anthony Hammond
#Extract the Thames AMAX daily mean flow and display the first six rows ThamesAM <- AnnualStat(ThamesPQ[,c(1,3)]) head(ThamesAM) #Extract the annual rainfall totals. ThamesAnnualP <- AnnualStat(ThamesPQ[,1:2], Stat = sum) #Extract maximum five day rainfall totals from the Thames rain ThamesP5DayAM <- AnnualStat(ThamesPQ[,1:2], Stat = sum, Sliding = TRUE, N = 5)#Extract the Thames AMAX daily mean flow and display the first six rows ThamesAM <- AnnualStat(ThamesPQ[,c(1,3)]) head(ThamesAM) #Extract the annual rainfall totals. ThamesAnnualP <- AnnualStat(ThamesPQ[,1:2], Stat = sum) #Extract maximum five day rainfall totals from the Thames rain ThamesP5DayAM <- AnnualStat(ThamesPQ[,1:2], Stat = sum, Sliding = TRUE, N = 5)
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.
ARF(Depth, Area, D)ARF(Depth, Area, D)
Depth |
depth of rainfall |
Area |
catchment area in km2 |
D |
duration in hours |
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.
the rainfall depth or rainfall return period
Anthony Hammond
#Derive the ARF for a depth of 30, an area of 500km2 and a duration of 12 hours ARF(30, 500, 12)#Derive the ARF for a depth of 30, an area of 500km2 and a duration of 12 hours ARF(30, 500, 12)
Calculates the baseflow index from a daily mean flow series
BFI(Q, x.lim = NULL, y.lim = NULL, PlotTitle = "Baseflow plot", Plot = TRUE)BFI(Q, x.lim = NULL, y.lim = NULL, PlotTitle = "Baseflow plot", Plot = TRUE)
Q |
the daily mean flow series. Numeric vector |
x.lim |
the x axis limits of the plot. Numeric vector of length two Default is the extents of the data |
y.lim |
the y axis limits of the plot. Numeric vector of length two. Default is the extents of the data |
PlotTitle |
the title of the plot. The default is "Baseflow plot" |
Plot |
a logical argument with a default of TRUE. If TRUE the daily flow is plotted with the baseflow highlighted. |
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)
the baseflow index and if Plot equals TRUE, a plot showing the flow time series (black) and the associated baseflow (red)
Anthony Hammond
# Calculate the BFI from daily discharge at Kingston upon Thames; # which is in column three of the ThamesPQ data BFI(ThamesPQ[,3])# Calculate the BFI from daily discharge at Kingston upon Thames; # which is in column three of the ThamesPQ data BFI(ThamesPQ[,3])
Resampling with replacement to approximate the sampling distribution of a statistic and quantify uncertainty.
Bootstrap(x, Stat, n = 500, Conf = 0.95, ReturnSD = FALSE, ...)Bootstrap(x, Stat, n = 500, Conf = 0.95, ReturnSD = FALSE, ...)
x |
a numeric vector. The sample of interest |
Stat |
the function (to calculate the statistic) to be applied to the bootstrapped samples. For example mean, max, or median. |
n |
the number of boostrapped samples (default 500). i.e. the size of the derived sampling distribution. |
Conf |
the confidence level of the intervals (default 0.95). Must be between 0 and 1. |
ReturnSD |
Logical argument with a default of FALSE. If true the bootstrapped sampling distribution is returned. |
... |
further arguments for the Stat function. For example if you use the GEVAM function you might want to add RP = 50 to derive a sampling distribution for the 50-year quantile. |
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.
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.
Anthony Hammond
#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 percent confidence level Bootstrap(AM.203018$Flow, Stat = sd, Conf = 0.90) #Return the sampling distribution of the mean and plot an associated histogram SampDist <- Bootstrap(AM.203018$Flow, Stat = mean, ReturnSD = TRUE) hist(SampDist)#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 percent confidence level Bootstrap(AM.203018$Flow, Stat = sd, Conf = 0.90) #Return the sampling distribution of the mean and plot an associated histogram SampDist <- Bootstrap(AM.203018$Flow, Stat = mean, ReturnSD = TRUE) hist(SampDist)
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
CDsXML(x)CDsXML(x)
x |
the xml file path |
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.
A data.frame with columns; Descriptor and Value.
Anthony Hammond
#Import catchment descriptors from a NRFA peakflows xml file and display in console ## Not run: CDs.4003 <- CDsXML("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/4003.xml") ## Not run: CDs.4003 #Import catchment descriptors from a FEH webserver xml file and display xml in the console ## Not run: CDs.MySite <- CDsXML("C:/Data/FEH_Catchment_384200_458200.xml")#Import catchment descriptors from a NRFA peakflows xml file and display in console ## Not run: CDs.4003 <- CDsXML("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/4003.xml") ## Not run: CDs.4003 #Import catchment descriptors from a FEH webserver xml file and display xml in the console ## Not run: CDs.MySite <- CDsXML("C:/Data/FEH_Catchment_384200_458200.xml")
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)
ConvertGridRef(x, fromBNG = TRUE, IGorLatLon = "LatLon")ConvertGridRef(x, fromBNG = TRUE, IGorLatLon = "LatLon")
x |
A vector of length 2. Either latitude and longitude (if fromBNG = FALSE) or BNG easting and northing (if fromBNG = TRUE). Or Irish easting and northing if IGorLatLon is set to IG and fromBNG = FALSE. |
fromBNG |
A logical argument with a default of TRUE. When TRUE it converts from BNG easting and northing to latitude and longitude (or to IG easting and northing if IGorLatLon is set to "IG"). When FALSE it converts the other way round. |
IGorLatLon |
This argument allows you to choose between Latitude & Longitude and Irish grid reference. The acceptable options are "LatLon" or "IG". If you choose "IG" you are converting between BNG and IG. If you choose "LatLon", you are converting between BNG and Lat Lon. |
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.
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".
Anthony Hammond
#Get Latitude and Longitude for a BNG numeric reference. ConvertGridRef(c(462899, 187850)) #Now we'll get easting and northing as a function of latitude and longitude ConvertGridRef(c(51.6, -1), fromBNG = FALSE)#Get Latitude and Longitude for a BNG numeric reference. ConvertGridRef(c(462899, 187850)) #Now we'll get easting and northing as a function of latitude and longitude ConvertGridRef(c(51.6, -1), fromBNG = FALSE)
Extracts results from a data frame imported using the DDFImport function
DDF(x, duration, RP = 100)DDF(x, duration, RP = 100)
x |
A data frame of DDF13 or DDF22 results imported using the DDFImport function |
duration |
the duration (hrs) for which a rainfall depth estimate is required |
RP |
the return period (years) for which a rainfall depth estimate is required |
The .xml files only provide a set number of durations and return periods for DDF13 and DDF22. This function optimises the GEV distribution on the results in order to interpolate across return periods. A linear interpolation is used between durations. The interpolation method may provide results that differ from the FEH webserver in the region of 0.1mm. The result is then rounded to an integer.
A DDF13 or DDF22 estimate of rainfall depth (mm)
Anthony Hammond
#Get DDF13 results from a the DDF ## Not run: DDF13.4003 <- DDFImport("C:/Data/NRFAPeakFlow_v9/Suitable for QMED/04003.xml") #Estimate the 20-year, 5 hour depth ## Not run: DDF(DDF13.4003, duration = 5, RP = 20)#Get DDF13 results from a the DDF ## Not run: DDF13.4003 <- DDFImport("C:/Data/NRFAPeakFlow_v9/Suitable for QMED/04003.xml") #Estimate the 20-year, 5 hour depth ## Not run: DDF(DDF13.4003, duration = 5, RP = 20)
Estimation of design rainfall depths, and the rarity of observed rainfall
DDF99(Duration, RP, pars, Depth = NULL, disc = NULL)DDF99(Duration, RP, pars, Depth = NULL, disc = NULL)
Duration |
numeric. The duration of interest (in hours) |
RP |
return period |
pars |
a numeric vector of length six. The six catchment parameters for the DDF model in the order of: c, d1, d2, d3, e, f |
Depth |
a user supplied rainfall depth for the duration under question |
disc |
converts from the sliding duration to fixed duration estimate. Choices are "hourly" or "daily" |
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.
the rainfall depth or rainfall return period
Anthony Hammond
#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))#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))
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
DDF99Pars(x)DDF99Pars(x)
x |
the xml file path |
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.
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.
Anthony Hammond
#Import DDF99 parameters from a NRFA peakflows xml file and display in console ## Not run: DDF99.4003 <- DDF99Pars("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/04003.xml") ## Not run: DDF99.4003 #Import DDF99 parameters from a FEH webserver xml file and display in the console ## Not run: DDF99.MySite <- DDF99Pars("C:/Data/FEH_Catchment_384200_458200.xml")#Import DDF99 parameters from a NRFA peakflows xml file and display in console ## Not run: DDF99.4003 <- DDF99Pars("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/04003.xml") ## Not run: DDF99.4003 #Import DDF99 parameters from a FEH webserver xml file and display in the console ## Not run: DDF99.MySite <- DDF99Pars("C:/Data/FEH_Catchment_384200_458200.xml")
Derive and plot rainfall Depth Duration Frequency curves from an input of rainfall data.
DDFExtract(x, Plot = TRUE, main = NULL)DDFExtract(x, Plot = TRUE, main = NULL)
x |
A data.frame with POSIXct in the first column and rainfall in the second. The data must have an hourly or sub-hourly sampling rate. |
Plot |
Logical argument with a default of TRUE. If TRUE, the DDF curves are plotted. |
main |
Title for the plot (character string). The default is no title. |
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.
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.
Anthony Hammond
# We'll extract all available 15 minute rainfall from the St Ives (Cambridgeshire) #rain gauge (WISKI ID = 179365). ## Not run: StIves <- GetDataEA_Rain(WISKI_ID = "179365", Period = "15Mins") #Then apply the DDF function. ## Not run: DDFExtract(StIves)# We'll extract all available 15 minute rainfall from the St Ives (Cambridgeshire) #rain gauge (WISKI ID = 179365). ## Not run: StIves <- GetDataEA_Rain(WISKI_ID = "179365", Period = "15Mins") #Then apply the DDF function. ## Not run: DDFExtract(StIves)
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
DDFImport(x, ARF = FALSE, Plot = TRUE, DDFVersion = 22)DDFImport(x, ARF = FALSE, Plot = TRUE, DDFVersion = 22)
x |
the xml file path |
ARF |
logical argument with a default of FALSE. If TRUE, the areal reduction factor is applied to the results. If FALSE, no area reduction factor is applied |
Plot |
logical argument with a default of TRUE. If TRUE the DDF curve is plotted for a few return periods |
DDFVersion |
Version of the DDF model (numeric). either 22 or 13. The default is 22. |
This function returns a data-frame of results. For further durations and return periods the separate DDF function can be applied with the data-frame as the argument/input. 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.
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.
Anthony Hammond
#Import DDF22 results from a NRFA peakflows xml file and display in console ## Not run: DDF22.4003 <- DDFImport("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/04003.xml") ## Not run: DDF22.4003 #Import DDF22 results from a FEH webserver xml file and display in the console ## Not run: DDF22.MySite <- DDFImport("C:/Data/FEH_Catchment_384200_458200.xml")#Import DDF22 results from a NRFA peakflows xml file and display in console ## Not run: DDF22.4003 <- DDFImport("C:/Data/NRFAPeakFlow_v11/Suitable for QMED/04003.xml") ## Not run: DDF22.4003 #Import DDF22 results from a FEH webserver xml file and display in the console ## Not run: DDF22.MySite <- DDFImport("C:/Data/FEH_Catchment_384200_458200.xml")
Extracts a mean hydrograph from a flow series
DesHydro( x, Threshold = 0.975, EventSep, N = 10, Exclude = NULL, Plot = TRUE, main = "Design Hydrograph" )DesHydro( x, Threshold = 0.975, EventSep, N = 10, Exclude = NULL, Plot = TRUE, main = "Design Hydrograph" )
x |
a dataframe with Date or POSIXct in the first coumn and the numeric vector of discharge in the second |
Threshold |
The threshold above which the peaks of the hydrograph are first identified. The default is 0.975. |
EventSep |
Number of timesteps to determine individual peak events during the extraction process. For the comparison and averaging process the start and end point of the hydrograph is Peak - EventSep and Peak + EventSep * 1.5. |
N |
number of event hydrographs from which to derive the mean hydrograph. Default is 10. Depending on the length of x, there may be fewer than 10 |
Exclude |
An index (single integer or vector of integers up to N) for which hydrographs to exclude if you so wish. This may require some trial and error. You may want to increase N for every excluded hydrograph. |
Plot |
logical argument with a default of TRUE. If TRUE, all the hydrographs from which the mean is derived are plotted along with the mean hydrograph. |
main |
Title for the plot |
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 EventsSep. 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.
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
Anthony Hammond
#Extract a design hydrograph from the Thames daily mean flow. Then print the resulting hydrograph ThamesDesHydro <- DesHydro(ThamesPQ[,c(1,3)], EventSep = 10, N = 10)#Extract a design hydrograph from the Thames daily mean flow. Then print the resulting hydrograph ThamesDesHydro <- DesHydro(ThamesPQ[,c(1,3)], EventSep = 10, N = 10)
Applies a linear detrend to a sample
DeTrend(x)DeTrend(x)
x |
a numeric vector |
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.
A numeric vector which is a linearly detrended version of x.
Anthony Hammond
# 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 (m3/s)") AM.Detrend <- DeTrend(AM.21025$Flow) plot(AM.Detrend, type = "h", ylab = "Discharge (m3/s)")# 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 (m3/s)") AM.Detrend <- DeTrend(AM.21025$Flow) plot(AM.Detrend, type = "h", ylab = "Discharge (m3/s)")
Provides 10 plots to compare the sites in the pooling group
DiagPlots(x, gauged = FALSE)DiagPlots(x, gauged = FALSE)
x |
pooling group derived from the Pool() function |
gauged |
logical argument with a default of FALSE. TRUE adds the top site in the pooling group to the plots in a different colour |
ten diagnostic plots for pooling groups
Anthony Hammond
#Form a gauged pooling group and plot the diagnostics with gauged = TRUE Pool.28015 <- Pool(GetCDs(28015)) DiagPlots(Pool.28015, gauged = TRUE) #Form an ugauged pooling group and plot the diagnostics Pool.28015 <- Pool(GetCDs(28015), exclude = 28015) DiagPlots(Pool.28015)#Form a gauged pooling group and plot the diagnostics with gauged = TRUE Pool.28015 <- Pool(GetCDs(28015)) DiagPlots(Pool.28015, gauged = TRUE) #Form an ugauged pooling group and plot the diagnostics Pool.28015 <- Pool(GetCDs(28015), exclude = 28015) DiagPlots(Pool.28015)
Provides donor adjustment candidates, descriptors, and results in order of the proximity to the centroid of the subject catchment.
DonAdj(CDs = NULL, x, y, QMEDscd = NULL, alpha = TRUE, rows = 10, d2 = NULL)DonAdj(CDs = NULL, x, y, QMEDscd = NULL, alpha = TRUE, rows = 10, d2 = NULL)
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
x |
catchment centroid easting (for when CDs isn't used) |
y |
catchment centroid northing (for when CDs isn't used) |
QMEDscd |
QMED estimate for the catchment of interest (for when CDs isn't used) |
alpha |
logical argument with a default of TRUE. If FALSE the exponent of the donor adjustment equation is set to one |
rows |
number of sites provided; default is 10 |
d2 |
a numeric vector of length two; the two site references for the donor catchments chosen for the two donor case |
When d2 is FALSE the results for single donor adjustment are in the final column headed 'QMED.adj' for each site. If alpha is set to FALSE, the results in this column are from the same donor equation but with an exponent of 1. The donor adjustment method is as outlined in Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation. The method for two donors is outlined in 'Kjeldsen, T. (2019). Adjustment of QMED in ungauged catchments using two donor sites. Circulation - The Newsletter of the British Hydrological Society, 4'. When two donors are used, only the result is returned, rather than donor candidates. The QMEDfse column provides the gauged factorial standard error for the median of the annual maximum sample. It is worth considering this when choosing a donor site (a high value indicates a poor donor). When choosing between two donors, the site with a lower QMEDfse would be an appropriate choice (all else being equal). The QMEDfse is calculated with the QMEDfseSS() function.
A data.frame with rownames of site references and columns of catchment descriptors, distance from subect site, and associated results. When two donors are used, only the resulting adjusted QMED is returned
Anthony Hammond
#Get some CDs and output candidate donor sites CDs.54022 <- GetCDs(54022) DonAdj(CDs.54022) #Get results with inputs of x,y, and QMEDscd DonAdj(x = 283261, y = 288067, QMEDscd = 17.931) #Get a result with two donors DonAdj(CDs.54022, d2 = c(54092, 54091))#Get some CDs and output candidate donor sites CDs.54022 <- GetCDs(54022) DonAdj(CDs.54022) #Get results with inputs of x,y, and QMEDscd DonAdj(x = 283261, y = 288067, QMEDscd = 17.931) #Get a result with two donors DonAdj(CDs.54022, d2 = c(54092, 54091))
Calculates the probability of experiencing at least n events with a given return period (RP), over a given number of years
EncProb(n, yrs, RP, dist = "Poisson")EncProb(n, yrs, RP, dist = "Poisson")
n |
number of events |
yrs |
number of years |
RP |
return period of the events |
dist |
choice of probability distribution. Either "Poisson" or "Binomial" |
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.
A probability
Anthony Hammond
#Calculate the probability of exceeding at least one 50-yr 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-yr RP events #over a 100 year period, using the binomial distribution. EncProb(n = 2, yrs = 100, RP = 100, dist = "Binomial")#Calculate the probability of exceeding at least one 50-yr 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-yr RP events #over a 100 year period, using the binomial distribution. EncProb(n = 2, yrs = 100, RP = 100, dist = "Binomial")
A plot to inspect the distribution of ordered data
ERPlot(x, dist = "GenLog", main = NULL, Pars = NULL, GF = NULL, ERType = 1)ERPlot(x, dist = "GenLog", main = NULL, Pars = NULL, GF = NULL, ERType = 1)
x |
numeric vector. A sample for inspection |
dist |
a choice of distribution. The choices are "GenLog" (the default), "GEV", "Kappa3,"Gumbel", and "GenPareto" |
main |
a character string to change the default title, which is the distribution choice. |
Pars |
a vector of parameters for the distribution. In the order of location, scale, & shape (ignoring the latter if Gumbel). If left null the parameters are estimated from x. |
GF |
a vector of length growth curve parameters, in the order of; Lcv, LSkew and Median (ignoring the LSkew if Gumbel). |
ERType |
Either 1, 2. If it is the defaul 1 then ranks are plotted on the x axis and percentage difference of modelled from observed is plotted on the y axis. |
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.
The extreme rank plot as described in the details
Anthony Hammond
#Get an AMAX sample and plot AM.27083 <- GetAM(27083) ERPlot(AM.27083$Flow) #Get a pooled estimate of Lcv & LSkew to use with the GF argument QuickResults(GetCDs(27083), gauged = TRUE) #Use the resulting Lcv, Lskew and RP2 for the GF argument and change the title ERPlot(AM.27083$Flow, main = "Site 27083 pooled comparison", GF = c(0.23, 0.17, 12))#Get an AMAX sample and plot AM.27083 <- GetAM(27083) ERPlot(AM.27083$Flow) #Get a pooled estimate of Lcv & LSkew to use with the GF argument QuickResults(GetCDs(27083), gauged = TRUE) #Use the resulting Lcv, Lskew and RP2 for the GF argument and change the title ERPlot(AM.27083$Flow, main = "Site 27083 pooled comparison", GF = c(0.23, 0.17, 12))
Plots the extreme value frequency curve or growth curve with observed sample points.
EVPlot( x, dist = "GenLog", scaled = TRUE, Title = "Extreme value plot", ylabel = NULL, LineName = NULL, Unc = TRUE )EVPlot( x, dist = "GenLog", scaled = TRUE, Title = "Extreme value plot", ylabel = NULL, LineName = NULL, Unc = TRUE )
x |
a numeric vector. The sample of interest |
dist |
a choice of distribution. "GEV", "GenLog", "Kappa3","Gumbel" or "GenPareto". The default is "GenLog" |
scaled |
logical argument with a default of TRUE. If TRUE the plot is a growth curve (scaled by the QMED). If FALSE, the plot is a frequency curve |
Title |
a character string. The user chosen plot title. The default is "Extreme value plot" |
ylabel |
a character string. The user chosen label for the y axis. The default is "Q/QMED" if scaled = TRUE and "Discharge (m3/s)" if scaled = FALSE |
LineName |
a character string. User chosen label for the plotted curve |
Unc |
logical argument with a default of TRUE. If TRUE, 95 percent uncertainty intervals are plotted. |
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.
An extreme value plot (frequency or growth curve) with intervals to quantify uncertainty
Anthony Hammond
#Get an AMAX sample and plot the growth curve with the GEV distribution AM.203018 <- GetAM(203018) EVPlot(AM.203018$Flow, dist = "GEV")#Get an AMAX sample and plot the growth curve with the GEV distribution AM.203018 <- GetAM(203018) EVPlot(AM.203018$Flow, dist = "GEV")
Functionality to add extra lines or points to an extreme value plot (derived from the EVPlot function).
EVPlotAdd( Pars, dist = "GenLog", Name = "Adjusted", MED = NULL, xyleg = NULL, col = "red", lty = 1, pts = NULL, ptSym = NULL )EVPlotAdd( Pars, dist = "GenLog", Name = "Adjusted", MED = NULL, xyleg = NULL, col = "red", lty = 1, pts = NULL, ptSym = NULL )
Pars |
a numeric vector of length two. The first is the Lcv (linear coefficient of variation) and the second is the Lskew (linear skewness). |
dist |
distribution name with a choice of "GenLog", "GEV", "GenPareto", "Kappa3", and "Gumbel" |
Name |
character string. User chosen name for points or line added (for the legend) |
MED |
The two year return level. Necessary In the case where the EV plot is not scaled |
xyleg |
a numeric vector of length two. They are the x and y position of the symbol and text to be added to the legend. |
col |
The colour of the points of line that have been added |
lty |
An integer. The type of line added |
pts |
A numeric vector. An annual maximum sample, for example. This is for points to be added |
ptSym |
An integer. The symbol of the points to be added |
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).
Additional, user specified line or points to an extreme value plot derived from the EVPlot function.
Anthony Hammond
#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 & red) for the generalised logistic distribution #first get the Lcv and Lskew using the LMoments 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 darkgreen 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")#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 & red) for the generalised logistic distribution #first get the Lcv and Lskew using the LMoments 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 darkgreen 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")
Plots the extreme value frequency curve or growth curve for gauged or ungauged pooled groups
EVPool( x, AMAX = NULL, gauged = FALSE, dist = "GenLog", QMED = NULL, Title = "Pooled growth curve", UrbAdj = FALSE, CDs )EVPool( x, AMAX = NULL, gauged = FALSE, dist = "GenLog", QMED = NULL, Title = "Pooled growth curve", UrbAdj = FALSE, CDs )
x |
pooling group derived from the Pool() function |
AMAX |
the AMAX sample to be plotted in the case of gauged. If NULL, & gauged equals TRUE, the AMAX from the first site in the pooling group is used |
gauged |
logical argument with a default of FALSE. If FALSE, the plot is the ungauged pooled curve accompanied by the single site curves of the group members. If TRUE, the plot is the gauged curve and single site curve with the observed points added |
dist |
a choice of distribution. Choices are "GEV", "GenLog", "Kappa3", or "Gumbel". The default is "GenLog" |
QMED |
a chosen QMED to convert the curve from a growth curve to the frequency curve |
Title |
a character string. The user chosen plot title. The default is "Pooled growth curve" |
UrbAdj |
a logical argument with a default of FALSE. If TRUE and urban adjustment is applied to the pooled growth curve |
CDs |
catchment descriptors derived from either GetCDs or CDsXML. Only necessary if UrbAdj is TRUE |
An extreme value plot for gauged or ungauged pooling groups
Anthony Hammond
#Get some CDs, form an ungauged pooling group and apply EVPlot. 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 105.5. PoolG.28015 <- Pool(CDs.28015) EVPool(PoolG.28015, gauged = TRUE, Title = "Gauged frequency curve - Site 28015", QMED = 9.8) #Pretend we have an extra AMAX for the gauge. Amend the pooling group Lcv and LSkew #for the site accordingly then apply EVPool with the updated AMAX. #Firstly, get the AMAX sample AM.28015 <- GetAM(28015) #Add an extra AMAX flow of 15m3/s Append28015 <- append(AM.28015$Flow, 15) #Amend the Lcv and Lskew in the pooling group PoolG.28015[1, c(16, 17)] <- c(Lcv(Append28015), LSkew(Append28015)) #Now plot gauged with the updated AMAX EVPool(PoolG.28015, AMAX = Append28015, gauged = TRUE)#Get some CDs, form an ungauged pooling group and apply EVPlot. 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 105.5. PoolG.28015 <- Pool(CDs.28015) EVPool(PoolG.28015, gauged = TRUE, Title = "Gauged frequency curve - Site 28015", QMED = 9.8) #Pretend we have an extra AMAX for the gauge. Amend the pooling group Lcv and LSkew #for the site accordingly then apply EVPool with the updated AMAX. #Firstly, get the AMAX sample AM.28015 <- GetAM(28015) #Add an extra AMAX flow of 15m3/s Append28015 <- append(AM.28015$Flow, 15) #Amend the Lcv and Lskew in the pooling group PoolG.28015[1, c(16, 17)] <- c(Lcv(Append28015), LSkew(Append28015)) #Now plot gauged with the updated AMAX EVPool(PoolG.28015, AMAX = Append28015, gauged = TRUE)
A function to plot flow duration curves for a single flow series or flow duration curves from multiple flow series.
FlowDurationCurve( x = NULL, main = "Flow duration curve", CompareCurves = NULL, LegNames = NULL, Cols = NULL, AddQs = NULL, ReturnData = FALSE )FlowDurationCurve( x = NULL, main = "Flow duration curve", CompareCurves = NULL, LegNames = NULL, Cols = NULL, AddQs = NULL, ReturnData = FALSE )
x |
a dataframe with date in the first column and numeric (flow) in the second. |
main |
A title for the plot. The default is Flow duration curve. |
CompareCurves |
A user supplied list where each element is a numeric vector (each a flow series). This is useful for when you want to compare curves from multiple flow series'. |
LegNames |
User supplied names for the legend. This only works when the CompareCurves argument is used. The default is Curve1, Curve2...CurveN. |
Cols |
User supplied vector of colours. This only works when the CompareCurves argument is used. The default is the Zissou 1 palette. |
AddQs |
The is to add additional flows and associated horizontal plot lines to the plot. It should be a single numeric value or a vector, for example c(25, 75, 100). |
ReturnData |
Logical argument with a default of FALSE. When TRUE, a dataframe is returned with the data from the plot. |
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. Or a list of flow series' (vectors) can be applied for a plot comparing the individual flow duration curves
If a dataframe of date in the first column and flow in the second is applied with the x argument a plot of a the flow duration curves for winter, summer and annual 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.
Anthony Hammond
# Plot a flow duration curve for the Thames at Kingston Oct 2000 - Sep 2015. FlowDurationCurve(ThamesPQ[,c(1,3)]) #Add a couple of extra 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"))# Plot a flow duration curve for the Thames at Kingston Oct 2000 - Sep 2015. FlowDurationCurve(ThamesPQ[,c(1,3)]) #Add a couple of extra 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"))
A function to separate baseflow from runoff.
FlowSplit( x, BaseQUpper = NULL, AdjUp = NULL, ylab = "Value", xlab = "Time index" )FlowSplit( x, BaseQUpper = NULL, AdjUp = NULL, ylab = "Value", xlab = "Time index" )
x |
A numeric vector (your flow series / hydrograph). |
BaseQUpper |
Numeric value which is an upper level of baseflow (i.e. the baseflow will not extend above this level). The default is the mean of x. It can be set arbitrarily high so that the baseflow joins all low points/troughs in the hydrograph. |
AdjUp |
A numeric value between 0 and 0.5. This allows the user to adjust the baseflow up the falling limb/s of the hydrogaph. With 0.05 being a small upward adjustment and 0.49 being a large upward adjustment. |
ylab |
Label for the yaxis (character string). The default is "value", |
xlab |
Label for the xaxis (character string). The default is "Time index". |
The function 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 hydrgraph/s can be raised using the AdjUp argument. The function works for any sampling frequency and arbitrary hydrograph length (although more suited for sub-annual to event scale). This function is not design for deriving long term baseflow index. It could be used for such a purpose but careful consideration is required for the BaseQUpper argument especially for comparison across river locations. If baseflow index is required the BFI function (with daily mean flow) may be more suitable.
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.
Anthony Hammond
# We'll extract a wet six month period on the Thames during the 2006-2007 hydrological year ThamesQ <- subset(ThamesPQ[,c(1,3)], Date >= "2006-11-04" & Date <= "2007-05-06") #Then apply the flow split with default settings QSplit <- FlowSplit(ThamesQ$Q) #Now do it with an upper baseflow level of 100m3/s QSplit <- FlowSplit(ThamesQ$Q, BaseQUpper = 100) # First we'll use the DesHydro function to pick out a reasonable looking event from the Thames flow. Q <- DesHydro(ThamesPQ[,c(1,3)], Plot = FALSE, EventSep = 15) Q <- Q$AllScaledHydrographs$hydro7 #Then we'll use our flow split function FlowSplit(Q) #Next we will get a single peaked "idealised" hydrograph using the ReFH function. QReFH <- ReFH(GetCDs(15006)) QReFH <- QReFH[[2]]$TotalFlow #Now use the function with and without an upward adjustment of the baseflow on the falling limb. QFlowSplit <- FlowSplit(QReFH) QFlowSplit <- FlowSplit(QReFH, AdjUp = 0.15)# We'll extract a wet six month period on the Thames during the 2006-2007 hydrological year ThamesQ <- subset(ThamesPQ[,c(1,3)], Date >= "2006-11-04" & Date <= "2007-05-06") #Then apply the flow split with default settings QSplit <- FlowSplit(ThamesQ$Q) #Now do it with an upper baseflow level of 100m3/s QSplit <- FlowSplit(ThamesQ$Q, BaseQUpper = 100) # First we'll use the DesHydro function to pick out a reasonable looking event from the Thames flow. Q <- DesHydro(ThamesPQ[,c(1,3)], Plot = FALSE, EventSep = 15) Q <- Q$AllScaledHydrographs$hydro7 #Then we'll use our flow split function FlowSplit(Q) #Next we will get a single peaked "idealised" hydrograph using the ReFH function. QReFH <- ReFH(GetCDs(15006)) QReFH <- QReFH[[2]]$TotalFlow #Now use the function with and without an upward adjustment of the baseflow on the falling limb. QFlowSplit <- FlowSplit(QReFH) QFlowSplit <- FlowSplit(QReFH, AdjUp = 0.15)
Estimated quantiles as a function of return period (RP) and vice versa, directly from the data
GenLogAM(x, RP = 100, q = NULL)GenLogAM(x, RP = 100, q = NULL)
x |
numeric vector (block maxima sample) |
RP |
return period (default = 100) |
q |
quantile (magnitude of variable) |
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. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'.
quantile as a function of RP or vice versa.
Anthony Hammond
#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 600m3/s discharge GenLogAM(AM.27090$Flow, q = 600)#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 600m3/s discharge GenLogAM(AM.27090$Flow, q = 600)
Estimated quantiles as function of return period (RP) and vice versa, from user input parameters
GenLogEst(loc, scale, shape, q = NULL, RP = 100)GenLogEst(loc, scale, shape, q = NULL, RP = 100)
loc |
location parameter |
scale |
scale parameter |
shape |
shape parameter |
q |
quantile. magnitude of the variable under consideration |
RP |
return period |
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.
quantile as a function of RP or vice versa
Anthony Hammond
#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) GenLogPars(AM.27090$Flow) #Store parameters in an object Pars <- as.numeric(GenLogPars(AM.27090$Flow)) #get estimate of 50-yr flow GenLogEst(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge GenLogEst(Pars[1], Pars[2], Pars[3], q = 600)#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) GenLogPars(AM.27090$Flow) #Store parameters in an object Pars <- as.numeric(GenLogPars(AM.27090$Flow)) #get estimate of 50-yr flow GenLogEst(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge GenLogEst(Pars[1], Pars[2], Pars[3], q = 600)
Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)
GenLogGF(lcv, lskew, RP)GenLogGF(lcv, lskew, RP)
lcv |
linear coefficient of variation |
lskew |
linear skewness |
RP |
return period |
Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.
Generalised logistic estimated growth factor
Anthony Hammond
#Estimate the 50-year growth factors from an Lcv and Lskew of 0.17 and 0.04, respectively. GenLogGF(0.17, 0.04, RP = 50)#Estimate the 50-year growth factors from an Lcv and Lskew of 0.17 and 0.04, respectively. GenLogGF(0.17, 0.04, RP = 50)
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)
GenLogPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)GenLogPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)
x |
numeric vector. The sample |
mle |
logical argument with a default of FALSE. If FALSE the parameters are estimated with Lmoments, if TRUE the parameters are estimated by maximum likelihood estimation. |
L1 |
first Lmoment |
LCV |
linear coefficient of variation |
LSKEW |
linear skewness |
The L-moment estimated parameters are by the method detailed in 'Hosking J. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'
Parameter estimates (location, scale, shape)
Anthony Hammond
#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GenLogPars(AM.27090$Flow) #Estimate parameters using MLE GenLogPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(AM.27090$Flow) #store linear moments in an object LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,5,6)] GenLogPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GenLogPars(AM.27090$Flow) #Estimate parameters using MLE GenLogPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(AM.27090$Flow) #store linear moments in an object LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,5,6)] GenLogPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])
Estimated quantiles as function of return period (RP) and vice versa, from user input parameters
GenParetoEst(loc, scale, shape, q = NULL, RP = 100, ppy = 1)GenParetoEst(loc, scale, shape, q = NULL, RP = 100, ppy = 1)
loc |
location parameter |
scale |
scale parameter |
shape |
shape parameter |
q |
quantile. magnitude of the variable under consideration |
RP |
return period |
ppy |
peaks per year. Default is one |
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.
quantile as a function of RP or vice versa
Anthony Hammond
#Get a POT sample, estimate the parameters, and estimate 50-year RP ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPars(ThamesPOT$peak) #Store parameters in an object Pars <- as.numeric(GenParetoPars(ThamesPOT$peak)) #get estimate of 50-yr flow GenParetoEst(Pars[1], Pars[2], Pars[3], ppy = 1.867, RP = 50) #Estimate the RP for a 600m3/s discharge GenParetoEst(Pars[1], Pars[2], Pars[3], ppy = 1.867, q = 600)#Get a POT sample, estimate the parameters, and estimate 50-year RP ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPars(ThamesPOT$peak) #Store parameters in an object Pars <- as.numeric(GenParetoPars(ThamesPOT$peak)) #get estimate of 50-yr flow GenParetoEst(Pars[1], Pars[2], Pars[3], ppy = 1.867, RP = 50) #Estimate the RP for a 600m3/s discharge GenParetoEst(Pars[1], Pars[2], Pars[3], ppy = 1.867, q = 600)
Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)
GenParetoGF(lcv, lskew, RP, ppy = 1)GenParetoGF(lcv, lskew, RP, ppy = 1)
lcv |
linear coefficient of variation |
lskew |
linear skewness |
RP |
return period |
ppy |
peaks per year |
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.
Generalised Pareto estimated growth factor
Anthony Hammond
#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-yr flood GenParetoGF(Lcv(TPOT$peak), LSkew(TPOT$peak), RP = 100, ppy = 1.867)*median(TPOT$peak)#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-yr flood GenParetoGF(Lcv(TPOT$peak), LSkew(TPOT$peak), RP = 100, ppy = 1.867)*median(TPOT$peak)
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)
GenParetoPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)GenParetoPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)
x |
numeric vector. The sample |
mle |
logical argument with a default of FALSE. If FALSE the parameters are estimated with Lmoments, if TRUE the parameters are estimated by maximum likelihood estimation |
L1 |
first Lmoment |
LCV |
linear coefficient of variation |
LSKEW |
linear skewness |
The L-moment estimated parameters are by the method detailed in 'Hosking J. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'
Parameter estimates (location, scale, shape)
Anthony Hammond
#Get a peaks over threshold sample and estimate the parameters using Lmoments ThamesPOT <- ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPars(ThamesPOT$peak) #Estimate parameters using MLE GenParetoPars(ThamesPOT$peak, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(ThamesPOT$peak) #store linear moments in an object LPars <- as.numeric(LMoments(ThamesPOT$peak))[c(1,5,6)] GenParetoPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])#Get a peaks over threshold sample and estimate the parameters using Lmoments ThamesPOT <- ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPars(ThamesPOT$peak) #Estimate parameters using MLE GenParetoPars(ThamesPOT$peak, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(ThamesPOT$peak) #store linear moments in an object LPars <- as.numeric(LMoments(ThamesPOT$peak))[c(1,5,6)] GenParetoPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])
Estimated quantiles as function of return period (RP) and vice versa, directly from the data
GenParetoPOT(x, ppy = 1, RP = 100, q = NULL)GenParetoPOT(x, ppy = 1, RP = 100, q = NULL)
x |
numeric vector (block maxima sample) |
ppy |
peaks per year |
RP |
return period (default = 100) |
q |
quantile (magnitude of variable) |
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. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'.
quantile as a function of RP or vice versa
Anthony Hammond
#Get a POT series and estimate the 50-year RP ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPOT(ThamesPOT$peak, ppy = 1.867, RP = 50) #Estimate the RP for a 600m3/s discharge GenParetoPOT(ThamesPOT$peak, ppy = 1.867, q = 600)#Get a POT series and estimate the 50-year RP ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) GenParetoPOT(ThamesPOT$peak, ppy = 1.867, RP = 50) #Estimate the RP for a 600m3/s discharge GenParetoPOT(ThamesPOT$peak, ppy = 1.867, q = 600)
Extracts the annual maximum peak flow sample and associated dates for the site of interest.
GetAM(ref)GetAM(ref)
ref |
the site reference of interest (numeric) |
A data.frame with columns; Date, Flow, and id
Anthony Hammond
#Get an AMAX sample and display it in the console GetAM(203018) #Save an AMAX sample as an object AM.203018 <- GetAM(203018)#Get an AMAX sample and display it in the console GetAM(203018) #Save an AMAX sample as an object AM.203018 <- GetAM(203018)
Extracts the catchment descriptors for a site of interest from the National River Flow Archive. If the site is considered suitable for QMED and pooling the CDs are extracted from the QMEDData data.frame. Otherwise they are extracted using the NRFA API.
GetCDs(x)GetCDs(x)
x |
the site reference of interest (numeric) |
A data.frame with columns; Descriptor and Value.
Anthony Hammond
#Get CDs and display in the console CDs.203018 <- GetCDs(203018) CDs.203018#Get CDs and display in the console CDs.203018 <- GetCDs(203018) CDs.203018
Function to extract flow or level data from the Environment Agency's Hydrology Data Explorer.
GetDataEA_QH( Lat = 54, Lon = -2.25, Range = 20, RiverName = NULL, WISKI_ID = NULL, From = NULL, To = NULL, Type = "flow", Period = "DailyMean" )GetDataEA_QH( Lat = 54, Lon = -2.25, Range = 20, RiverName = NULL, WISKI_ID = NULL, From = NULL, To = NULL, Type = "flow", Period = "DailyMean" )
Lat |
Latitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
Lon |
Longitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
Range |
Radius of search when using latitude and longitude inputs (km). |
RiverName |
Name of the river along which you want to search for gauges. Character string. |
WISKI_ID |
The WISKI ID for the gauge from which you want to obtain data (character string). Note that sometimes a preceding zero, which is not returned via the API, is needed. If the data extraction fails, this may be the cause and you can resolve it by including the preceding zero in the WISKI_ID. |
From |
Date for start of data extraction in the format of "2015-12-02". If NULL the start date of the available data is used. |
To |
Date for the end of data extraction in the format of "2015-12-02". If NULL end date of the available data is used. |
Type |
The variable to extract, either "flow" or "level" |
Period |
The sampling rate of the data you want. Either "DailyMax", "DailyMean", "Hourly", "15Mins". |
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.
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.
Anthony Hammond
#Find gauges on the river Tame ## Not run: GetDataEA_QH(RiverName = "Tame") #Find gauges within 10km of a latlon grid reference somewhere near the #centre of Dartmoor ## Not run: GetDataEA_QH(Lat = 50.6, Lon = -3.9, Range = 10) #Get all available daily maximum flow data from the Bellever gauge on the #East Dart River. ## Not run: BelleverMax <- GetDataEA_QH(WISKI_ID = "SX67F051") #Get 15-minute data from the Bellever for the Novermeber 2024 event ## Not run: BelleverNov2024 <- GetDataEA_QH(WISKI_ID = "SX67F051", From = "2024-11-23", To = "2024-11-25", Period = "15Mins") ## End(Not run)#Find gauges on the river Tame ## Not run: GetDataEA_QH(RiverName = "Tame") #Find gauges within 10km of a latlon grid reference somewhere near the #centre of Dartmoor ## Not run: GetDataEA_QH(Lat = 50.6, Lon = -3.9, Range = 10) #Get all available daily maximum flow data from the Bellever gauge on the #East Dart River. ## Not run: BelleverMax <- GetDataEA_QH(WISKI_ID = "SX67F051") #Get 15-minute data from the Bellever for the Novermeber 2024 event ## Not run: BelleverNov2024 <- GetDataEA_QH(WISKI_ID = "SX67F051", From = "2024-11-23", To = "2024-11-25", Period = "15Mins") ## End(Not run)
Extract rainfall data from the Environment Agency's Hydrology Data Explorer.
GetDataEA_Rain( Lat = 54, Lon = -2, Range = 10, WISKI_ID = NULL, Period = "Daily", From = NULL, To = NULL )GetDataEA_Rain( Lat = 54, Lon = -2, Range = 10, WISKI_ID = NULL, Period = "Daily", From = NULL, To = NULL )
Lat |
Latitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
Lon |
Longitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
Range |
The radius (km) from the point of interest (Lat, Lon) for which the user wants rain gauge information (currently it only seems to work to just over 20km). |
WISKI_ID |
The WISKI identification (as "character") for the rain gauge of interest |
Period |
The sampling rate of the rainfall in hours. Either "Daily", "15Mins", "Hourly". |
From |
The start date of the data extraction in the form of "YYYY-MM-DD". To get data from the first date available leave as NULL. |
To |
The end date of the data extraction in the form of "YYYY-MM-DD". To get data from the first date available leave as NULL. |
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.
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.
Anthony Hammond
#Get information about available rain gauges. #within a 10km radius of Lat = 54.5, Lon = -3.2 ## Not run: GetDataEA_Rain(Lat = 54.5, Lon = -3.2) #Now we'll use the WISKI reference for the Honister rain gauge # to get some hourly rain data for the Dec 2015 ## Not run: HonisterDec2015 <- GetDataEA_Rain(WISKI_ID = "592463", Period = "Hourly", From = "2015-12-01", To = "2015-12-31") ## End(Not run) #Now we'll have a look at the top of the data and plot it ## Not run: head(HonisterDec2015) ## Not run: plot(HonisterDec2015, type = "h", ylab = "Rainfall (mm)")#Get information about available rain gauges. #within a 10km radius of Lat = 54.5, Lon = -3.2 ## Not run: GetDataEA_Rain(Lat = 54.5, Lon = -3.2) #Now we'll use the WISKI reference for the Honister rain gauge # to get some hourly rain data for the Dec 2015 ## Not run: HonisterDec2015 <- GetDataEA_Rain(WISKI_ID = "592463", Period = "Hourly", From = "2015-12-01", To = "2015-12-31") ## End(Not run) #Now we'll have a look at the top of the data and plot it ## Not run: head(HonisterDec2015) ## Not run: plot(HonisterDec2015, type = "h", ylab = "Rainfall (mm)")
Extracts regional mean temperature or rainfall from the met office UK & regional series. The total duration of bright sunshine is also available.
GetDataMetOffice(Variable, Region)GetDataMetOffice(Variable, Region)
Variable |
Either Tmean, Rainfall, or Sunshine |
Region |
One of "UK", "England", "Wales", "Scotland", "Northern_Ireland", "England_and_Wales", "England_N", "England_S", "Scotland_N", "Scotland_E", "Scotland_W", "England_E_and_NE", "England_NW_and_N_Wales", "Midlands", "East_Anglia", "England_SW_and_S_Wales", "England_SE_and_Central_S". |
The function returns time series data from the 19th century through to the present month.
A data.frame with 18 columns; year, months, seasons, and annual. Rows then represent each year of the timeseries.
Anthony Hammond
#Get the Rainfall series for the UK ## Not run: UKRain <- GetDataMetOffice(Variable = "Rainfall", Region = "UK") #Now we'll get mean temperature data for East Anglia ## Not run: TempEastAnglia <- GetDataMetOffice(Variable = "Tmean", Region = "East_Anglia")#Get the Rainfall series for the UK ## Not run: UKRain <- GetDataMetOffice(Variable = "Rainfall", Region = "UK") #Now we'll get mean temperature data for East Anglia ## Not run: TempEastAnglia <- GetDataMetOffice(Variable = "Tmean", Region = "East_Anglia")
Extracts NRFA data using the API.
GetDataNRFA(ID, Type = "Q")GetDataNRFA(ID, Type = "Q")
ID |
ID number of the gauge of interest. |
Type |
Type of data required. One of "Q", "P", "PQ", "Gaugings", "AMAX", "POT", or "Catalogue". |
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" it will return a dataframe of all the NRFA gauges, associated details, comments, and descriptors.
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).
Anthony Hammond
#Get the concurrent rainfall and mean flow series for the Tay at Ballathie (site 15006). ## Not run: BallathiePQ <- GetDataNRFA(15006, "PQ") #Now we'll get the gaugings ## Not run: BallathieGaugings <- GetDataNRFA(15006, "Gaugings")#Get the concurrent rainfall and mean flow series for the Tay at Ballathie (site 15006). ## Not run: BallathiePQ <- GetDataNRFA(15006, "PQ") #Now we'll get the gaugings ## Not run: BallathieGaugings <- GetDataNRFA(15006, "Gaugings")
Function to extract flow or level data from SEPA.
GetDataSEPA_QH( Lat = NULL, Lon = NULL, RiverName = NULL, Type = "Flow", StationID = NULL, Range = 20, From = NULL, To = NULL, Period = "Daily" )GetDataSEPA_QH( Lat = NULL, Lon = NULL, RiverName = NULL, Type = "Flow", StationID = NULL, Range = 20, From = NULL, To = NULL, Period = "Daily" )
Lat |
Latitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
Lon |
Longitude (as a decimal) for the centre of the search for gauges. You can convert BNG to Lat and Lon using the GridRefConvert function. |
RiverName |
Name of the river along which you want to search for gauges. Character string. |
Type |
The variable to extract, either "flow" or "level" |
StationID |
The ID for the gauge from which you want to obtain data (character string) |
Range |
Radius of search when using latitude and longitude inputs (km). |
From |
Date for start of data extraction in the format of "2015-12-02". If NULL the first date of the available data is used. |
To |
Date for the end of data extraction in the format of "2015-12-02". If NULL the present date is used (and most recent available data is returned). |
Period |
The sampling rate of the data you want. Either "Daily", "Hourly", or "15Mins". |
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. The lat and lon option will provide all flow or level gauges within a specified range (default of 50km). This provides gauged details including the StationID. You can get data from specific gauges using the StationID. Note that flow gauges also have level data available. You can get data from a date range using the From and To arguments. If the From and To arguments are left as NULL the full range of available data are returned.
If searching for gauge details with lat and lon or river name, then a dataframe is returned with necessary information to obtain flow or level data. When extracting flow or level data with a station 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.
Anthony Hammond
#Find gauges on the river Spey ## Not run: GetDataSEPA_QH(RiverName = "Spey") #Find gauges within 20km of a latlon grid reference somewhere near the centre of Scotland ## Not run: GetDataSEPA_QH(lat = 56, lon = -4, Range = 20) #Get all available daily mean flow data from the Boat o Brig gauge on the Spey ## Not run: SpeyDaily <- GetDataSEPA_QH(StationID = "37174") #Get 15-minute data from the Boat o Brig for the highest recorded peak ## Not run: BoatOBrigAug1970 <- GetDataSEPA_QH(StationID = "37174", From = "1970-08-16", To = "1970-08-19", Period = "15Mins") ## End(Not run)#Find gauges on the river Spey ## Not run: GetDataSEPA_QH(RiverName = "Spey") #Find gauges within 20km of a latlon grid reference somewhere near the centre of Scotland ## Not run: GetDataSEPA_QH(lat = 56, lon = -4, Range = 20) #Get all available daily mean flow data from the Boat o Brig gauge on the Spey ## Not run: SpeyDaily <- GetDataSEPA_QH(StationID = "37174") #Get 15-minute data from the Boat o Brig for the highest recorded peak ## Not run: BoatOBrigAug1970 <- GetDataSEPA_QH(StationID = "37174", From = "1970-08-16", To = "1970-08-19", Period = "15Mins") ## End(Not run)
Extract hourly rainfall data from SEPA's API.
GetDataSEPA_Rain( Lat = NULL, Lon = NULL, Range = 30, StationName, From = NULL, To = NULL )GetDataSEPA_Rain( Lat = NULL, Lon = NULL, Range = 30, StationName, From = NULL, To = NULL )
Lat |
Latitude of the point of interest. Provided when the user wants information about available local rain gauges |
Lon |
Longitude of the point of interest. Provided when the user wants information about available local rain gauges |
Range |
The radius (km) from the point of interest (Lat, Lon) for which the user wants a list of rain gauges (default is 30). |
StationName |
The name of the station for which you want rainfall. If you type something other than one of the available stations, the list of stations will be returned. |
From |
A start date for the data in the form of "YYYY-MM-DD". Default of NULL means the earliest available date is used |
To |
An end date for the data in the form of "YYYY-MM-DD". The default is the most recent date available. |
You can download data using the gauge name and you can find gauges within a given range using the latitude and longitude. 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.
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.
Anthony Hammond
#Get the list of available stations ## Not run: GetDataSEPA_Rain(StationName = "AnythingButAStationName") #Now we'll get rain from the Bannockburn station ## Not run: Bannockburn <- GetDataSEPA_Rain(StationName = "Bannockburn", From = "1998-10-01", To = "1998-10-31") ## End(Not run) #Now we'll have a look at the top of the data and plot it ## Not run: head(Bannockburn) ## Not run: plot(Bannockburn, type = "h", ylab = "Rainfall (mm)")#Get the list of available stations ## Not run: GetDataSEPA_Rain(StationName = "AnythingButAStationName") #Now we'll get rain from the Bannockburn station ## Not run: Bannockburn <- GetDataSEPA_Rain(StationName = "Bannockburn", From = "1998-10-01", To = "1998-10-31") ## End(Not run) #Now we'll have a look at the top of the data and plot it ## Not run: head(Bannockburn) ## Not run: plot(Bannockburn, type = "h", ylab = "Rainfall (mm)")
Provides QMED (median annual maximum flow) from a site suitable for QMED, using the site reference.
GetQMED(x)GetQMED(x)
x |
the gauged reference |
the median annual maximum
Anthony Hammond
#Get the observed QMED from sites 55002 GetQMED(55002)#Get the observed QMED from sites 55002 GetQMED(55002)
Estimated quantiles as function of return period (RP) and vice versa, directly from the data
GEVAM(x, RP = 100, q = NULL)GEVAM(x, RP = 100, q = NULL)
x |
numeric vector (block maxima sample) |
RP |
return period (default = 100) |
q |
quantile (magnitude of variable) |
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. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'.
quantile as a function of RP or vice versa.
Anthony Hammond
#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 600m3/s discharge GEVAM(AM.27090$Flow, q = 600)#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 600m3/s discharge GEVAM(AM.27090$Flow, q = 600)
Estimated quantiles as function of return period (RP) and vice versa, from user input parameters
GEVEst(loc, scale, shape, q = NULL, RP = 100)GEVEst(loc, scale, shape, q = NULL, RP = 100)
loc |
location parameter |
scale |
scale parameter |
shape |
shape parameter |
q |
quantile. magnitude of the variable under consideration |
RP |
return period |
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.
quantile as a function of RP or vice versa
Anthony Hammond
#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) GEVPars(AM.27090$Flow) #Store parameters in an object Pars <- as.numeric(GEVPars(AM.27090$Flow)) #get estimate of 50-yr flow GEVEst(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge GEVEst(Pars[1], Pars[2], Pars[3], q = 600)#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) GEVPars(AM.27090$Flow) #Store parameters in an object Pars <- as.numeric(GEVPars(AM.27090$Flow)) #get estimate of 50-yr flow GEVEst(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge GEVEst(Pars[1], Pars[2], Pars[3], q = 600)
Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)
GEVGF(lcv, lskew, RP)GEVGF(lcv, lskew, RP)
lcv |
linear coefficient of variation |
lskew |
linear skewness |
RP |
return period |
Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.
Generalised extreme value estimated growth factor
Anthony Hammond
#Estimate the 50-year growth factors from Lcv = 0.17 and Lskew = 0.04 GEVGF(0.17, 0.04, RP = 50)#Estimate the 50-year growth factors from Lcv = 0.17 and Lskew = 0.04 GEVGF(0.17, 0.04, RP = 50)
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)
GEVPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)GEVPars(x = NULL, mle = FALSE, L1, LCV, LSKEW)
x |
numeric vector. The sample |
mle |
logical argument with a default of FALSE. If FALSE the parameters are estimated with Lmoments, if TRUE the parameters are estimated by maximum likelihood estimation |
L1 |
first Lmoment |
LCV |
linear coefficient of variation |
LSKEW |
linear skewness |
The L-moment estimated parameters are by the method detailed in 'Hosking J. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'
Parameter estimates (location, scale, shape)
Anthony Hammond
#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GEVPars(AM.27090$Flow) #Estimate parameters using MLE GEVPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(AM.27090$Flow) #store linear moments in an object LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,5,6)] GEVPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GEVPars(AM.27090$Flow) #Estimate parameters using MLE GEVPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1, Lcv and Lskew LMoments(AM.27090$Flow) #store linear moments in an object LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,5,6)] GEVPars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])
compares the RMSE of four distribution fits for a single AMAX sample.
GoFCompare(x)GoFCompare(x)
x |
a numeric vector (your AMAX sample) |
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.
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.
Anthony Hammond
# Get an AMAX sample then compare the fit.. AM15006 <- GetAM(15006) GoFCompare(AM15006$Flow)# Get an AMAX sample then compare the fit.. AM15006 <- GetAM(15006) GoFCompare(AM15006$Flow)
compares the RMSE of four distribution fits for a pooling group.
GoFComparePool(x)GoFComparePool(x)
x |
a numeric vector (your AMAX sample) |
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.
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.
Anthony Hammond
# Get a pooling group then compare the fit.. Pool60009 <- Pool(GetCDs(60009)) GoFComparePool(Pool60009)# Get a pooling group then compare the fit.. Pool60009 <- Pool(GetCDs(60009)) GoFComparePool(Pool60009)
Estimated quantiles as a function of return period (RP) and vice versa, directly from the data
GumbelAM(x, RP = 100, q = NULL)GumbelAM(x, RP = 100, q = NULL)
x |
numeric vector (block maxima sample) |
RP |
return period (default = 100) |
q |
quantile (magnitude of variable) |
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. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'.
quantile as a function of RP or vice versa.
Anthony Hammond
#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 600m3/s discharge GumbelAM(AM.27090$Flow, q = 600)#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 600m3/s discharge GumbelAM(AM.27090$Flow, q = 600)
Estimated quantiles as function of return period (RP) and vice versa, from user input parameters
GumbelEst(loc, scale, q = NULL, RP = 100)GumbelEst(loc, scale, q = NULL, RP = 100)
loc |
location parameter |
scale |
scale parameter |
q |
quantile. magnitude of the variable under consideration |
RP |
return period |
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.
quantile as a function of RP or vice versa
Anthony Hammond
#Get an annual maximum sample, estimate the parameters and estimate 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 600m3/s discharge GumbelEst(Pars[1], Pars[2], q = 600)#Get an annual maximum sample, estimate the parameters and estimate 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 600m3/s discharge GumbelEst(Pars[1], Pars[2], q = 600)
Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)
GumbelGF(lcv, RP)GumbelGF(lcv, RP)
lcv |
linear coefficient of variation |
RP |
return period |
Growth factors are calculated by the method outlined in the Flood Estimation Handbook, volume 3, 1999.
Gumbel estimated growth factor
Anthony Hammond
#Estimate the 50-year growth factors from an Lcv of 0.17. GumbelGF(0.17, RP = 50)#Estimate the 50-year growth factors from an Lcv of 0.17. GumbelGF(0.17, RP = 50)
Estimated parameters from a sample (with Lmoments or maximum likelihood estimation) or from L1 (first L-moment), Lcv (linear coefficient of variation)
GumbelPars(x = NULL, mle = FALSE, L1, LCV)GumbelPars(x = NULL, mle = FALSE, L1, LCV)
x |
numeric vector. The sample |
mle |
logical argument with a default of FALSE. If FALSE the parameters are estimated with Lmoments, if TRUE the parameters are estimated by maximum likelihood estimation |
L1 |
first Lmoment |
LCV |
linear coefficient of variation |
The L-moment estimated parameters are by the method detailed in 'Hosking J. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'
Parameter estimates (location, scale)
Anthony Hammond
#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GumbelPars(AM.27090$Flow) #Estimate parameters using MLE GumbelPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1 and Lcv Pars <- as.numeric(LMoments(AM.27090$Flow)[c(1,5)]) GumbelPars(L1 = Pars[1], LCV = Pars[2])#Get an annual maximum sample and estimate the parameters using Lmoments AM.27090 <- GetAM(27090) GumbelPars(AM.27090$Flow) #Estimate parameters using MLE GumbelPars(AM.27090$Flow, mle = TRUE) #calculate Lmoments and estimate the parmeters with L1 and Lcv Pars <- as.numeric(LMoments(AM.27090$Flow)[c(1,5)]) GumbelPars(L1 = Pars[1], LCV = Pars[2])
Quantifies the heterogeneity of a pooled group
H2(x, H1 = FALSE)H2(x, H1 = FALSE)
x |
pooling group derived from the Pool() function |
H1 |
logical with a default of FALSE. If TRUE, the function applies the 'H1' version of the test (see Hosking & Wallis 1997 reference). If FALSE, the default H2 version is applied. |
The H2 measure was developed by Hosking & Wallis and can be found in their book 'Regional Frequency Analysis: an approach based on LMoments (1997). It was also adopted for use by the Flood Estimation Handbook (1999) and is described in volume 3.
A vector of two characters; the first representing the H2 score and the second stating a qualitative measure of heterogeneity.
Anthony Hammond
#Get CDs, form a pooling group and calculate H2 CDs.203018 <- GetCDs(203018) Pool.203018 <- Pool(CDs.203018) H2(Pool.203018)#Get CDs, form a pooling group and calculate H2 CDs.203018 <- GetCDs(203018) Pool.203018 <- Pool(CDs.203018) H2(Pool.203018)
Plots concurrent precipitation and discharge with precipitation along the top and discharge along the bottom
HydroPlot( x, main = "Concurrent Rainfall & Discharge", ylab = "Discharge (m3/s)", From = NULL, To = NULL, adj.y = 1.5, plw = 1, qlw = 1.8, Return = FALSE )HydroPlot( x, main = "Concurrent Rainfall & Discharge", ylab = "Discharge (m3/s)", From = NULL, To = NULL, adj.y = 1.5, plw = 1, qlw = 1.8, Return = FALSE )
x |
a data.frame with three columns in the order of date (or POSIXct), precipitation, and discharge |
main |
a character string. The user chosen plot title. The default is "Concurrent Rainfall & Discharge" |
ylab |
User choice for the y label of the plot. THe default is "Discharge (m3/s)". |
From |
a starting time for the plot. In the form of a date or POSIXct object. The default is the first row of x |
To |
an end time for the plot. In the form of a date or POSIXct object. The default is the last row of x |
adj.y |
a numeric value to adjust the closeness of the preciptation and discharge in the plot. Default is 1.5. A lower value brings them closer and a larger value further apart |
plw |
a numeric value to adjust the width of the precipitation lines. Default is one. A larger value thickens them and vice versa |
qlw |
a numeric value to adjust the width of the discharge line. Default is 1.8. A larger value thickens them and vice versa |
Return |
a logical argument with a default of FALSE. If TRUE the data-frame of time, precipitation, and flow is returned |
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.
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.
Anthony Hammond
#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)#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)
Estimated quantiles as a function of return period (RP) and vice versa, directly from the data
Kappa3AM(x, RP = 100, q = NULL)Kappa3AM(x, RP = 100, q = NULL)
x |
numeric vector (block maxima sample) |
RP |
return period (default = 100) |
q |
quantile (magnitude of variable) |
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. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'. The Kappa3 distribution is as defined by 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.
quantile as a function of RP or vice versa.
Anthony Hammond
#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 600m3/s discharge Kappa3AM(AM.27090$Flow, q = 600)#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 600m3/s discharge Kappa3AM(AM.27090$Flow, q = 600)
Estimated quantiles as function of return period (RP) and vice versa, from user input parameters
Kappa3Est(loc, scale, shape, q = NULL, RP = 100)Kappa3Est(loc, scale, shape, q = NULL, RP = 100)
loc |
location parameter |
scale |
scale parameter |
shape |
shape parameter |
q |
quantile. magnitude of the variable under consideration |
RP |
return period |
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.
quantile as a function of RP or vice versa
Anthony Hammond
#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) #Get parameters and Store as an object Pars <- as.numeric(Kappa3Pars(AM.27090$Flow)) #get estimate of 50-yr flow Kappa3Est(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge Kappa3Est(Pars[1], Pars[2], Pars[3], q = 600)#Get an annual maximum sample, estimate the parameters and estimate 50-year RP AM.27090 <- GetAM(27090) #Get parameters and Store as an object Pars <- as.numeric(Kappa3Pars(AM.27090$Flow)) #get estimate of 50-yr flow Kappa3Est(Pars[1], Pars[2], Pars[3], RP = 50) #Estimate the RP for a 600m3/s discharge Kappa3Est(Pars[1], Pars[2], Pars[3], q = 600)
Estimated growth factors as a function of return period, with inputs of Lcv & LSkew (linear coefficient of variation & linear skewness)
Kappa3GF(lcv, lskew, RP)Kappa3GF(lcv, lskew, RP)
lcv |
linear coefficient of variation |
lskew |
linear skewness |
RP |
return period |
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
Kappa3 distribution estimated growth factor
Anthony Hammond
#Get a ungauged pooled Lcv and LSkew for catchment 15006 PooledRes <- as.numeric(QuickResults(GetCDs(15006), plot = FALSE)[[2]]) #Calculate Kappa growth factor for the 100-year flood Kappa3GF(PooledRes[1], PooledRes[2], RP = 100)#Get a ungauged pooled Lcv and LSkew for catchment 15006 PooledRes <- as.numeric(QuickResults(GetCDs(15006), plot = FALSE)[[2]]) #Calculate Kappa growth factor for the 100-year flood Kappa3GF(PooledRes[1], PooledRes[2], RP = 100)
Estimated parameters from a sample (using Lmoments) or from user supplied L1 (first L-moment), Lcv (linear coefficient of variation), and LSkew (linear skewness)
Kappa3Pars(x = NULL, L1, LCV, LSKEW)Kappa3Pars(x = NULL, L1, LCV, LSKEW)
x |
numeric vector. The sample |
L1 |
first Lmoment |
LCV |
linear coefficient of variation |
LSKEW |
linear skewness |
The L-moment estimated parameters are by the method detailed in 'Hosking J. Wallis J. 1997 Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, New York'. The Kappa3 distribution is as defined by 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.
Parameter estimates (location, scale, shape)
Anthony Hammond
#Get an annual maximum sample and estimate the parameters. AM.27090 <- GetAM(27090) Kappa3Pars(AM.27090$Flow) #calculate Lmoments and estimate the parmeters with L1, L2, Lcv, and Lskew LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,2,5,6)] Kappa3Pars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])#Get an annual maximum sample and estimate the parameters. AM.27090 <- GetAM(27090) Kappa3Pars(AM.27090$Flow) #calculate Lmoments and estimate the parmeters with L1, L2, Lcv, and Lskew LPars <- as.numeric(LMoments(AM.27090$Flow))[c(1,2,5,6)] Kappa3Pars(L1 = LPars[1], LCV = LPars[2], LSKEW = LPars[3])
Calculates the Lcv from a sample of data
Lcv(x)Lcv(x)
x |
a numeric vector. The sample of interest |
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.
Numeric. The Lcv of a sample.
Anthony Hammond
#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) Lcv(AM.27051$Flow)#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) Lcv(AM.27051$Flow)
Urbanises or de-urbanises the Lcv using the methods outlined in the guidance by Wallingford HydroSolutions: 'WINFAP 4 Urban Adjustment Procedures'
LcvUrb(lcv, URBEXT2000, DeUrb = FALSE)LcvUrb(lcv, URBEXT2000, DeUrb = FALSE)
lcv |
the Lcv (numeric) |
URBEXT2000 |
quantiication of urban and suburbanisation for the subject catchment |
DeUrb |
logical argument with a default of FALSE. If set to TRUE, de-urbanisation adjustment is performed, if FALSE, urbanisation adjustment is performed |
The method for de-urbanisation isn't explicitly provided in 'WINFAP 4 Urban Adjustment Procedures', but the procedure is a re-arrangment of the urbanisation equation, solving for Lcv rather than Lcv-urban.
The urban adjust Lcv or the de-urbanised Lcv
Anthony Hammond
#Choose an urban site (site 53006) from the NRFA data then apply a de-urban #adjustment using the Lcv and URBEXT2000 displayed NRFAData[which(rownames(NRFAData) == 53006),] LcvUrb(0.21, 0.1138, DeUrb = TRUE) #Get the pooled Lmoment ratios results for catchment 53006 and apply the #urban adjustment using the pooled Lcv, and the URBEXT2000 for site 53006. CDs.53006 <- GetCDs(53006) QuickResults(CDs.53006)[[2]] LcvUrb(0.196, 0.1138)#Choose an urban site (site 53006) from the NRFA data then apply a de-urban #adjustment using the Lcv and URBEXT2000 displayed NRFAData[which(rownames(NRFAData) == 53006),] LcvUrb(0.21, 0.1138, DeUrb = TRUE) #Get the pooled Lmoment ratios results for catchment 53006 and apply the #urban adjustment using the pooled Lcv, and the URBEXT2000 for site 53006. CDs.53006 <- GetCDs(53006) QuickResults(CDs.53006)[[2]] LcvUrb(0.196, 0.1138)
Calculates the LKurtosis from a sample of data
LKurt(x)LKurt(x)
x |
a numeric vector. The sample of interest |
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.
Numeric. The LSkew of a sample.
Anthony Hammond
#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LKurt(AM.27051$Flow)#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LKurt(AM.27051$Flow)
Calculates the Lmoments and Lmoment ratios from a sample of data
LMoments(x)LMoments(x)
x |
a numeric vector. The sample of interest |
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.
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.
Anthony Hammond
#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LMoments(AM.27051$Flow)#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LMoments(AM.27051$Flow)
Adjusts the linear coefficient of variation (Lcv) and the linear skewness (LSkew) for a chosen site in a pooling group
LRatioChange(x, SiteID, lcv, lskew)LRatioChange(x, SiteID, lcv, lskew)
x |
pooling group derived with the Pool function |
SiteID |
the identification number of the site in the pooling group that is to be changed (character or integer) |
lcv |
The user supplied Lcv. numeric |
lskew |
The user supplied LSkew. numeric |
Pooling groups are formed from the NRFAData 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 permeable 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.
A new pooling group, the same as x except for the user adjusted Lcv and Lskew for the user selected site.
Anthony Hammond
# Get some catchment descriptors and create a pooling group. CDs.39001 <- GetCDs(39001) Pool.39001 <- Pool(CDs.39001, iug = TRUE) # apply the function to create a new adjusted pooling group, #changing the subject site lcv and lskew to 0.187 and 0.164, respectively Pool.39001Adj <- LRatioChange(Pool.39001, SiteID = 39001, lcv = 0.187, lskew = 0.164)# Get some catchment descriptors and create a pooling group. CDs.39001 <- GetCDs(39001) Pool.39001 <- Pool(CDs.39001, iug = TRUE) # apply the function to create a new adjusted pooling group, #changing the subject site lcv and lskew to 0.187 and 0.164, respectively Pool.39001Adj <- LRatioChange(Pool.39001, SiteID = 39001, lcv = 0.187, lskew = 0.164)
Calculates the LSkew from a sample of data
LSkew(x)LSkew(x)
x |
a numeric vector. The sample of interest |
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.
Numeric. The LSkew of a sample.
Anthony Hammond
#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LSkew(AM.27051$Flow)#Get an AMAX sample and calculate the Lmoments AM.27051 <- GetAM(27051) LSkew(AM.27051$Flow)
Urbanises or de-urbanises the LSkew using the methods outlined in the guidance by Wallingford HydroSolutions: 'WINFAP 4 Urban Adjustment Procedures'
LSkewUrb(lskew, URBEXT2000, DeUrb = FALSE)LSkewUrb(lskew, URBEXT2000, DeUrb = FALSE)
lskew |
the LSkew (numeric) |
URBEXT2000 |
quantiication of urban and suburbanisation for the subject site |
DeUrb |
logical argument with a default of FALSE. If set to TRUE, de-urbanisation adjustment is performed, if FALSE, urbanisation adjustment is performed |
The method for de-urbanisation isn't explicitly provided in 'WINFAP 4 Urban Adjustment Procedures', but the procedure is a re-arrangment of the urbanisation equation, solving for LSkew rather than LSkew-urban.
The urban adjust Lcv or the de-urbanised Lcv
Anthony Hammond
#Choose an urban site (site 53006) from the NRFA data then apply a de-urban #adjustment using the Lcv and URBEXT2000 displayed NRFAData[which(rownames(NRFAData) == 53006),] LSkewUrb(0.124, 0.1138, DeUrb = TRUE) #Get the pooled Lmoment ratios results for catchment 53006 and apply the urban #Get the CDS & adjustment using the pooled LSkew, and the URBEXT2000 for site 53006. CDs.53006 <- GetCDs(53006) QuickResults(CDs.53006)[[2]] LSkewUrb(0.194, 0.1138)#Choose an urban site (site 53006) from the NRFA data then apply a de-urban #adjustment using the Lcv and URBEXT2000 displayed NRFAData[which(rownames(NRFAData) == 53006),] LSkewUrb(0.124, 0.1138, DeUrb = TRUE) #Get the pooled Lmoment ratios results for catchment 53006 and apply the urban #Get the CDS & adjustment using the pooled LSkew, and the URBEXT2000 for site 53006. CDs.53006 <- GetCDs(53006) QuickResults(CDs.53006)[[2]] LSkewUrb(0.194, 0.1138)
Derives monthly statistics from a data.frame with Dates or POSIXct in the first column and variable of interest in the second
MonthlyStats( x, Stat, AggStat = NULL, TS = FALSE, Plot = FALSE, ylab = "Magnitude", main = "Monthly Statistics", col = "grey" )MonthlyStats( x, Stat, AggStat = NULL, TS = FALSE, Plot = FALSE, ylab = "Magnitude", main = "Monthly Statistics", col = "grey" )
x |
a data.frame with Dates or POSIXct in the first column and numeric vector in the second. |
Stat |
A user chosen function to calculate the statistic of interest; mean or sum for example. Could be a user developed function. |
AggStat |
the aggregating statistic. The default is mean. The function applied must have an na.rm argument (base R stat functions such as mean, max, and sum all have an na.rm argument.). |
TS |
A logical statement with a default of FALSE. If TRUE, instead of a dataframe of monthly statistics and average statistics, a monthly time series is returned. |
Plot |
logical argument with a default of FALSE. If TRUE the monthly statistics are plotted. |
ylab |
A label for the y axis of the plot. The default is "Magnitude" |
main |
A title for the plot. The default is "Monthly Statistics" |
col |
A choice of colour for the bar plot. A single colour or a vector (a colour for each bar). |
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.
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.
Anthony Hammond
# Get the mean flows for each month for the Thames at Kingston QMonThames <- MonthlyStats(ThamesPQ[,c(1,3)], Stat = mean, ylab = "Discharge (m3/s)", main = "Thames at Kingston monthly mean flow", Plot = TRUE) # Get the monthly sums of rainfall for the Thames at Kingston PMonThames <- MonthlyStats(ThamesPQ[,c(1,2)], Stat = sum, ylab = "Rainfall (mm)", main = "Thames as Kingston monthly rainfall", Plot = TRUE)# Get the mean flows for each month for the Thames at Kingston QMonThames <- MonthlyStats(ThamesPQ[,c(1,3)], Stat = mean, ylab = "Discharge (m3/s)", main = "Thames at Kingston monthly mean flow", Plot = TRUE) # Get the monthly sums of rainfall for the Thames at Kingston PMonThames <- MonthlyStats(ThamesPQ[,c(1,2)], Stat = sum, ylab = "Rainfall (mm)", main = "Thames as Kingston monthly rainfall", Plot = TRUE)
Calculates the euclidean distance between two british national grid reference points using the pythagorean method
NGRDist(i, j)NGRDist(i, j)
i |
a numeric vector of length two. The first being the easting and the second being the northing of the first site |
j |
a numeric vector of length two. The first being the easting and the second being the northing of the second site |
Note, that the result is converted to km when six digits are used for easting and northing, when six digits would usually provide a result in metres.
A distance in kilometres (if six digits for easting and northing are used)
Anthony Hammond
#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 view the eastings and northings GetCDs(10001) GetCDs(39001) NGRDist(i = c(381355, 839183), j = c(462899, 187850))#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 view the eastings and northings GetCDs(10001) GetCDs(39001) NGRDist(i = c(381355, 839183), j = c(462899, 187850))
Adjusts the linear coefficient of variation (Lcv) and the linear skewness to account for non-flood years
NonFloodAdj(x)NonFloodAdj(x)
x |
The annual maximum sample. Numeric vector |
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.
A list is returned. The first element of the list is a dataframe with one row and two columns. 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.
Anthony Hammond
# Get an anuual maximum sample with a BFIHOST above 0.65 and with some # annual maximums lower than median(AM)/2. And then apply the function. NonFloodAdj(GetAM(44013)[,2])# Get an anuual maximum sample with a BFIHOST above 0.65 and with some # annual maximums lower than median(AM)/2. And then apply the function. NonFloodAdj(GetAM(44013)[,2])
Applies the NonFloodAdj function to adjust the LCV and LSKEW of one or more sites in a pooling group.
NonFloodAdjPool(x, Index = NULL, AutoP = NULL, ReturnStats = FALSE)NonFloodAdjPool(x, Index = NULL, AutoP = NULL, ReturnStats = FALSE)
x |
A pooling group, derived from the Pool() or PoolSmall() functions. |
Index |
An vector of indices (row numbers) of sites to be adjusted. If Index = NULL (the default) the function is applied to all sites. |
AutoP |
A percentage (numeric) of non flood years. Any sites in the group exceeding this value will be adjusted. This is an automated approach so that the user doesn't need to specify Index. If no sites are above AutoP, the function is applied to all sites. |
ReturnStats |
Logical with a default of FALSE. If set to TRUE, a dataframe of non-flood year stats is returned (see 'Value' section below) instead of the adjusted Pooling group. |
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.
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 forst the number of non-flood years, the second is the number of years, and the third is the associated percentage.
Anthony Hammond
# Set up a pooling group for site 44013. Then apply the function. Pool44013 <- Pool(GetCDs(44013)) PoolNF <- NonFloodAdjPool(Pool44013) #return the non flood stats for the pooling group NonFloodAdjPool(Pool44013, ReturnStats = TRUE)# Set up a pooling group for site 44013. Then apply the function. Pool44013 <- Pool(GetCDs(44013)) PoolNF <- NonFloodAdjPool(Pool44013) #return the non flood stats for the pooling group NonFloodAdjPool(Pool44013, ReturnStats = TRUE)
A data.frame of catchment descriptors, Lmoments, Lmoment ratios, sample size and median annual maximum flow (QMED). NRFA Peak Flow Dataset - Version 13.0.2.
NRFADataNRFAData
A data frame with 543 rows and 27 variables
The functions for pooling group formation and estimation rely on this dataframe. 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. Although, usually, in the latter case, such changes will be more appropriately applied to the formed pooling group. 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.
https://nrfa.ceh.ac.uk/peak-flow-dataset
Estimates the parameters of the Generalised extreme value, generalised logistic, Kappa3, or Gumbel distribution from known return period estimates
OptimPars(x, dist = "GenLog")OptimPars(x, dist = "GenLog")
x |
a data.frame with RPs in the first column and associated estimates in the second column |
dist |
a choice of distribution for the estimates. The choices are "GenLog", "GEV", "Kappa3", or "Gumbel" - the generalised logistic, generalised extreme value, Kappa3, and Gumbel distribution, respectively. The default is "GenLog" |
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.
The parameters of one of four user chosen distributions; Generalised logistic, generalised extreme value, Gumbel, and Kappa3.
Anthony Hammond
#Get some catchment descriptors and some quick results. Then estmate the GenLog parameters Results <- QuickResults(GetCDs(27051), plot = FALSE)[[1]] OptimPars(Results[,1:2])#Get some catchment descriptors and some quick results. Then estmate the GenLog parameters Results <- QuickResults(GetCDs(27051), plot = FALSE)[[1]] OptimPars(Results[,1:2])
Function to develop a pooling group based on catchment descriptors
Pool( CDs = NULL, AREA, SAAR, FARL, FPEXT, N = 500, exclude = NULL, iug = FALSE, UrbMax = 0.03, DeUrb = FALSE )Pool( CDs = NULL, AREA, SAAR, FARL, FPEXT, N = 500, exclude = NULL, iug = FALSE, UrbMax = 0.03, DeUrb = FALSE )
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
AREA |
catchment area in km2 |
SAAR |
catchment standard average annual rainfall (1961-1990) in mm |
FARL |
catchment flood attenuation from reservoirs & lakes |
FPEXT |
catchment floodplain extent. The proportion of the catchment that is estimated to be inundated by a 100-year flood |
N |
minimum Number of total gauged record years for the pooling group |
exclude |
sites to exclude from the pooling group. Either a single site reference or a vector of site references (numeric) |
iug |
iug stands for 'include urban gauge' - which refers to a gauged subject site if it's > UrbMax. It's a logical argument with default of FALSE. TRUE will over-ride the default and add the closest site in catchment descriptor space (should be the gauge of interest) to the pooling group if it has URBEXT2000 >= UrbMax |
UrbMax |
Maximum URBEXT2000 level with a default of 0.03. Any catchment with URBEXT2000 above this level will be excluded from the pooling group |
DeUrb |
logical argument with a default of FALSE. If true, the Lcv and LSkew of any site in the pooling group with URBEXT2000 > 0.03 will be de-urbanised |
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). Sites with URBEXT2000 (urban extent) > 0.03 are excluded by default and this can be adjusted with UrbMax. If a gauged assessment is required and the site of interest is > UrbMax it can be included by setting iug = TRUE. De-urbanise the Lcv and Lskew (L-moment ratios) for sites with URBEXT2000 > UrbMax by setting DeUrb = TRUE. 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. An example of doing so is provided below. The pooling method is outlined in Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation.
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.
Anthony Hammond
#Get some catchment descriptors CDs.73005 <- GetCDs(73005) #Set up a pooling group object called Pool.73005 excluding sites 79005 & 71011. #Then print the group to the console Pool.73005 <- Pool(CDs.73005, exclude = c(79005, 71011)) Pool.73005 #Form a pooling group, called PoolGroup, with the catchment descriptors specifically PoolGroup <- Pool(AREA = 1000, SAAR = 800, FARL = 1, FPEXT = 0.01) #Form a pooling group using an urban catchment which is intended for enhanced #single site estimation - by including it in the group. CDs.39001 <- GetCDs(39001) Pool.39001 <- Pool(CDs.39001, iug = TRUE, DeUrb = TRUE) #Change the Lcv and LSkew of the top site in the pooling group to 0.19 & 0.18, #respectively. PoolUpdate <- LRatioChange(Pool.39001, SiteID = 39001, 0.19, 0.18)#Get some catchment descriptors CDs.73005 <- GetCDs(73005) #Set up a pooling group object called Pool.73005 excluding sites 79005 & 71011. #Then print the group to the console Pool.73005 <- Pool(CDs.73005, exclude = c(79005, 71011)) Pool.73005 #Form a pooling group, called PoolGroup, with the catchment descriptors specifically PoolGroup <- Pool(AREA = 1000, SAAR = 800, FARL = 1, FPEXT = 0.01) #Form a pooling group using an urban catchment which is intended for enhanced #single site estimation - by including it in the group. CDs.39001 <- GetCDs(39001) Pool.39001 <- Pool(CDs.39001, iug = TRUE, DeUrb = TRUE) #Change the Lcv and LSkew of the top site in the pooling group to 0.19 & 0.18, #respectively. PoolUpdate <- LRatioChange(Pool.39001, SiteID = 39001, 0.19, 0.18)
Provides pooled results from a pooling group - gauged, ungauged and with urban adjustment if necessary.
PoolEst( x, gauged = FALSE, QMED, dist = "GenLog", RP = c(2, 5, 10, 20, 50, 75, 100, 200, 500, 1000), UrbAdj = FALSE, CDs = NULL, URBEXT = NULL, fseQMED = 1.46 )PoolEst( x, gauged = FALSE, QMED, dist = "GenLog", RP = c(2, 5, 10, 20, 50, 75, 100, 200, 500, 1000), UrbAdj = FALSE, CDs = NULL, URBEXT = NULL, fseQMED = 1.46 )
x |
pooling group derived from the Pool function |
gauged |
logical argument with a default of FALSE. TRUE for gauged results and FALSE for ungauged |
QMED |
estimate of the median annual maximum flow |
dist |
a choice of distribution for the estimates. The choices are "GenLog", "GEV", "Kappa3", or "Gumbel"; the generalised logistic, generalised extreme value, Kappa3, and Gumbel distribution, respectively. The default is "GenLog" |
RP |
return period of interest. By default the following RPs are provided: 2, 5, 10, 20, 50, 75, 100, 200, 500, 1000 |
UrbAdj |
logical argument with a default of FALSE. When TRUE, an urban adjustment is applied to the pooled Lcv and LSkew |
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
URBEXT |
the catchment URBEXT2000, to be supplied if UrbAdj is TRUE and if CDs have not been |
fseQMED |
factorial standard error of the median annual maximum (QMED) estimate, used for quantifying ungauged uncertainty. Default is 1.46 |
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. If the catchment of interest is urban, the UrbAdj argument can be set to TRUE. If this is done, either URBEXT (urban extent) needs to be provided or the catchment descriptors, derived from CDsXML or GetCDs. The methods for estimating pooled growth curves are according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation. The methods for estimating the L-moments and growth factors are outlined in the Flood Estimation Handbook (1999), volume 3. The methods for quantifying uncertainty are detailed in Hammond, A. (2022). Easy methods for quantifying the uncertainty of FEH pooling analysis. Circulation - The Newsletter of the British Hydrological Society (152). When UrbAdj = TRUE, urban adjustment is applied to the QMED estimate according to the method outlined in the guidance by Wallingford HydroSolutions: 'WINFAP 4 Urban Adjustment Procedures'.
If RP is default then a list of length 4. Element one is a data frame with columns; return period (a range from 2 - 1000), peak flow estimates (Q), growth factor estimates (GF), lower and upper intervals of uncertainty (68 percent intervals for ungauged and 95 percent for gauged). The second element is the estimated Lcv and Lskew. The third provides distribution parameters for the growth curve. The fourth provides distribution parameters for the frequency curve. If RP is not the default only the first two elements are returned.
Anthony Hammond
#Get some catchment descriptors and form a pooling group. It's urban and #therefore the site of interest is not included. CDs.27083 <- GetCDs(27083) Pool.27083 <- Pool(CDs.27083) #Get results for the ungauged case, with urban adjustment PoolEst(Pool.27083, QMED = 12, UrbAdj = TRUE, CDs = CDs.27083) #Form the group again with the urban gauge included & undertake a gauged estimate #with urban adjustment. QMED in this example is estimated as the median of the annual #maximum series for site 27083. PoolG.27083 <- PoolG.27083 <- Pool(CDs.27083, iug = TRUE, DeUrb = TRUE) PoolEst(PoolG.27083, QMED = 12.5, UrbAdj = TRUE, CDs = CDs.27083)#Get some catchment descriptors and form a pooling group. It's urban and #therefore the site of interest is not included. CDs.27083 <- GetCDs(27083) Pool.27083 <- Pool(CDs.27083) #Get results for the ungauged case, with urban adjustment PoolEst(Pool.27083, QMED = 12, UrbAdj = TRUE, CDs = CDs.27083) #Form the group again with the urban gauge included & undertake a gauged estimate #with urban adjustment. QMED in this example is estimated as the median of the annual #maximum series for site 27083. PoolG.27083 <- PoolG.27083 <- Pool(CDs.27083, iug = TRUE, DeUrb = TRUE) PoolEst(PoolG.27083, QMED = 12.5, UrbAdj = TRUE, CDs = CDs.27083)
Function to develop a small catchments pooling group based on catchment descriptors
PoolSmall( CDs = NULL, AREA, SAAR, N = 500, exclude = NULL, iug = FALSE, UrbMax = 0.03, DeUrb = FALSE )PoolSmall( CDs = NULL, AREA, SAAR, N = 500, exclude = NULL, iug = FALSE, UrbMax = 0.03, DeUrb = FALSE )
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
AREA |
catchment area in km2 |
SAAR |
catchment standard average annual rainfall (1961-1990) in mm |
N |
minimum Number of total gauged record years for the pooling group |
exclude |
sites to exclude from the pooling group. Either a single site reference or a vector of site references (numeric) |
iug |
iug stands for 'include urban gauge' - which refers to a gauged subject site if it's > UrbMax. It's a logical argument with default of FALSE. TRUE will over-ride the default and add the closest site in catchment descriptor space (should be the gauge of interest) to the pooling group if it has URBEXT2000 >= UrbMax |
UrbMax |
Maximum URBEXT2000 level with a default of 0.03. Any catchment with URBEXT2000 above this level will be excluded from the pooling group |
DeUrb |
logical argument with a default of FALSE. If true, the Lcv and LSkew of any site in the pooling group with URBEXT2000 > 0.03 will be de-urbanised |
A pooling group is created from a CDs object, derived from GetCDs or CDsXML, or specifically with the necessary 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). Sites with URBEXT2000 (urban extent) > 0.03 are excluded by default and this can be adjusted with the UrbMax argument. If a gauged assessment is required and the site of interest is > UrbMax it can be included by setting iug = TRUE. De-urbanise the Lcv and Lskew (L-moment ratios) of sites with URBEXT2000 > 0.03 by setting DeUrb = TRUE. 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.
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.
Anthony Hammond
#Get some catchment descriptors CDs.21001 <- GetCDs(21001) #Set up a pooling group object called Pool.21001 excluding site 206006 #Then print the group to the console Pool.21001 <- PoolSmall(CDs.21001, exclude = 206006) Pool.21001 #Form a pooling group, called PoolGroup, with the catchment descriptors specifically PoolGroup <- PoolSmall(AREA = 22, SAAR = 1702)#Get some catchment descriptors CDs.21001 <- GetCDs(21001) #Set up a pooling group object called Pool.21001 excluding site 206006 #Then print the group to the console Pool.21001 <- PoolSmall(CDs.21001, exclude = 206006) Pool.21001 #Form a pooling group, called PoolGroup, with the catchment descriptors specifically PoolGroup <- PoolSmall(AREA = 22, SAAR = 1702)
Extracts independent peaks over a threshold from a sample
POTextract( x, div = NULL, TimeDiv = NULL, thresh = 0.975, Plot = TRUE, ylab = "Magnitude", xlab = "Time", main = "Peaks over threshold" )POTextract( x, div = NULL, TimeDiv = NULL, thresh = 0.975, Plot = TRUE, ylab = "Magnitude", xlab = "Time", main = "Peaks over threshold" )
x |
either a numeric vector or dataframe with date (or POSIXct) in the first column and hydrological variable in the second |
div |
numeric percentile (between 0 and thresh), either side of which two peaks over the threshold are considered independent. Default is the mean of the sample. |
TimeDiv |
Number of timesteps to define independence (supplements the div argument). As a default this is NULL and only 'div' defines independence. Currently this is only applicable for data.frames. |
thresh |
user chosen threshold. Default is 0.975 |
Plot |
logical argument with a default of TRUE. When TRUE, the full hydrograph with the peaks over the threshold highlighted is plotted |
ylab |
Label for the plot yaxis. Default is "Magnitude" |
xlab |
Label (character) for the plot x axis. Default is "Time". |
main |
Title for the plot. Default is "Peaks over threshold" |
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. In which case, with hourly data for example, TimeDiv could be set to 120 to ensure each peak is separated by five days either side as well as at least one hour with 0 rainfall. When plotted, the blue line is the threshold, and the green line is the independence line (div).
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.
Anthony Hammond
#Extract POT data from Thames mean daily flow 2000-10-01 to 2015-09-30 with #div = mean and threshold = 0.95. Then display the first six rows ThamesQPOT <- POTextract(ThamesPQ[, c(1,3)], thresh = 0.9) head(ThamesQPOT) #Extract Thames POT from only the numeric vector of flows and display the #first six rows ThamesQPOT <- POTextract(ThamesPQ[, 3], thresh = 0.9) head(ThamesQPOT) #Extract the Thames POT precipitation with a div of 0, the default #threshold, and 5 timesteps (days) either side of the peak. Then display the first six rows ThamesPPOT <- POTextract(ThamesPQ[, c(1,2)], div = 0, TimeDiv = 5) head(ThamesPPOT)#Extract POT data from Thames mean daily flow 2000-10-01 to 2015-09-30 with #div = mean and threshold = 0.95. Then display the first six rows ThamesQPOT <- POTextract(ThamesPQ[, c(1,3)], thresh = 0.9) head(ThamesQPOT) #Extract Thames POT from only the numeric vector of flows and display the #first six rows ThamesQPOT <- POTextract(ThamesPQ[, 3], thresh = 0.9) head(ThamesQPOT) #Extract the Thames POT precipitation with a div of 0, the default #threshold, and 5 timesteps (days) either side of the peak. Then display the first six rows ThamesPPOT <- POTextract(ThamesPQ[, c(1,2)], div = 0, TimeDiv = 5) head(ThamesPPOT)
Extracts independent peaks over a threshold from a sample, using time as the independence criteria.
POTt( x, threshold = 0.975, div, Plot = TRUE, PlotType = "l", main = "Peaks over threhsold", ylab = "Magnitude", xlab = "Time" )POTt( x, threshold = 0.975, div, Plot = TRUE, PlotType = "l", main = "Peaks over threhsold", ylab = "Magnitude", xlab = "Time" )
x |
either a numeric vector or dataframe with date (or POSIXct) in the first column and hydrological variable in the second |
threshold |
user chosen threshold. Default is 0.975 |
div |
number of time steps between peaks to ensure independence. |
Plot |
logical argument with a default of TRUE. When TRUE, the full hydrograph with the peaks over the threshold highlighted is plotted |
PlotType |
Type of plot with a default of "l" for line graph. For rainfall type "h" for bars could be used. |
main |
Title for the plot. Default is "Peaks over threshold" |
ylab |
Label (character) for the plot y axis. Default is "Magnitude" |
xlab |
Label (character) for the plot x axis. Default is "Time". |
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.
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.
Anthony Hammond
#Extract POT data from Thames catchment daily rainfall 2000-10-01 to 2015-09-30 with #div = 14 (14 days) and threshold = 0.975. Then display the first six rows ThamesPPOT <- POTt(ThamesPQ[, c(1,2)], div = 14) head(ThamesPPOT) #Extract Thames rainfall POT from only the numeric vector of rainfall, with threshold #set to 0.95 and div set to 14. Then display the first six rows ThamesPPOT <- POTt(ThamesPQ[, 2], threshold = 0.95, div = 14) head(ThamesPPOT)#Extract POT data from Thames catchment daily rainfall 2000-10-01 to 2015-09-30 with #div = 14 (14 days) and threshold = 0.975. Then display the first six rows ThamesPPOT <- POTt(ThamesPQ[, c(1,2)], div = 14) head(ThamesPPOT) #Extract Thames rainfall POT from only the numeric vector of rainfall, with threshold #set to 0.95 and div set to 14. Then display the first six rows ThamesPPOT <- POTt(ThamesPQ[, 2], threshold = 0.95, div = 14) head(ThamesPPOT)
Estimated median annual maximum flow from catchment descriptors and donor sites
QMED( CDs = NULL, Don1 = NULL, Don2 = NULL, UrbAdj = FALSE, uef = FALSE, DonUrbAdj = FALSE, AREA, SAAR, FARL, BFIHOST, URBEXT2000 = NULL, Easting = NULL, Northing = NULL )QMED( CDs = NULL, Don1 = NULL, Don2 = NULL, UrbAdj = FALSE, uef = FALSE, DonUrbAdj = FALSE, AREA, SAAR, FARL, BFIHOST, URBEXT2000 = NULL, Easting = NULL, Northing = NULL )
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
Don1 |
numeric site reference for the a single donor (for donor candidates see DonAdj function) |
Don2 |
vector of two site references for two donors (for donor candidates see DonAdj function) |
UrbAdj |
logical argument with a default of FALSE. True applies an urban adjustment |
uef |
logical argument with a default of FALSE. If true an urban expansion factor is applied to the URBEXT2000 value - using the current year. |
DonUrbAdj |
logical argument with a default of FALSE. If TRUE, an urban adjustement is applied to the donor/s QMEDcds estimate. |
AREA |
catchment area in km2 |
SAAR |
standard average annual rainfall (mm) |
FARL |
flood attenuation from reservoirs and lakes |
BFIHOST |
baseflow index calculated from the catchment hydrology of soil type classification |
URBEXT2000 |
measure of catchment urbanisation |
Easting |
Easting. A six digit Easting (British national grid reference). |
Northing |
Northing. A six digit Northing (British national grid reference). |
QMED is estimated from catchment descriptors: QMED = 8.3062*AREA^0.8510 0.1536^(1000/SAAR) FARL^3.4451 0.0460^(BFIHOST^2) as derived in Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation. The single donor method is from the same paper. The method for two donors is outlined in 'Kjeldsen, T. (2019). Adjustment of QMED in ungauged catchments using two donor sites. Circulation - The Newsletter of the British Hydrological Society, 4'. When UrbAdj = TRUE, urban adjustment is applied to the QMED estimate according to the method outlined in the guidance by Wallingford HydroSolutions: 'WINFAP 4 Urban Adjustment Procedures'. Urban donors should be avoided, but in the case that the subject catchment is rural, and the donor is urban, the QMEDcd estimate of the donor (or donors) can be urban adjusted by setting the DonUrbAdj argument to TRUE. For flexibility there is the option to input the relevant catchment descriptors directly rather than a CDs object.
An estimate of QMED from catchment descriptors. If two donors are used the associated weights are also returned
Anthony Hammond
#Get some catchment descriptors and calculate QMED as if it was ungauged, with #no donors, one donor, and two donors CDs.55004 <- GetCDs(55004) QMED(CDs.55004) QMED(CDs.55004, Don1 = 55012) QMED(CDs.55004, Don2 = c(55012, 60007)) #Get CDs for urban gauge and calculate QMED with urban adjustment CDs.27083 <- GetCDs(27083) QMED(CDs.27083, UrbAdj = TRUE)#Get some catchment descriptors and calculate QMED as if it was ungauged, with #no donors, one donor, and two donors CDs.55004 <- GetCDs(55004) QMED(CDs.55004) QMED(CDs.55004, Don1 = 55012) QMED(CDs.55004, Don2 = c(55012, 60007)) #Get CDs for urban gauge and calculate QMED with urban adjustment CDs.27083 <- GetCDs(27083) QMED(CDs.27083, UrbAdj = TRUE)
A data.frame of catchment & data descriptors relating to the median annual maximum flow (QMED). NRFA Peak Flow Dataset - Version 13.0.2
QMEDDataQMEDData
A data frame with 897 rows and 26 variables
The functions for QMED estimation and retrieval of catchment descriptors rely on this dataframe. 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.
https://nrfa.ceh.ac.uk/peak-flow-dataset
Applies a donor adjustment to the median annual maximum flow (QMED) estimate
QMEDDonEq( AREA, SAAR, FARL, BFIHOST, QMEDgObs, QMEDgCds, xSI, ySI, xDon, yDon, alpha = TRUE )QMEDDonEq( AREA, SAAR, FARL, BFIHOST, QMEDgObs, QMEDgCds, xSI, ySI, xDon, yDon, alpha = TRUE )
AREA |
catchment area in km2 |
SAAR |
standardised average annual rainfall in mm |
FARL |
flood attenuation from reservoirs and lakes |
BFIHOST |
the baseflow index as a function of soil type |
QMEDgObs |
the observed QMED at the donor site |
QMEDgCds |
the QMED equation derived QMED at the donor site |
xSI |
the catchment centroid easting for the site of interest |
ySI |
the catchment centroid northing for the site of interest |
xDon |
the catchment centroid easting for the donor site |
yDon |
the catchment centroid northing for the donor site |
alpha |
a logical argument with a default of TRUE. When FALSE the exponent in the donor equation is set to one. Otherwise it is determined by the distance between the donor and the subject site |
Although a single donor adjustment can be applied with the DonAdj() function and the QMED(), this is provided for flexibility. The method is that of Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation (2008).
Anthony Hammond
#Get observed QMED for site 15006 Qob <- median(GetAM(15006)[,2]) #Get QMED equation estimated QMED for the donor site QCD <- QMED(CDs = GetCDs(15006)) #display CDs for site 27051 & note the easting and northing GetCDs(27051) #display CDs for site 15006 & note the easting and northing GetCDs(15006) #Apply the QMEDDonEq function with the information gained QMEDDonEq(194, 1096, 0.955, 0.297, Qob, QCD, xSI = 289289,ySI = 947523,xDon = 280908,yDon = 953653)#Get observed QMED for site 15006 Qob <- median(GetAM(15006)[,2]) #Get QMED equation estimated QMED for the donor site QCD <- QMED(CDs = GetCDs(15006)) #display CDs for site 27051 & note the easting and northing GetCDs(27051) #display CDs for site 15006 & note the easting and northing GetCDs(15006) #Apply the QMEDDonEq function with the information gained QMEDDonEq(194, 1096, 0.955, 0.297, Qob, QCD, xSI = 289289,ySI = 947523,xDon = 280908,yDon = 953653)
Estimates the median annual maximum flow (QMED) factorial standard error (FSE) by bootstrapping the sample
QMEDfseSS(x)QMEDfseSS(x)
x |
a numeric vector. The sample of interest |
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.
The factorial standard error for the median of a sample.
Anthony Hammond
#Extract an AMAX sample and estimate the QMED factorial standard error AM.203018 <- GetAM(203018) QMEDfseSS(AM.203018$Flow)#Extract an AMAX sample and estimate the QMED factorial standard error AM.203018 <- GetAM(203018) QMEDfseSS(AM.203018$Flow)
Estimates the median annual maximum flow (QMED) from non-flood flows
QMEDLink(Q5dmf, Q10dmf, DPSBAR, BFI)QMEDLink(Q5dmf, Q10dmf, DPSBAR, BFI)
Q5dmf |
numeric. The daily mean flow that is exceeded 5 percent of the time |
Q10dmf |
numeric. The daily mean flow that is exceeded 10 percent of the time |
DPSBAR |
a catchment descriptor. The average drainage path slope of the catchment |
BFI |
the baseflow index of the gauged flow |
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 desciptor; 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.
Anthony Hammond
#Calculate the QMED for site 1001 (Wick at Tarroul) QMEDLink(10.14, 7.352, 29.90, 0.39)#Calculate the QMED for site 1001 (Wick at Tarroul) QMEDLink(10.14, 7.352, 29.90, 0.39)
Estimates the median annual maximum flow (QMED) from peaks over threshold data
QMEDPOT(x, ppy)QMEDPOT(x, ppy)
x |
numerical vector. POT data |
ppy |
number of peaks per year in the POT data |
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).
Anthony Hammond
#Extract some POT data and estimate QMED ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) QMEDPOT(ThamesPOT$peak, ppy = 1.867263)#Extract some POT data and estimate QMED ThamesPOT <- POTextract(ThamesPQ[,c(1,3)], thresh = 0.90) QMEDPOT(ThamesPOT$peak, ppy = 1.867263)
Provides pooled gauged, ungauged, or fake ungauged results, directly from the catchment descriptors
QuickResults( CDs, gauged = FALSE, dons = 2, Qmed = NULL, FUngauged = FALSE, plot = TRUE, dist = "GenLog" )QuickResults( CDs, gauged = FALSE, dons = 2, Qmed = NULL, FUngauged = FALSE, plot = TRUE, dist = "GenLog" )
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
gauged |
logical argument with a default of FALSE. TRUE for gauged results and FALSE for ungauged |
dons |
number of donors required with a choice of 0, 1, or 2 |
Qmed |
user supplied QMED which overrides the default QMED estimate |
FUngauged |
logical argument with a default of FALSE. TRUE provides an ungauged estimate whilst excluding the gauged site (the site with the most similar CDs) |
plot |
logical argument with a default of TRUE. TRUE provides an extreme value plot. FALSE prevents the plot |
dist |
a choice of distribution for the estimates. The choices are "GenLog", "GEV", "Kappa3", or "Gumbel; the generalised logistic, generalised extreme value, Kappa3,and Gumbel distributions, respectively. The default is "GenLog" |
The quick results function provides results with a default pooling group. If gauged = FALSE the median annual maximum flood (QMED) is estimated from catchment descriptors using the QMED equation and then adjusted with two of the closest un-urban gauged sites (can be changed to 0 or 1 donors). If the site is urban, an urban adjustment is made to the QMED and to the pooled growth curve. If gauged = TRUE QMED is the median of the gauged annual maxima and the growth curve is formed with the gauged weighting procedure (often known as enhanced single site). If the gauged catchment is urban, it's included in the pooling group and deurbanised before an urban adjustment is made to the final growth curve. If FUngauged = TRUE, the top site in the pooling group is excluded and the estimate is performed henceforth in the manner of gauged = FALSE. If the CDs are from a gauged site that is not in the list of sites that are considered suitable for pooling, it won't be included in the pooling group. In which case, if gauged = TRUE, the result will be erroneous.
A list of length two. Element one is a data frame with columns; return period (RP), peak flow estimates (Q) and growth factor estimates (GF). Two additional columns quantify the uncertainty. The second element is the estimated Lcv and Lskew (linear coefficient of variation and skewness). By default an extreme value plot is also returned
Anthony Hammond
#Get some catchment descriptors CDs.73005 <- GetCDs(73005) #Get default ungauged results QuickResults(CDs.73005) #Get gauged results with a GEV distribution QuickResults(CDs.73005, gauged = TRUE, dist = "GEV") #Get fake ungauged results with one donor QuickResults(CDs.73005, FUngauged = TRUE, dons = 1)#Get some catchment descriptors CDs.73005 <- GetCDs(73005) #Get default ungauged results QuickResults(CDs.73005) #Get gauged results with a GEV distribution QuickResults(CDs.73005, gauged = TRUE, dist = "GEV") #Get fake ungauged results with one donor QuickResults(CDs.73005, FUngauged = TRUE, dons = 1)
Optimises a power law rating equation from observed discharge and stage
Rating(x, a = NULL)Rating(x, a = NULL)
x |
a data.frame with discharge in the first column and stage in the second |
a |
a user defined stage correction |
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.
A list with three elements. The first is a vetor 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.
Anthony Hammond
# Make up Some data: Q <- 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(Q, Stage) #apply the rating function: Rating(Observations) #Hold the stage correction at zero Rating(Observations, a = 0)# Make up Some data: Q <- 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(Q, Stage) #apply the rating function: Rating(Observations) #Hold the stage correction at zero Rating(Observations, a = 0)
Provides outputs of the ReFH model from catchment descriptors or user defined inputs
ReFH( CDs = NULL, Depth = NULL, duration = NULL, timestep = NULL, scaled = NULL, PlotTitle = NULL, RPa = NULL, alpha = TRUE, season = NULL, AREA = NULL, TP = NULL, BR = NULL, BL = NULL, Cmax = NULL, Cini = NULL, BFini = NULL, Rain = NULL )ReFH( CDs = NULL, Depth = NULL, duration = NULL, timestep = NULL, scaled = NULL, PlotTitle = NULL, RPa = NULL, alpha = TRUE, season = NULL, AREA = NULL, TP = NULL, BR = NULL, BL = NULL, Cmax = NULL, Cini = NULL, BFini = NULL, Rain = NULL )
CDs |
catchment descriptors derived from either GetCDs or ImportCD |
Depth |
a numeric value. The depth of rainfall used as input in the estimation of a design hydrograph. The default, when Depth = NULL, is a two year rainfall. |
duration |
a numeric value. A duration (hrs) for the design rainfall |
timestep |
a numeric value. A user defined data interval. The default changes depending on the estimated time to peak to formulate a sensible looking result |
scaled |
a numeric value of peak flow in m3/s |
PlotTitle |
a character string. A user defined title for the ReFH plot |
RPa |
return period for alpha adjustment. This is only for the purposes of the alpha adjustment, it doesn't change the rainfall input |
alpha |
a logical argument with default TRUE. If TRUE the alpha adjustment is applied based on RPa. If FALSE, no alpha adjustment is made |
season |
a choice of "summer" or "winter". The default is "summer" in urban catchments (URBEXT2000 > 0.03) and "winter" in rural catchments |
AREA |
numeric. Catchment area in km2. |
TP |
numeric. Time to peak parameter (hours) |
BR |
numeric. Baseflow recharge parameter |
BL |
numeric. Baseflow lag parameter (hours) |
Cmax |
numeric. Maximum soil moisture capacity parameter (mm) |
Cini |
numeric. Initial soil moisture content (mm) |
BFini |
numeric. Initial baseflow (m3/s) |
Rain |
numeric. User input rainfall (hourly). A numeric vector |
The ReFH is described in the Flood Estimation Handbook Supplementary Report No.1 (2007). The method to derive design rainfall profiles is described in the Flood Estimation Handbook (1999), volume 2. 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 the two year rainfall for the critical duration. The latter is based on a quadratic interpolation of the catchment descriptors RMED1H, RMED1D, and RMED2D (then a seasonal correction factor is applied). Parameters and initial conditions can also be individually input by the user. If a parameter argument is used for one or more of the parameters, then these overwrite the CD derived parameters. If a value for the scaled argument is provided (m3/s), a scaled hydrograph is returned. The RPa argument doesn't change the rainfall input and is only needed for the alpha adjustment (see the FEH supplement report no.1).
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.
Anthony Hammond
#Get CDs and apply the ReFH function CDs.203018 <- GetCDs(203018) ReFH(CDs.203018) #Apply the ReFH function, scale to a 100-year flow estimate and change the plot title accordingly ReFH(CDs.203018, scaled = 182, PlotTitle = "100-Year Design Hydrograph - Site 203018") #Apply the ReFH function with a user defined initial baseflow ReFH(CDs.203018, BFini = 6)#Get CDs and apply the ReFH function CDs.203018 <- GetCDs(203018) ReFH(CDs.203018) #Apply the ReFH function, scale to a 100-year flow estimate and change the plot title accordingly ReFH(CDs.203018, scaled = 182, PlotTitle = "100-Year Design Hydrograph - Site 203018") #Apply the ReFH function with a user defined initial baseflow ReFH(CDs.203018, BFini = 6)
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
SCF(SAAR, duration)SCF(SAAR, duration)
SAAR |
standardised average annual rainfall. Numeric |
duration |
duration in hours. Numeric |
The SCF and it's 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.
A data.frame of one row and two columns: SCFSummer and SCFWinter.
Anthony Hammond
#Derive the SCFs for a SAAR of 1981 and a duration of 6.5 SCF(1981, 6.5)#Derive the SCFs for a SAAR of 1981 and a duration of 6.5 SCF(1981, 6.5)
Simulation of a random sample from the generalised extreme value, generalised logistic, Gumbel, Kappa3, or generalised Pareto distributions
SimData(n, pars = NULL, dist = "GenLog", GF = NULL)SimData(n, pars = NULL, dist = "GenLog", GF = NULL)
n |
sample size to be simulated |
pars |
vector of parameters in the order of location, scale, shape (only location and shape for Gumbel) |
dist |
choice of distribution. Either "GEV", "GenLog", "Gumbel", "Kappa3", or "GenPareto" |
GF |
vector of GF inputs in the order of Lcv, LSkew, QMED (only Lcv and QMED if dist = "Gumbel") |
The simulated sample can be generated using distribution parameters, or the growth factor (GF) inputs; linear coefficient of variationn (Lcv), linear skewness (LSkew) & the median annual maximum (QMED).
A random sample of size n for the chosen distribution.
Anthony Hammond
#Simulate a sample of size 30 using parameters GenLog and 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")#Simulate a sample of size 30 using parameters GenLog and 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")
A data.frame of four columns; Date, Precipitation (P), & daily mean flow (Q)
ThamesPQThamesPQ
A data frame with 5478 rows and 4 columns:
Date
Precipitation, in mm
Daily mean discharge, in m3/s
https://nrfa.ceh.ac.uk/data/station/meanflow/39001
A hypothesis test for trend
TrendTest(x, method = "mk", alternative = "two.sided")TrendTest(x, method = "mk", alternative = "two.sided")
x |
a numeric vector or a data.frame with dates in the first column and chronologically ordered variable in the second. |
method |
a choice of test method. Choices are "mk" (Mann Kendall - the default), "pearson", and "spearman". |
alternative |
the alternative hypothesis. Options are "less", "greater", and "two.sided". The default is "two.sided". |
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.
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.
Anthony Hammond
#Get 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 method with dates included and not TrendTest(AM.27083, method = "pearson") TrendTest(AM.27083$Flow, method = "pearson") #Apply the default mk test for positive trend TrendTest(AM.27083$Flow, alternative = "greater")#Get 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 method with dates included and not TrendTest(AM.27083, method = "pearson") TrendTest(AM.27083$Flow, method = "pearson") #Apply the default mk test for positive trend TrendTest(AM.27083$Flow, alternative = "greater")
UAF and PRUAF from catchment descriptors for QMED estimation in ungauged urban catchments
UAF(CDs = NULL, URBEXT2000, BFIHOST)UAF(CDs = NULL, URBEXT2000, BFIHOST)
CDs |
catchment descriptors derived from either GetCDs or CDsXML |
URBEXT2000 |
quantification of catchment urbanisation (used when CDs is not) |
BFIHOST |
baseflow index as a function of hydrological soil type of the catchment (used when CDs is not) |
a data.frame with columns UAF and PRUAF
Anthony Hammond
#Get some catchment descriptors for an urban catchment calculate the UAF & PRUAF CDs.53006 <- GetCDs(53006) UAF(CDs.53006) #Calculate UAF and PRUAF using a user input URBEXT2000 and BFIHOST UAF(URBEXT2000 = 0.1138, BFIHOST = 0.3620)#Get some catchment descriptors for an urban catchment calculate the UAF & PRUAF CDs.53006 <- GetCDs(53006) UAF(CDs.53006) #Calculate UAF and PRUAF using a user input URBEXT2000 and BFIHOST UAF(URBEXT2000 = 0.1138, BFIHOST = 0.3620)
This function provides a coefficient to multiply by URBEXT2000 to adjust it to a given year
UEF(Year)UEF(Year)
Year |
The year for consideration. Numeric |
The urban expansion factor is detailed in Bayliss, A. Black, K. Fava-Verde, A. Kjeldsen, T. (2006). URBEXT2000 - A new FEH catchment descriptor: Calculation, dissemination and application. R&D Technical Report FD1919/TR, DEFRA, CEH Wallingford
A numeric urban expansion factor.
Anthony Hammond
# Get an expansion factor for the year 2023 UEF(2023)# Get an expansion factor for the year 2023 UEF(2023)
Easting and northing national grid reference points around the coast of the UK
UKOutlineUKOutline
A data frame with 3867 rows and 2 variables
Easting, British national grid reference
Northing, British national grid reference
https://environment.data.gov.uk/
Quantification of uncertainty for pooling results for the gauged and ungauged case
Uncertainty( x, Gauged = FALSE, qmed = NULL, Dist = "GenLog", Conf = 0.95, fseQMED = 1.55, UrbAdj = FALSE, URBEXT = NULL, Plot = TRUE, IncAMest = TRUE, Parametric = TRUE )Uncertainty( x, Gauged = FALSE, qmed = NULL, Dist = "GenLog", Conf = 0.95, fseQMED = 1.55, UrbAdj = FALSE, URBEXT = NULL, Plot = TRUE, IncAMest = TRUE, Parametric = TRUE )
x |
the pooled group derived from the Pool() or PoolSmall function |
Gauged |
a logical argument with a default of FALSE. If FALSE the uncertainty intervals are calculated for the ungauged case. If TRUE they are calculated for the gauged case |
qmed |
the QMED estimate for the ungauged case. It is derived from the observed AMAX if Gauged equals TRUE. |
Dist |
a choice of distribution to use for the estimates. Choices are "GEV", "GenLog" ,"Gumbel", or "Kappa3". The default is "GenLog" |
Conf |
the confidence level of the uncertainty intervals. Default is 0.95. Must be between 0 and 1. |
fseQMED |
The factorial standard error of the QMED estimate for an ungauged assessment. The default is 1.55. |
UrbAdj |
applies an urban adjustment to the growth curves |
URBEXT |
URBEXT value for the site of interest. This is necessary if a UrbAdj equals TRUE. |
Plot |
logical argument with a default of TRUE. If TRUE a return level plot with results and margin of error is plotted. If FALSE, it is not. |
IncAMest |
logical argument with a default of TRUE. Sometimes when doing gauged (enhanced single site analysis), the central estimate of the single site estimate is outside the intervals of the ESS estimate. When this argument is true the confidence interval is expanded to include the central estimate for the single site. If FALSE, it is not. |
Parametric |
logical argument with a default of TRUE. If TRUE, the bootsrapping is done by simulation with the distribution of choice. If FALSE the bootsrapping is done by resampling with replacement. |
Uncertainty for both the gauged (enhanced single site) and ungauged case are quantified according to the bootstrapping procedures, which account for weights in the pooling group, detailed in Hammond, A. (2021). Sampling uncertainty of UK design flood estimation. Hydrology Research. 1357-1371. 52 (6).
A dataframe with 10 rows and four columns. Return period in the first column, central estimate in the second, lower in the third, and upper in the fourth. If Plot = TRUE, a return level plot is also returned.
Anthony Hammond
#Get an ungauged pooling group: Pool.203018 <- Pool(GetCDs(203018), exclude = 203018) #Use the function to quantify the central estimate and uncertainty Uncertainty(Pool.203018, qmed = QMED(GetCDs(203018))) #Form a pooling group with subject site included. Pool.203018 <- Pool(GetCDs(203018)) #Quantify the central estimate and uncertainty Uncertainty(Pool.203018, Gauged = TRUE)#Get an ungauged pooling group: Pool.203018 <- Pool(GetCDs(203018), exclude = 203018) #Use the function to quantify the central estimate and uncertainty Uncertainty(Pool.203018, qmed = QMED(GetCDs(203018))) #Form a pooling group with subject site included. Pool.203018 <- Pool(GetCDs(203018)) #Quantify the central estimate and uncertainty Uncertainty(Pool.203018, Gauged = TRUE)
Provides the gauged Lcv weights for each site in a pooling group
WeightsGLcv(x)WeightsGLcv(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
A data.frame with site references in the first column and associated weights in the second
Anthony Hammond
#Get some CDs, form a gauged pooling group, and estimate gauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WeightsGLcv(Pool.27051)#Get some CDs, form a gauged pooling group, and estimate gauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WeightsGLcv(Pool.27051)
Provides the gauged LSkew weights for each site in a pooling group
WeightsGLSkew(x)WeightsGLSkew(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
A data.frame with site references in the first column and associated weights in the second
Anthony Hammond
#Get some CDs, form a gauged pooling group, and estimate gauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WeightsGLSkew(Pool.27051)#Get some CDs, form a gauged pooling group, and estimate gauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WeightsGLSkew(Pool.27051)
Provides the ungauged Lcv weights for each site in a pooling group
WeightsUnLcv(x)WeightsUnLcv(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
A data.frame with site references in the first column and associated weights in the second
Anthony Hammond
#Get some CDs, form an ungauged pooling group, and estimate ungauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WeightsUnLcv(Pool.27051)#Get some CDs, form an ungauged pooling group, and estimate ungauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WeightsUnLcv(Pool.27051)
Provides the ungauged LSkew weights for each site in a pooling group
WeightsUnLSkew(x)WeightsUnLSkew(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
A data.frame with site references in the first column and associated weights in the second
Anthony Hammond
#Get some CDs, form an ungauged pooling group, and estimate ungauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WeightsUnLSkew(Pool.27051)#Get some CDs, form an ungauged pooling group, and estimate ungauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WeightsUnLSkew(Pool.27051)
Calculates the gauged weighted Lcv from a pooling group (enhanced single site)
WGaugLcv(x)WGaugLcv(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
the gauged weighted Lcv from a pooling group
Anthony Hammond
#Get some CDs, form a gauged pooling group, and estimate gauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WGaugLcv(Pool.27051)#Get some CDs, form a gauged pooling group, and estimate gauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WGaugLcv(Pool.27051)
Calculates the gauged weighted LSkew from a pooling group (enhanced single site)
WGaugLSkew(x)WGaugLSkew(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
the gauged weighted LSkew from a pooling group
Anthony Hammond
#Get some CDs, form a gauged pooling group, and estimate gauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WGaugLSkew(Pool.27051)#Get some CDs, form a gauged pooling group, and estimate gauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051) WGaugLSkew(Pool.27051)
Calculates the ungauged weighted Lcv from a pooling group
WungLcv(x)WungLcv(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
the ungauged weighted Lcv from a pooling group
Anthony Hammond
#Get some CDs, form an ungauged pooling group, and estimate ungauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WungLcv(Pool.27051)#Get some CDs, form an ungauged pooling group, and estimate ungauged Lcv CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WungLcv(Pool.27051)
Calculates the ungauged weighted LSkew from a pooling group
WungLSkew(x)WungLSkew(x)
x |
pooling group derived with the Pool() function |
Weighting method as according to Science Report: SC050050 - Improving the FEH statistical procedures for flood frequency estimation
the ungauged weighted LSkew from a pooling group
Anthony Hammond
#Get some CDs, form an ungauged pooling group, and estimate ungauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WungLSkew(Pool.27051)#Get some CDs, form an ungauged pooling group, and estimate ungauged LSkew CDs.27051 <- GetCDs(27051) Pool.27051 <- Pool(CDs.27051, exclude = 27051) WungLSkew(Pool.27051)
Calculates the goodness of fit score for pooling groups.
Zdists(x)Zdists(x)
x |
pooling group derived from the Pool() function |
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).
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.
Anthony Hammond
#Get CDs, form a pooling group and calculate the Zdist CDs.203018 <- GetCDs(203018) Pool.203018 <- Pool(CDs.203018) Zdists(Pool.203018)#Get CDs, form a pooling group and calculate the Zdist CDs.203018 <- GetCDs(203018) Pool.203018 <- Pool(CDs.203018) Zdists(Pool.203018)