| Title: | Analyte Flux and Load from Estimates of Concentration and Discharge |
|---|---|
| Description: | Flux (mass per unit time) and Load (mass) are computed from timeseries estimates of analyte concentration and discharge. Concentration timeseries are computed from regression between surrogate and user-provided analyte. Uncertainty in calculations is estimated using bootstrap resampling. Code for the processing of acoustic backscatter from horizontally profiling acoustic Doppler current profilers is provided. All methods detailed in Livsey et al (2020) <doi:10.1007/s12237-020-00734-z>, Livsey et al (2023) <doi:10.1029/2022WR033982>, and references therein. |
| Authors: | Daniel Livsey [aut, cre, cph] |
| Maintainer: | Daniel Livsey <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.1 |
| Built: | 2025-02-21 04:31:40 UTC |
| Source: | https://github.com/dlivsey/realtimeloads |
Processes acoustic backscatter from horizontally profiling ADCP (hADCP). Returns attenuation of sound due to water and suspended-sediment. Applies all corrections to acoustic backscatter detailed in the guideline.
acoustic_backscatter_processing( Site, ADCP, Height, Sonde, Echo_Intensity_Beam_1, Echo_Intensity_Beam_2, Instrument_Noise_Level = NULL, Include_Rayleigh = FALSE, Include_near_field_correction = TRUE )acoustic_backscatter_processing( Site, ADCP, Height, Sonde, Echo_Intensity_Beam_1, Echo_Intensity_Beam_2, Instrument_Noise_Level = NULL, Include_Rayleigh = FALSE, Include_near_field_correction = TRUE )
Site |
Data frame with site, local vertical datum, and ADCP elevation information
|
ADCP |
Data frame with various readings from ADCP
|
Height |
Data frame with timeseries of river height
|
Sonde |
Data frame with timeseries of conductivity, temperature, and depth from sonde
|
Echo_Intensity_Beam_1 |
Data frame of acoustic backscatter measurements from beam 2
|
Echo_Intensity_Beam_2 |
Data frame of acoustic backscatter measurements from beam 2
|
Instrument_Noise_Level |
Estimate of noise level, recommended if ambient noise level is not recorded (counts) |
Include_Rayleigh |
Logical to include data within Rayleigh Distance for processing of acoustic backsactter |
Include_near_field_correction |
Logical to include near-field correction of Downing et al (1995) |
List with processed data, all variable names and units are written-out in list items, see Livsey (in review) for details of each variable
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
InputData <- realTimeloads::ExampleData Site <- InputData$Site ADCP <- InputData$ADCP Height <- InputData$Height Sonde <- InputData$Sonde EIa <- InputData$Echo_Intensity # example code assumes backscatter is equal across beams EIb <- InputData$Echo_Intensity Output <- acoustic_backscatter_processing(Site,ADCP,Height,Sonde,EIa,EIb)InputData <- realTimeloads::ExampleData Site <- InputData$Site ADCP <- InputData$ADCP Height <- InputData$Height Sonde <- InputData$Sonde EIa <- InputData$Echo_Intensity # example code assumes backscatter is equal across beams EIb <- InputData$Echo_Intensity Output <- acoustic_backscatter_processing(Site,ADCP,Height,Sonde,EIa,EIb)
Computes attenuation of sound in water per Ainslie and McColm (1998)
attenuation_of_sound_by_water(freq, temp, sal)attenuation_of_sound_by_water(freq, temp, sal)
freq |
frequency of sound (Hz) |
temp |
Water temperature (degrees C) |
sal |
Salinity (PSU) |
attenuation of sound in water (dB/m), divide by 20*log10(exp(1)) to convert to Nepers/m
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Ainslie, M. A., & McColm, J. G. (1998). A simplified formula for viscous and chemical absorption in sea water. The Journal of the Acoustical Society of America, 103(3), 1671-1672.
Author modified Matlab code from David Schoellhamer
InputData <- realTimeloads::ExampleData freq <- InputData$ADCP$Accoustic_Frequency_kHz*1000 cond <-InputData$Sonde$Conductivity_uS_per_cm temp <- InputData$Sonde$Water_Temperature_degC dbar <- InputData$Sonde$Pressure_dbar sal <- ctd2sal(cond,temp,dbar) aw <- attenuation_of_sound_by_water(freq,temp,sal) # dB/m awNp <- attenuation_of_sound_by_water(freq,temp,sal)/(20*log10(exp(1))) # Np/mInputData <- realTimeloads::ExampleData freq <- InputData$ADCP$Accoustic_Frequency_kHz*1000 cond <-InputData$Sonde$Conductivity_uS_per_cm temp <- InputData$Sonde$Water_Temperature_degC dbar <- InputData$Sonde$Pressure_dbar sal <- ctd2sal(cond,temp,dbar) aw <- attenuation_of_sound_by_water(freq,temp,sal) # dB/m awNp <- attenuation_of_sound_by_water(freq,temp,sal)/(20*log10(exp(1))) # Np/m
Computes uncertainty in regression parameters of y(x) after Rustomji and Wilkinson (2008)
bootstrap_regression(Calibration, fit_eq, fit_glm = FALSE)bootstrap_regression(Calibration, fit_eq, fit_glm = FALSE)
Calibration |
data frame with surrogate(s) followed by analyte in last column |
fit_eq |
equation used to fit y(x), string (e.g, "y ~ x + x2", "y ~ x", "log10(y)~x') |
fit_glm |
logical to use Generalized Linear Models for models with factor (i.e., categorical) predictors |
list with bootstrap regression parameters and list output from stats::lm()
User should inspect regression residuals and relevant statistics to ensure model form is reasonable, suggested reading: regression diagnostics in Statistical Methods in Water Resources (https://doi.org/10.3133/tm4a3).
One can call plot(fit) to view various regression diagnostic plots
Bias Correction Factor (BCF) is only relevant when analyte is transformed to log units, see https://doi.org/10.3133/tm4a3 to convert a model that used log(analyte) back to linear units use: analyte = 10^(f(surrogates)) x BCF
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Rustomji, P., & Wilkinson, S. N. (2008). Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curves. Water resources research, 44(9).https://doi.org/10.1029/2007WR006088
Helsel, D.R., Hirsch, R.M., Ryberg, K.R., Archfield, S.A., and Gilroy, E.J., 2020, #' Statistical methods in water resources: U.S. Geological Survey Techniques and Methods, book 4, chap. A3, 458 p. https://doi.org/10.3133/tm4a3
# linear model x <- 1:10 y <- 0.5*x + 10 boot <- bootstrap_regression(data.frame(x,y),"y~x") # polynomial model, call to I() needed for squaring x in equation string x <- 1:10 y <- x + x^2 boot <- bootstrap_regression(data.frame(x,y),"y ~ x+I(x^2)") # power law model # BCF returned since y is transformed to log units x <- 1:10 y <- x^0.3 boot <- bootstrap_regression(data.frame(x,y),"log10(y)~log10(x)") # multivariate model a <- 1:10 b <- a*2 c <- a^2*b^3 boot <- bootstrap_regression(data.frame(a,b,c),"log10(c)~log10(a)+log10(b)")# linear model x <- 1:10 y <- 0.5*x + 10 boot <- bootstrap_regression(data.frame(x,y),"y~x") # polynomial model, call to I() needed for squaring x in equation string x <- 1:10 y <- x + x^2 boot <- bootstrap_regression(data.frame(x,y),"y ~ x+I(x^2)") # power law model # BCF returned since y is transformed to log units x <- 1:10 y <- x^0.3 boot <- bootstrap_regression(data.frame(x,y),"log10(y)~log10(x)") # multivariate model a <- 1:10 b <- a*2 c <- a^2*b^3 boot <- bootstrap_regression(data.frame(a,b,c),"log10(c)~log10(a)+log10(b)")
Applies a Butterworth filter with a 30-hour stop period and a 40-hour pass period
butterworth_tidal_filter(time, x)butterworth_tidal_filter(time, x)
time |
time for x (time, POSIXct) |
x |
any quantity, for example discharge (double) |
non-tidal signal in x with data affected by filter ringing removed
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Ruhl, C. A., & Simpson, M. R. (2005). Computation of discharge using the index-velocity method in tidally affected areas (Vol. 2005). Denver: US Department of the Interior, US Geological Survey. https://pubs.usgs.gov/sir/2005/5004/sir20055004.pdf
time <- realTimeloads::ExampleData$Height$time x <- realTimeloads::ExampleData$Height$Height_m xf <- butterworth_tidal_filter(time,x)time <- realTimeloads::ExampleData$Height$time x <- realTimeloads::ExampleData$Height$Height_m xf <- butterworth_tidal_filter(time,x)
Compute load with uncertainty on concentration estimates from bootstrap regression after Rustomji and Wilkinson (2008)
compute_load(Surrogate, Discharge, Regression, period = NULL)compute_load(Surrogate, Discharge, Regression, period = NULL)
Surrogate |
data frame with time (PosixCt) and surrogate(s) (x,...) |
Discharge |
data frame with time (PosixCt) and discharge in cubic meters per second |
Regression |
data frame from bootstrap_regression() that determines analyte(surrogate) |
period |
two element vector time (PosixCt) indicating period over which load is computed |
list with data frames of estimated concentration and flux used to compute load (i.e., the sum of flux)
Surrogate and Discharge time series can be on different time steps
If period is NULL, computes load over time in Surrogate
Discharge should be in cubic meters per second
Analyte concentration estimated from surrogate should be in milligrams per second
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Rustomji, P., & Wilkinson, S. N. (2008). Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curves. Water resources research, 44(9).https://doi.org/10.1029/2007WR006088
Helsel, D.R., Hirsch, R.M., Ryberg, K.R., Archfield, S.A., and Gilroy, E.J., 2020, #' Statistical methods in water resources: U.S. Geological Survey Techniques and Methods, book 4, chap. A3, 458 p. https://doi.org/10.3133/tm4a3
Turbidity_FNU <- realTimeloads::ExampleData$Sonde$Turbidity TSS_mg_per_l <- realTimeloads::ExampleData$Sediment_Samples$SSCpt_mg_per_liter Discharge <- realTimeloads::ExampleData$Discharge Calibration <- data.frame(Turbidity_FNU,TSS_mg_per_l) time <- realTimeloads::ExampleData$Sonde$time Surrogate <- data.frame(time,Turbidity_FNU) Regression = bootstrap_regression(Calibration,'TSS_mg_per_l~Turbidity_FNU') period <- c(as.POSIXct("2000-02-16 AEST"),as.POSIXct("2000-03-16 AEST")) Output <- compute_load(Surrogate,Discharge,Regression,period)Turbidity_FNU <- realTimeloads::ExampleData$Sonde$Turbidity TSS_mg_per_l <- realTimeloads::ExampleData$Sediment_Samples$SSCpt_mg_per_liter Discharge <- realTimeloads::ExampleData$Discharge Calibration <- data.frame(Turbidity_FNU,TSS_mg_per_l) time <- realTimeloads::ExampleData$Sonde$time Surrogate <- data.frame(time,Turbidity_FNU) Regression = bootstrap_regression(Calibration,'TSS_mg_per_l~Turbidity_FNU') period <- c(as.POSIXct("2000-02-16 AEST"),as.POSIXct("2000-03-16 AEST")) Output <- compute_load(Surrogate,Discharge,Regression,period)
Computes salinity from conductivity, water temperature, and depth.
ctd2sal(cond, temp, dbar)ctd2sal(cond, temp, dbar)
cond |
Conductance (microS/cm) |
temp |
Water temperature (degrees C) |
dbar |
Pressure (dBar) or water depth (m) |
Salinity in PSU
If specific conductivity is returned from the sonde, the temperature at which specific conductivity is computed should be utilized
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Fofonoff, N. P., & Millard Jr, R. C. (1983). Algorithms for the computation of fundamental properties of seawater.
Chen, C. T. A., & Millero, F. J. (1986). Thermodynamic properties for natural waters covering only the limnological range 1. Limnology and Oceanography, 31(3), 657-662.
Hill, K., Dauphinee, T., & Woods, D. (1986). The extension of the Practical Salinity Scale 1978 to low salinities. IEEE Journal of Oceanic Engineering, 11(1), 109-112.
Author modified Matlab code from David Schoellhamer
Sonde <- realTimeloads::ExampleData$Sonde sal <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar)Sonde <- realTimeloads::ExampleData$Sonde sal <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar)
Compute uncertainty on timeseries from bootstrap regression after Rustomji and Wilkinson (2008)
estimate_timeseries(Surrogate, Regression)estimate_timeseries(Surrogate, Regression)
Surrogate |
data frame with time (PosixCt) and surrogate(s) (x,...) |
Regression |
data frame from bootstrap_regression() that determines analyte(surrogate) |
list with inputs and uncertainty on timeseries estimated from Regression
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Rustomji, P., & Wilkinson, S. N. (2008). Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curves. Water resources research, 44(9).https://doi.org/10.1029/2007WR006088
Helsel, D.R., Hirsch, R.M., Ryberg, K.R., Archfield, S.A., and Gilroy, E.J., 2020, #' Statistical methods in water resources: U.S. Geological Survey Techniques and Methods, book 4, chap. A3, 458 p. https://doi.org/10.3133/tm4a3
Turbidity_FNU <- realTimeloads::ExampleData$Sonde$Turbidity TSS_mg_per_l <- realTimeloads::ExampleData$Sediment_Samples$SSCpt_mg_per_liter Calibration <- data.frame(Turbidity_FNU,TSS_mg_per_l) time <- realTimeloads::ExampleData$Sonde$time Surrogate <- data.frame(time,Turbidity_FNU) Regression = bootstrap_regression(Calibration,'TSS_mg_per_l~Turbidity_FNU') Output <- estimate_timeseries(Surrogate,Regression)Turbidity_FNU <- realTimeloads::ExampleData$Sonde$Turbidity TSS_mg_per_l <- realTimeloads::ExampleData$Sediment_Samples$SSCpt_mg_per_liter Calibration <- data.frame(Turbidity_FNU,TSS_mg_per_l) time <- realTimeloads::ExampleData$Sonde$time Surrogate <- data.frame(time,Turbidity_FNU) Regression = bootstrap_regression(Calibration,'TSS_mg_per_l~Turbidity_FNU') Output <- estimate_timeseries(Surrogate,Regression)
Computes sediment load per guideline from ExampleData
ExampleCode()ExampleCode()
list with data frames of estimated concentration and flux along with data used in regression and surrogate timeseries
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
realTimeloads Package help file
Output <- ExampleCode()Output <- ExampleCode()
Computes sediment load per guideline from optical and acoustic backscatter measurements combined to the "Sediment Composition Index" (SCI) per Livsey et al (2023)
ExampleCodeSCI()ExampleCodeSCI()
total load with uncertainty computed from estimates of concentration from SCI
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
realTimeloads Package help file
Output <- ExampleCodeSCI()Output <- ExampleCodeSCI()
Synthetic dataset from modeled sediment transport and acoustic scattering detailed in the Appendices of Livsey (in review) Following dataframes are provided in list
ExampleDataExampleData
Site , Site, site datum, and ADCP elevation informationSite name (string)
Unique site code (string)
Elevation of the ADCP above the bed (m)
Elevation of the ADCP above local gauge datum (m)
Distance from local gauge datum to lower point in cross-section (m)
Installation date of ADCP (time, POSIXct)
Date if/when ADCP is moved vertically (time, POSIXct)
User comment (string)
ADCP , ADCP readings except acoustic backscatterUnique site code (string)
Date and time (time, POSIXct)
Measurment ensemble number (integer)
Acoustic frequency of ADCP (kHz)
Radius of ADCP transducer (m)
Angle of beam relative to normal (degrees)
Ratio of beam radius to beam length (-)
Number of measurement cells along beam (integer)
Cell width measured normal to ADCP (m)
Blanking distance measured normal to ADCP (m)
Serial number of ADCP instrument (string)
Serial number of ADCP CPU (string)
Ambient noise level for beam 1 (counts)
Ambient noise level for beam 2 (counts)
Reported distance normal to ADCP to midpoint of bin/cell (m)
Speed of sound used by ADCP in the field (m/s)
Temperature recorded by ADCP (degrees C)
Pressure recorded by ADCP (dBar)
Distance to water surface reported by vertical beam of ADCP (m)
Power to ADCP (V)
Echo_Intensity , Acoustic backscatter measurements from ADCPUnique site code (string)
Date and time (time, POSIXct)
Acoustic backscatter in nth cell (counts)
Sonde , Conductivity, temperature, and depth from sondeDate and time (time, POSIXct)
Temperature (degrees C)
Conductivity (microS/cm)
Pressure (dbar)
Turbidity (FNU)
Unique site code (string)
Height , River height in meters referenced to gauge datumDate and time (time, POSIXct)
Water surface elevation above gauge datum (m)
Unique site code (string)
Discharge, Discharge timeseries in cubic meters per secondDate and time (time, POSIXct)
Dischage (cubic meters per second)
Unique site code (string)
Sediment_Samples , Measured sediment concentration in milligrams per liter (SSC or TSS)Date and time (time, POSIXct)
Concentration of suspended-sediment in milligrams per liter, depth-averaged and velocity weighted average for cross-section
Concentration of suspended-sediment in milligrams per liter, measured at-a-point at elevation of hADCP
Unique site code (string)
data(ExampleData) # lazy-load ony, unable to inspect contents in Rstudio
names(ExampleData) # load data for inspection in Rstudio and view names of items in list
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
Computes sediment load per guideline from user data in list "InputData" generated by function import_data()
hADCPLoads(InputData)hADCPLoads(InputData)
InputData |
List generated by import_data.R |
list with data frames of estimated concentration and flux along with data used in regression and surrogate timeseries
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
import_data Import data from files in user-specified folder
# loads example data in package folder extdata InputData <- import_data() # import_data(path) can be used to import user data Output <- hADCPLoads(InputData)# loads example data in package folder extdata InputData <- import_data() # import_data(path) can be used to import user data Output <- hADCPLoads(InputData)
Imports csv files to R, file names, variable names (and units) in csv text files must match variable names used in ExampleData.rda
import_data(data_folder)import_data(data_folder)
data_folder |
file path to folder containing .txt csv files with format that matches files in extdata package folder |
list with data frames used in package code, see ?ExampleData for list format
Synthetic data used in ExampleData only has backscatter for one beam ("ADCP_Echo_Intensity.txt"), for user data, one should have backscatter for two beams with following names: "ADCP_Echo_Intensity_Beam_1.txt" and "ADCP_Echo_Intensity_Beam_2.txt"
Package arguments require variable names and units to match the names and variable units provided (see ?ExampleData, or .txt files in extdata folder)
Suggest saving all csv files in .txt format to ensure time format is not changed when editing/saving csv in Excel
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Livsey, D.N. (in review). National Industry Guidelines for hydrometric monitoring–Part 12: Application of acoustic Doppler velocity meters to measure suspended-sediment load. Bureau of Meteorology. Melbourne, Australia.
hADCPLoads Process acoustic backscatter from hADCP and compute load using InputData from import_Data()
InputData <- import_data() # loads text files provided in package folder "extdata"InputData <- import_data() # loads text files provided in package folder "extdata"
Returns x with gaps imputed using ARIMA and Decision Trees with option to use harmonic model as predictors for x in decision tree algorithm. Uncertainty on imputed data is estimated using using Monte Carlo (MC) resampling adapting methods of Rustomji and Wilkinson (2008)
impute_data( time, x, Xreg = NULL, ti = NULL, hfit = NULL, harmonic = FALSE, only_use_Xreg = FALSE, MC = 1, ptrain = 1 )impute_data( time, x, Xreg = NULL, ti = NULL, hfit = NULL, harmonic = FALSE, only_use_Xreg = FALSE, MC = 1, ptrain = 1 )
time |
time for x (time, POSIXct) |
x |
any quantity (double) |
Xreg |
additional predictors for decision tree, required if harmonic is FALSE (rows = time, or if given, ti) |
ti |
time vector for interpolation (time, POSIXct) |
hfit |
model object from TideHarmonics::ftide |
harmonic |
logical if x exhibits tidal or diurnal variability |
only_use_Xreg |
logical for using Xreg only in decision tree |
MC |
number of Monte Carlo simulations for uncertainty estimation |
ptrain |
proportion of data used for training and testing model |
list with x imputed at time or ti, if given. Uncertainty estimated from Monte Carlo simulations
If MC == 1, uncertainty is not evaluated. If ptrain == 1, uncertainty and validation accuracy are not computed
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Rustomji, P., & Wilkinson, S. N. (2008). Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curves. Water resources research, 44(9).
van Buuren S, Groothuis-Oudshoorn K (2011). “mice: Multivariate Imputation by Chained Equations in R.” Journal of Statistical Software, 45(3), 1-67. doi:10.18637/jss.v045.i03.
Stephenson AG (2016). Harmonic Analysis of Tides Using TideHarmonics. https://CRAN.R-project.org/package=TideHarmonics.
Moritz S, Bartz-Beielstein T (2017). “imputeTS: Time Series Missing Value Imputation in R.” The R Journal, 9(1), 207–218. doi:10.32614/RJ-2017-009.
# Impute non-tidal data time <- realTimeloads::ExampleData$Sediment_Samples$time xo <- realTimeloads::ExampleData$Sediment_Samples$SSCxs_mg_per_liter Q <- realTimeloads::ExampleData$Discharge$Discharge_m_cubed_per_s idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples flow_concentrtion_ratio <- imputeTS::na_interpolation(Q/x) Xreg <- cbind(Q,flow_concentrtion_ratio) Output <- impute_data(time,x,Xreg,MC = 10,ptrain = 0.8) # Impute tidal data time <-TideHarmonics::Portland$DateTime[1:(24*90)] xo <-TideHarmonics::Portland$SeaLevel[1:(24*90)] idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples Output <- impute_data(time,x,harmonic = TRUE,MC = 10,ptrain = 0.8)# Impute non-tidal data time <- realTimeloads::ExampleData$Sediment_Samples$time xo <- realTimeloads::ExampleData$Sediment_Samples$SSCxs_mg_per_liter Q <- realTimeloads::ExampleData$Discharge$Discharge_m_cubed_per_s idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples flow_concentrtion_ratio <- imputeTS::na_interpolation(Q/x) Xreg <- cbind(Q,flow_concentrtion_ratio) Output <- impute_data(time,x,Xreg,MC = 10,ptrain = 0.8) # Impute tidal data time <-TideHarmonics::Portland$DateTime[1:(24*90)] xo <-TideHarmonics::Portland$SeaLevel[1:(24*90)] idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples Output <- impute_data(time,x,harmonic = TRUE,MC = 10,ptrain = 0.8)
Returns x with gaps imputed using forecast::auto.arima with optional uncertainty estimation using Monte Carlo resampling. Uncertainty on imputed data is estimated using using Monte Carlo (MC) resampling adapting methods of Rustomji and Wilkinson (2008). forecast::auto.arima searches for best ARIMA model and allows for step-changes in the ARIMA model of x(Xreg)
impute_data_with_arima(time, x, Xreg = NULL, ti = NULL, MC = 1, ptrain = 0.8)impute_data_with_arima(time, x, Xreg = NULL, ti = NULL, MC = 1, ptrain = 0.8)
time |
time for x (time, POSIXct) |
x |
any quantity (double) |
Xreg |
additional predictors for ARIMA |
ti |
time vector for interpolation (time, POSIXct) |
MC |
number of Monte Carlo simulations for uncertainty estimation |
ptrain |
proportion of data used for training and testing model |
list with x imputed at time or ti, if given. Uncertainty estimated from Monte Carlo simulations
For infilling missing suspended-sediment concentration data (e.g., TSS or SSC) I suggest using the Froude number and discharge for Xreg. For tidally affected sites, consider using tidally filtered discharge.
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Rustomji, P., & Wilkinson, S. N. (2008). Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curves. Water resources research, 44(9).
Hyndman R, Athanasopoulos G, Bergmeir C, Caceres G, Chhay L, O'Hara-Wild M, Petropoulos F, Razbash S, Wang E, Yasmeen F (2023). forecast: Forecasting functions for time series and linear models. R package version 8.21.1, https://pkg.robjhyndman.com/forecast/.
Hyndman RJ, Khandakar Y (2008). “Automatic time series forecasting: the forecast package for R.” Journal of Statistical Software, 26(3), 1–22. doi:10.18637/jss.v027.i03.
# Impute non-tidal data time <- realTimeloads::ExampleData$Sediment_Samples$time xo <- realTimeloads::ExampleData$Sediment_Samples$SSCxs_mg_per_liter Q <- realTimeloads::ExampleData$Discharge$Discharge_m_cubed_per_s idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples Output <- impute_data_with_arima(time,x,Xreg = Q)# Impute non-tidal data time <- realTimeloads::ExampleData$Sediment_Samples$time xo <- realTimeloads::ExampleData$Sediment_Samples$SSCxs_mg_per_liter Q <- realTimeloads::ExampleData$Discharge$Discharge_m_cubed_per_s idata <- sample(1:length(xo),round(length(xo)*0.5),replace=FALSE) x <- rep(NA,length(xo)) x[idata] <- xo[idata] # simulated samples Output <- impute_data_with_arima(time,x,Xreg = Q)
Linear interpolation limited by time since previous or following reading
linear_interpolation_with_time_limit(time, x, ti, threshold)linear_interpolation_with_time_limit(time, x, ti, threshold)
time |
time for x (time, POSIXct) |
x |
any quantity, for example discharge (double) |
ti |
time where time(x) will be interpolated to (time, POSIXct) |
threshold |
maximum duration where interpolation is allowed (hours) |
a data frame with time (ti), x interpolated from time(x) onto ti, and logical (ibad) if interpolation exceeded threshold
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Dowle M, and others (2023). data.table: Extension of 'data.frame'. https://cran.r-project.org/web/packages/data.table
InputData <- realTimeloads::ExampleData ADCP <- InputData$ADCP Height <- InputData$Height # Interpolate river height to ADCP time time <- realTimeloads::ExampleData$Height$time x <- realTimeloads::ExampleData$Height$Height_m ti <-realTimeloads::ExampleData$ADCP$time threshold <- 1 Output<- linear_interpolation_with_time_limit(time,x,ti,threshold)InputData <- realTimeloads::ExampleData ADCP <- InputData$ADCP Height <- InputData$Height # Interpolate river height to ADCP time time <- realTimeloads::ExampleData$Height$time x <- realTimeloads::ExampleData$Height$Height_m ti <-realTimeloads::ExampleData$ADCP$time threshold <- 1 Output<- linear_interpolation_with_time_limit(time,x,ti,threshold)
Computes dimensionless near-field correction
near_field_correction(freq, c, r, at)near_field_correction(freq, c, r, at)
freq |
Frequency of sound (Hz) |
c |
Speed of sound in water (m/s) |
r |
range to cell center measured along-beam (m) |
at |
Radius of ADCP transducer (m) |
Near-field correction (dimensionless)
See various references cautioning use of near-field correction (e.g., https://doi.org/10.1002/2016WR019695)
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Downing, A., Thorne, P. D., & Vincent, C. E. (1995). Backscattering from a suspension in the near field of a piston transducer. The Journal of the Acoustical Society of America, 97(3), 1614-1620.
InputData <- realTimeloads::ExampleData Sonde<- InputData$Sonde freq <- InputData$ADCP$Accoustic_Frequency_kHz[1]*1000 S <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) c <- speed_of_sound(S,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) at <- InputData$ADCP$Transducer_radius_m r <- seq(0.1,10,0.1) psi <- near_field_correction(freq,c[1],r,at[1])InputData <- realTimeloads::ExampleData Sonde<- InputData$Sonde freq <- InputData$ADCP$Accoustic_Frequency_kHz[1]*1000 S <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) c <- speed_of_sound(S,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) at <- InputData$ADCP$Transducer_radius_m r <- seq(0.1,10,0.1) psi <- near_field_correction(freq,c[1],r,at[1])
Computes speed of sound in water per Del grosso (1974)
speed_of_sound(sal, temp, depth)speed_of_sound(sal, temp, depth)
sal |
Salinity (PSU) |
temp |
Water temperature (degrees C) |
depth |
Water depth (m) or pressure (dBar) |
Speed of sound in water (m/s)
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
Del Grosso, V. A. (1974). New equation for the speed of sound in natural waters (with comparisons to other equations). The Journal of the Acoustical Society of America, 56(4), 1084-1091. Author modified matlab code from David Schoellhamer
InputData <- realTimeloads::ExampleData Sonde<- InputData$Sonde sal <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) c <- speed_of_sound(sal,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar)InputData <- realTimeloads::ExampleData Sonde<- InputData$Sonde sal <- ctd2sal(Sonde$Conductivity_uS_per_cm,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar) c <- speed_of_sound(sal,Sonde$Water_Temperature_degC,Sonde$Pressure_dbar)
Interpolate timeseries x(tx) onto y(ty) with temporal threshold on interpolation
surrogate_to_analyte_interpolation(tx, x, ty, y, threshold)surrogate_to_analyte_interpolation(tx, x, ty, y, threshold)
tx |
time for x "surrogate" (time, POSIXct) |
x |
quantity used to estimate y, for example, accoustic backscatter |
ty |
time for y "analyte" (time, POSIXct) |
y |
measured quantity, for example, an analyte such as suspended-sediment concentration |
threshold |
maximum duration where interpolation is allowed (minutes) |
a data frame with surrogate (x) interpolated onto timestep of analyte (y), interpolated values exceeding threshold are excluded from the output
Daniel Livsey (2023) ORCID: 0000-0002-2028-6128
tx <- as.POSIXct(seq(0,24*60^2,60*1), origin = "2000-01-01",tz = "Australia/Brisbane") x <- sin(1:length(tx)) ty <- as.POSIXct(seq(0,24*60^2,60*15), origin = "2000-01-01",tz = "Australia/Brisbane") y <- seq(0,24*60^2,60*15) threshold <- 10 calibration <- surrogate_to_analyte_interpolation(tx,x,ty,y,threshold)tx <- as.POSIXct(seq(0,24*60^2,60*1), origin = "2000-01-01",tz = "Australia/Brisbane") x <- sin(1:length(tx)) ty <- as.POSIXct(seq(0,24*60^2,60*15), origin = "2000-01-01",tz = "Australia/Brisbane") y <- seq(0,24*60^2,60*15) threshold <- 10 calibration <- surrogate_to_analyte_interpolation(tx,x,ty,y,threshold)