Nothing
R/GPPFourier.R
In GPPFourier: Calculate Gross Primary Production (GPP) from O2 Time Series
Defines functions
GPPFourierPreprocess
GPPFourier
WindowGPPFourier.gts
WindowGPPFourier
demod
GPPFourier_t
Documented in
demod
GPPFourier
GPPFourierPreprocess
GPPFourier_t
WindowGPPFourier
WindowGPPFourier.gts
#
# Core functions to calculate GPP from O2 time series making use of Fourier method
#
#'
#' @title Pre-process O2 data for GPP calculation
#' @description Removes longer term trend of time series and confine series to integer
#' number of days
#' @param x A vector containing consecutive O2 concentrations, sampled at regular intervals (time step = dt).
#' @param dt Sampling time step. Class {difftime}, or \code{numeric} when units are provided.
#' @param units Units of sampling time step
#' @param Detrend Toggle time series detrending
#' @param filter Toggle time series filtering
#' @param Nfilt Moving average filter width (see documentation of
#' \code{filter()})
#' @param circular Moving average boundary condition (see documentation of
#' \code{filter()})
#' @param sides Moving average central or one sided (see documentation of
#' \code{filter()})
#' @details GPPFourierPreprocess() is called by GPPFourier() to pre-process the time series before calculating the Fourier amplitude at diel frequency and estimate GPP.
#' @return A list with components \item{filt}{The long term trend}
#' \item{res}{The original series without the long term trend, confined to an
#' integer number of days} \item{indices}{The indices of the subset of the
#' series that are retained, after filtering and confining to an integer number
#' of days.}
#' @seealso \code{\link{GPPFourier}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier}},
#' \code{\link{WindowGPPFourier.gts}}
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#' @keywords utilities
#' @export GPPFourierPreprocess
GPPFourierPreprocess <- function(x, # O2 time series from which GPP is to be calculated
dt=NULL, # Sampling time step
units=c("days", "hours", "mins"), # Time units of sampling time step
Detrend=TRUE, # Toggle time series detrending
filter = TRUE, # Toggle time series filtering
Nfilt=NULL, # Moving average filter width (see documentation of filter())
circular = FALSE, # Moving average boundary condition (see documentation of filter())
sides = 2) # Moving average central or one sided (see documentation of filter())
{
# Removes longer term trends and confines time series to integer number of periods
#
#
unit <- match.arg(units)
if (is.null(dt)) stop("Provide dt")
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}}
if (filter&is.null(Nfilt)) stop("Provide filter length (Nfilt)")
if (Nfilt>length(x)) stop("Filter longer than time series")
dt <- as.double(dt, units="days")
# confine series to integer number of periods (T=1 day = 1/dt samples)
Nperiod <- floor(length(x)*dt)/dt
indices <- 1:Nperiod
x <- x[indices]
if (Detrend){
detrendx <- detrend(x, returntrend=TRUE)
x <- detrendx$detrended
trendx <- detrendx$trend
} else {
trendx <- rep(mean(x), length(x))
}
if (filter){
cfilt <- rep(1/Nfilt,Nfilt)
filtx <- filter(x,cfilt,sides=sides,circular=circular)
resx <- x - filtx
indices <- indices[!is.na(filtx)]
} else {
filtx <- rep(mean(x),length(x))
resx <- x - mean(x)
}
# confine filtered series to integer number of periods (T=1 day = 1/dt samples)
filtx <- filtx[!is.na(filtx)]
resx <- resx[!is.na(resx)]
Nperiod <- floor(length(filtx)*dt)/dt
indices <- indices[1:Nperiod]
filtx <- filtx[1:Nperiod]
resx <- resx[1:Nperiod]
return(list(filt = filtx, res=resx, trend=trendx, indices=indices))
}
#'
#' @title Calculate GPP from O2 time series
#'
#' @description This is the workhorse to calculate aquatic Gross Primary Production from
#' high frequency dissolved oxygen data. It will often be more convenient to use the higher
#' level functions \code{\link{WindowGPPFourier}} and
#' \code{\link{WindowGPPFourier.gts}}
#'
#' @param x Regularly sampled O2 time series. x can be a dataframe with time (POSIXt) in first column and O2 concentrations in
#' second column or a vector of concentrations sampled with time step = dt
#' @param dt Sampling time step. Not used when x is a data frame, then dt is calculated #' from the time spacing between
#' the first two samples. Either a difftime object or a numerical value. When dt is given as a numerical value, \code{units}
#' have to be provided.
#' @param units Unit of sampling time step
#' @param Nfilt Moving average filter width
#' @param fDL Relative fraction of light hours during the day. When \code{x} is a data frame, and no value for fDL is provided,
#' fDL is calculated with \code{\link{fDLfun}}
#' @param filtcorrect Logical controlling whether GPP estimate
#' is corrected for signal falsely removed by filtering
#' @param ... Other parameters to be passed to \code{\link{GPPFourierPreprocess}}
#' @param phi Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
#' @param lambda Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
#' @return Volume specific average gross primary production in units of oxgen as in \code{x}.
#' per time unit (unit of \code{dt}, given by \code{units}).
#' Specifically: if O2 concentrations are in mg/L, and \code{units = days}, calculated GPP has units mg/L/day.
#' If surface specific primary production (e.g. per m2) is desired, this number has to be multiplied by the depth of the water body.
#' @details no details
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier}},
#' \code{\link{WindowGPPFourier.gts}}
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#'
#' @keywords utilities
#' @examples
#' DO <- Kruibeke[Kruibeke$time>="2010-06-03"&Kruibeke$time<="2010-06-13",]
#' dt <- as.numeric(difftime(DO$time[2],DO$time[1],units="days"))
#' DL <- fDLfun(DO$time[1], phi=51.176, lambda=4.326)
#' GPP <- GPPFourier(DO$O2,dt=dt,Nfilt=100,fDL=DL)
#'
#' @export GPPFourier
GPPFourier <- function( x, # O2 time series from which GPP is to be calculated
dt=NULL, # Sampling time step
units=c("days","hours","mins"), # Time units of sampling time step
Nfilt=NULL, # Moving average filter width
fDL=NULL, # Relative fraction of light hours during the day
filtcorrect = FALSE, # Logical controlling whether GP estimate is corrected for signal removed by moving average filtering
phi=NULL, # Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
lambda=NULL, # Optional: Longitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
...) # Other parameters to be passed to GPPFourierprocess
{
# Calculate average gross production from O2 time series
# Units: unit of oxgen as in ts per days
#
#
# Error handling
#
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (inherits(t, "POSIXt")){
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
} else {
if (any(0!=round(diff(diff(t)), digits=5))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- diff(c(t[1],t[2]))
}
if (is.null(fDL)){
if (is.null(phi)|is.null(lambda)) stop("Provide longitude and lattitude") else {
fDL <- fDLfun(t, phi=phi, lambda=lambda)
fDL <- mean(fDL) # average fDL of time series is used
}
}
O2 <- ts[,2]
}
} else {
if (is.null(dt)) stop("Provide dt")
if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
if (any(is.na(ts))) stop("NAs not allowed")
O2 <- ts
# fDL <- rep(fDL, length(O2))
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
# if (!inherits(dt,"difftime")){
# if (is.null(units)) {
# stop("Either provide dt as difftime object, or provide units of dt")
# } else {
# dt <- as.difftime(dt, units=unit)
# }}
# if (is.null(Nfilt)) stop("Provide moving average filter width (Nfilt)")
# if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
# if (as.numeric(length(x) < Nfilt)) stop("Filter longer than time series")
# angular frequency (rad per day) associated with diurnal periodicity, i.e. T = 1 day
w <- 2*pi/1
# Remove longer term trends and confines time series to integer number of periods
xpreprocessed <- GPPFourierPreprocess(x, dt=dt, Nfilt=Nfilt, ...)
filtx <- xpreprocessed$filt
resx <- xpreprocessed$res
indices <- xpreprocessed$indices
dt <- as.double(dt,units="days")
# direct calculation of A24
times <- (1:length(resx))*dt
A24 <- 2*sum(exp(-times*1i*w)*resx)/length(resx)
# correction for truncated sinusoidal form of GPP
theta1 <- pi*(1-fDL)
truncFac <- (sin(theta1) + (pi-theta1) * cos(theta1))/(pi-theta1 + 0.5 * sin(2*theta1))
A24 <- A24*truncFac
# Calculate time derivative in frequency domain
A24 <- 1i*w*A24
# Correct for frequency response of moving average filter
if (filtcorrect){
if (is.null(Nfilt)){
warning("Nfilt=NULL. No filter correction calculated")
} else {
L <- Nfilt
maxfreq <- (length(resx) - 1)/(length(resx))/dt
w24h <- 2*pi/maxfreq
A24 <- A24/(1-Hw(w24h,Nfilt))
}}
return(Mod(A24))
}
#'
#' @title Calculate GPP from O2 time series with gaps, in consecutive time-windows of prescribed
#' length
#'
#' @description Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length. High level function to handle time series with gaps. Short gaps are
#' interpolated, series is split at large gaps.
#'
#'
#' @param x Dataframe containing time (POSIXt) in first column and O2 concentrations in the second column.
#' @param dt Sampling time step. Either a difftime object or a numerical value.
#' When \code{dt} is given as a numerical value, \code{units} have to be provided (defaults to \code{days}). When omitted \code{dt} is calculated from the time spacing between the first two samples in the data frame.
#' @param units Unit of sampling time step
#' @param gapMaxN Minimum number of missing data points to be considered a gap in the time series. Gaps smaller than gapMaxN are interpolated. Defaults to the number of samples corresponding to a 4 hour gap.
#' @param Width Width [days] of the time-windows for which GPP is calculated
#' @param filtWidth Length of moving average filter [hours] to filter O2
#' @param ... Other parameters to be passed to WindowGPPFourier()
#' @return Average gross primary production in units of oxgen as in \code{x}
#' per day.
#' @importFrom Tides gapsts
#' @details This high level function analyses multiple regular O2 time series at once, each separated by a gap larger than gapMaxN (the number of . This is convenient e.g. for regular O2 data from a single location, where some data is missing due to sensor maintenance or replacement. \code{x} is assumed to contain all these data, pasted together in a single data frame (e.g. by \code{\link{cbind}}ing all individual regular time series or just reading from a single data file). Additionally, missing values within the regular series are allowed and are interpolated with \code{\link{gapfill}}. Uses \code{\link{WindowGPPFourier}}.
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier_t}}, \code{\link{SunRiseSet}},
#' \code{\link{WindowGPPFourier}}
#' @keywords utilities
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#' @examples
#' par(mfrow=c(2,1))
#' plot(Kruibeke, pch=20,xlim=as.POSIXct(c("2010-01-01", "2010-12-31")))
#' GPPAll_4 <- WindowGPPFourier.gts(Kruibeke,
#' gapMaxN = 10,
#' Width = 10,
#' filtWidth=1*24,
#' phi=51.176,
#' lambda=4.326,
#' Detrend=TRUE,
#' filter=TRUE,
#' filtcorrect=FALSE)
#' plot(GPPAll_4$time,GPPAll_4$GPP*9.3,
#' col="black",
#' pch=19,
#' type="b",
#' xlab="time",
#' ylab="GPP",
#' xlim=as.POSIXct(c("2010-01-01", "2010-12-31")))
#'
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @export WindowGPPFourier.gts
WindowGPPFourier.gts <- function(x,
dt = difftime(x[2,1],x[1,1],units="days"),
units = c("days", "hours", "mins"),
gapMaxN = 4/24/dt,
Width = 14,
filtWidth=16,
...) {
# x: data frame with first column t, and second column O2
# gapMaxN: number of missing samples that will be interpolated. More than gapMaxN missing samples will be handled as a gap
# with dt = 10 min, N = 6 amounts to 1h-gaps being linearly interpolated, 6*4 = 4h gaps interpolated
ts <- x # internal alias
ts$O2interpol <- ts[,2]
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
dt <- as.double(dt, units="days")
# Remove all gaps > dt*gapMaxN
gaps2 <- gapsts(ts[,1], dtMax=dt*gapMaxN, unit="days")
ts$N <- 0
ts$N[is.element(ts[,1],gaps2$t2)] <- 1
ts$N <- cumsum(ts$N)
# interpolate all gaps < dt*gapsMaxN
tslist <- split(ts[,1:2], ts$N)
tsfilledAll <- lapply(tslist, gapfill)
tsfilled <- do.call(rbind, tsfilledAll)
# re-add N
gaps2 <- gapsts(tsfilled[,1], dtMax=dt*gapMaxN, unit="days")
tsfilled$N <- 0
tsfilled$N[is.element(tsfilled$t,gaps2$t2)] <- 1
tsfilled$N <- cumsum(tsfilled$N)
# Windowed GP calculation on 'continuous' subsets of 10d or more, with small gaps filled
tslist <- split(tsfilled[,c("t","y")], tsfilled$N)
dtlist <- lapply(tslist, function(x) {x[1,1]- x[2,1]})
sublengths <- tapply(tsfilled$t, tsfilled$N, length)
tslist <- tslist[sublengths > Width/dt]
GPAll <- lapply(tslist, WindowGPPFourier, Width=Width, filtWidth=filtWidth, ...)
GPAll <- do.call(rbind, GPAll)
return(GPAll)
}
#'
#' @title Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length
#'
#' @description Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length.
#'
#' @param x Vector or dataframe containing consecutive O2 concentrations, sampled at regular intervals (time step = dt). In case of irregular time intervals, consider \code{\link{resample}}.
#' When x is a data frame, time (POSIXt) must be in first column and O2 in second column; when x is a vector, dt must be provided.
#' @param dt Sampling time step. Either a difftime object or a numerical value. When \code{dt} is given as a numerical value,
#' \code{units} have to be provided. Can be omitted when \code{x} is a data frame, then \code{dt} is calculated from the
#' time spacing between the first two samples.
#' @param units Unit of sampling time step
#' @param Width Width [days] of the time-windows for which GPP is calculated
#' @param Nblocks Number of consecutive subsections on which to calculate GPP. By default the number of windows of Width fitting into the range of x[,1]
#' Default the largest number of blocks ow width=Width that
#' @param phi Latitude of location where O2 series was recorded. Used to
#' calculate relative fraction of light hours during the day, during which
#' production takes place
#' @param lambda Longitude of location where O2 series was recorded
#' @param filtWidth [hours] Lenght of moving average filter to filter O2
#' series. See GPPFourierPreprocess()
#' @param fDL Optional. Relative fraction of daylight. When \code{x} is a data frame, fDL is calculated with \code{\link{fDLfun}}
#' @param ... Other parameters to be passed to \code{GPPFourier()}
#' @return Average gross primary production in units of oxgen as in \code{x}
#' per day.
#' @details no details
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier.gts}}, \code{\link{SunRiseSet}}
#' @keywords utilities
#'
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#'
#'
#' @export WindowGPPFourier
WindowGPPFourier <- function(x,
dt=NULL,
units=c("days","hours","mins"),
Width=14,
Nblocks=floor(unclass(difftime(range(ts[,1])[2], range(ts[,1])[1],units="days"))[1]/Width),
fDL=NULL,
phi,
lambda,
filtWidth=16,
...){
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
O2 <- ts[,2]
if (is.null(fDL)){
fDL <- fDLfun(t, phi=phi, lambda=lambda)
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
}
} else {
if (is.null(dt)) stop("Provide dt")
if (any(is.na(ts))) stop("NAs not allowed")
if (is.null(fDL)) {
stop("Provide daylight fraction of day (fDL)")
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
O2 <- ts
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
dtnumeric <- as.double(dt, units="days")
Nfilt <- filtWidth/dtnumeric/24
blockst1 <- t[1] + (0:(Nblocks-1))*86400*Width
blockmids <- t[1] + Width/2*86400 + (0:(Nblocks-1))*86400*Width
blocknum <- approx(x=blockst1,y=0:(Nblocks-1),xout=t,method="constant", rule=2)$y
blocklength <- tapply(t, blocknum, length)
# calculate GPP on all consecutive subsets of Width
GPAll <- NULL
for (i in unique(blocknum))
{
O2i <- O2[blocknum==i]
fDLi <- mean(fDL[blocknum==i])
ti <- t[blocknum==i]
if (length(O2i) < 1.5*Nfilt){
GPAll <- c(GPAll,NA)
} else {
# RiseSet <- SunRiseSet(ti,phi=phi,lambda=lambda) + 2 #summertime = GMT+2
# DL <- diff(RiseSet)
# fDL <- DL/24
GPO2 <- GPPFourier(x=O2i,dt=dt,Nfilt= Nfilt, fDL=fDLi, ...) # Unit: mmol O2/m3/d
GPAll <- c(GPAll,GPO2)
}
}
return(data.frame(time=blockmids[unique(blocknum+1)],GPP=GPAll, t1=blockst1[unique(blocknum+1)], t2=blockst1[unique(blocknum+1)] + Width*86400))
}
#'
#' @title Complex demodulation
#' @description Retreive the time varying amplitude of a signal with fixed carrier frequency by complex demodulation
#' @param x Vector containing time series observations at constant sampling interval
#' @param dt Time step of time series. Arbitrary units. Class \code{Numeric} or \code{difftime}
#' @param w1 Carrier frequency. Inverse unit of \code{dt}.
#' @param filttype Filter type to be applied for low pass filtering. "MA" for moving average; "bw" for butterworth filter.
#' @param Nf Moving average filter width
#' @param nf Number of iterative applications of moving average filter. To remove filter bandwidth.
#' @param freq Cut-off frequency of Buttersworth filter
#' @details no details
#' @importFrom mFilter bwfilter
#' @export demod
demod <- function( x, # Time series to be demodulated
dt, # Time series time step
w1 = 2*pi/1, # Carrier frequency
filttype=c("MA","bw"), # Filter type
Nf, # Moving average filter width
nf = 1, # Number of moving averaging filter passes
freq) # Cut-off frequency of Buttersworth filter
{
filttype <- match.arg(filttype)
if (Nf>length(x)) {
warning("Demodulation filter longer than time series")
return(rep(NA, length(x)))
}
if (inherits(dt, "difftime")) dt <- unclass(dt)
times <- (1:length(x))*dt
demodfac <- exp(-1i*w1*times)
xdemod <- (x - mean(x))*demodfac
if (filttype=="MA"){
demodfiltRe <- Re(xdemod)
demodfiltIm <- Im(xdemod)
for (i in 1:nf){
demodfiltRe <- filter(demodfiltRe,rep(1/Nf,Nf))
demodfiltIm <- filter(demodfiltIm,rep(1/Nf,Nf))
}
demodfilt <- demodfiltRe + demodfiltIm*1i
} else {
demodfilt <- bwfilter(xdemod,freq=freq)
demodfilt <- demodfilt$trend
}
return(2*demodfilt)
}
#'
#' @title Calculate GPP(t) from by complex demodulation of O2 time series
#' @description Calculate GPP(t) from by complex demodulation of O2 time series
#' @param x Regularly sampled O2 time series. x can be a dataframe with time (POSIXt) in first column and O2 concentrations in second column
#' or a vector of concentrations sampled with time step = dt
#' @param dt Sampling time step. Either a difftime object or a numerical value.
#' When \code{dt} is given as a numerical value, the unit is assumed to be
#' days, unless \code{units} is provided. Can be omitted when \code{x} is a data frame, then \code{dt} is calculated from the time spacing between the first two samples.
#' @param units Unit of sampling time step.
#' @param Nfilt Moving average filter width for detrending time series. See \code{\link{GPPFourierPreprocess}}
#' @param NLowPass Moving average filter width for low pass filtering demodulated O2 series. See \code{\link{demod}}
#' @param fDL Optional. Relative fraction of light hours during the day. When \code{x} is a data frame, and no value for fDL is
#' provided, fDL is calculated with \code{\link{fDLfun}}
#' @param MAcorrect correction factor for moving average filter
#' @param nf Number of iterative applications of moving average filter, to reduce filter width
#' @param phi Latitude
#' @param lambda Longitude
#' @param ... Other parameters to be passed to GPPFourierPreprocess
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier}}, \code{\link{WindowGPPFourier.gts}}, \code{\link{SunRiseSet}}
#' @return Data frame with time and GPPt column.
# #' @details no details
#' @examples
#' # Calculate GPP(t) by complex demodulation of simulated water column time series
#' # More examples in vignette("GPPFourier").
#'
#' dt <- as.numeric(difftime(watercolumn$time[2], watercolumn$time[1] , units="days"))
#' Nfilt <- 1/dt
#' GPPt <- GPPFourier_t(watercolumn[,c("time", "O2")],
#' dt=dt,
#' units="days",
#' Detrend=TRUE,
#' filter=TRUE,
#' Nfilt=Nfilt,
#' NLowPass=Nfilt,
#' fDL=NULL,
#' circular=FALSE,
#' sides=2,
#' nf=1)
#'
#'
#' par(mfrow=c(2,1), cex=1.2)
#' plot(watercolumn[,c("time", "O2")], type="l", xlab="", ylab=expression(paste(O[2], " [", mu, "M]")))
#' title(main="Water column")
#' plot(GPPt, type="l", lwd=3, ylim=c(0,30), xlab="", ylab="GPP")
#' lines(watercolumn[,c("time","GPP")], col="red")
#' legend("topleft",
#' lty=1,
#' col=c("red", "black"),
#' legend=c( "Simulated GPP", expression(paste("Complex demodulated ", O[2], " series"))),
#' bty="n")
#'
#'
#' @references Cox et al (2017). Tune in on 11.57 muHz and listen to primary production. Biogeosciences Discussions. doi:10.5194/bg-2017-81
#' @export GPPFourier_t
GPPFourier_t <- function( x,
dt=NULL, # Sample time step
units=c("days", "hours", "mins", "secs"), # Time units of sample time step
Nfilt=NULL, # Moving average filter width for detrending time series
nf=1, # Number of moving average filter passes
NLowPass = NULL, # Moving average filter width for low pass filtering demodulated signal
fDL=NULL, # Relative fraction of light hours during the day
MAcorrect = FALSE, # Correction factor for moving average filter
phi= 51.176,
lambda=4.326,
...) # Other parameters to be passed to O22GPpreprocess
{
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (inherits(t, "POSIXt")){
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
} else {
if (any(0!=round(diff(diff(t)), digits=5))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- diff(c(t[1],t[2]))
}
if (is.null(fDL)){
if (is.null(phi)|is.null(lambda)) {stop("Provide longitude and lattitude")} else {
fDL <- fDLfun(t, phi=phi, lambda=lambda)
fDLfilt <- filter(fDL, rep(1/Nfilt, Nfilt))
}
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
O2 <- ts[,2]
}
} else {
if (is.null(dt)) stop("Provide dt")
if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
if (any(is.na(ts))) stop("NAs not allowed")
O2 <- ts
fDL <- rep(fDL, length(O2))
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
if (inherits(dt, "difftime")){
dtnumeric <- as.double(dt, units="days")
} else {
dtnumeric <- dt
}
if (is.null(Nfilt)) stop("Provide moving average filter width (Nfilt)")
if (Nfilt>length(O2)) stop("Filter longer than time series")
# angular frequency (cycles per day) associated with diurnal periodicity, i.e. T = 1 day
w <- 2*pi/1
xpreprocessed <- GPPFourierPreprocess(O2, dt=dt, Nfilt=Nfilt,...)
xfilt <- xpreprocessed$filt
xres <- xpreprocessed$res
timessub <- t[xpreprocessed$indices]
fDLsub <- fDL[xpreprocessed$indices]
# demodulate time series
xdemod <- demod(xres,dt=dtnumeric,Nf=NLowPass, nf=nf)
# correct for truncated sinusoid form of GPP(t)
theta1 <- pi*(1-fDLsub)
truncFac <- (sin(theta1) + (pi-theta1) * cos(theta1))/(pi-theta1 + 0.5 * sin(2*theta1))
if (MAcorrect){
# correct for moving average frequency response
L <- Nfilt
maxfreq <- (length(xdemod) - 1)/(length(xdemod))/dt
w24h <- 2*pi/maxfreq
MAcorrection <- 1/(1-Hw(w24h,L))
} else{
MAcorrection <- 1
}
# collect all correction factors
xdemodcorrect <- xdemod*truncFac*Mod(MAcorrection)
# time derivative in the frequency domain
GPt <- w*Mod(xdemodcorrect)
return(data.frame(times=timessub,GPPt=GPt))
}
Try the GPPFourier package in your browser
Any scripts or data that you put into this service are public.
GPPFourier documentation built on Sept. 22, 2017, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions GPPFourierPreprocess GPPFourier WindowGPPFourier.gts WindowGPPFourier demod GPPFourier_t
Documented in demod GPPFourier GPPFourierPreprocess GPPFourier_t WindowGPPFourier WindowGPPFourier.gts
#
# Core functions to calculate GPP from O2 time series making use of Fourier method
#
#'
#' @title Pre-process O2 data for GPP calculation
#' @description Removes longer term trend of time series and confine series to integer
#' number of days
#' @param x A vector containing consecutive O2 concentrations, sampled at regular intervals (time step = dt).
#' @param dt Sampling time step. Class {difftime}, or \code{numeric} when units are provided.
#' @param units Units of sampling time step
#' @param Detrend Toggle time series detrending
#' @param filter Toggle time series filtering
#' @param Nfilt Moving average filter width (see documentation of
#' \code{filter()})
#' @param circular Moving average boundary condition (see documentation of
#' \code{filter()})
#' @param sides Moving average central or one sided (see documentation of
#' \code{filter()})
#' @details GPPFourierPreprocess() is called by GPPFourier() to pre-process the time series before calculating the Fourier amplitude at diel frequency and estimate GPP.
#' @return A list with components \item{filt}{The long term trend}
#' \item{res}{The original series without the long term trend, confined to an
#' integer number of days} \item{indices}{The indices of the subset of the
#' series that are retained, after filtering and confining to an integer number
#' of days.}
#' @seealso \code{\link{GPPFourier}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier}},
#' \code{\link{WindowGPPFourier.gts}}
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#' @keywords utilities
#' @export GPPFourierPreprocess
GPPFourierPreprocess <- function(x, # O2 time series from which GPP is to be calculated
dt=NULL, # Sampling time step
units=c("days", "hours", "mins"), # Time units of sampling time step
Detrend=TRUE, # Toggle time series detrending
filter = TRUE, # Toggle time series filtering
Nfilt=NULL, # Moving average filter width (see documentation of filter())
circular = FALSE, # Moving average boundary condition (see documentation of filter())
sides = 2) # Moving average central or one sided (see documentation of filter())
{
# Removes longer term trends and confines time series to integer number of periods
#
#
unit <- match.arg(units)
if (is.null(dt)) stop("Provide dt")
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}}
if (filter&is.null(Nfilt)) stop("Provide filter length (Nfilt)")
if (Nfilt>length(x)) stop("Filter longer than time series")
dt <- as.double(dt, units="days")
# confine series to integer number of periods (T=1 day = 1/dt samples)
Nperiod <- floor(length(x)*dt)/dt
indices <- 1:Nperiod
x <- x[indices]
if (Detrend){
detrendx <- detrend(x, returntrend=TRUE)
x <- detrendx$detrended
trendx <- detrendx$trend
} else {
trendx <- rep(mean(x), length(x))
}
if (filter){
cfilt <- rep(1/Nfilt,Nfilt)
filtx <- filter(x,cfilt,sides=sides,circular=circular)
resx <- x - filtx
indices <- indices[!is.na(filtx)]
} else {
filtx <- rep(mean(x),length(x))
resx <- x - mean(x)
}
# confine filtered series to integer number of periods (T=1 day = 1/dt samples)
filtx <- filtx[!is.na(filtx)]
resx <- resx[!is.na(resx)]
Nperiod <- floor(length(filtx)*dt)/dt
indices <- indices[1:Nperiod]
filtx <- filtx[1:Nperiod]
resx <- resx[1:Nperiod]
return(list(filt = filtx, res=resx, trend=trendx, indices=indices))
}
#'
#' @title Calculate GPP from O2 time series
#'
#' @description This is the workhorse to calculate aquatic Gross Primary Production from
#' high frequency dissolved oxygen data. It will often be more convenient to use the higher
#' level functions \code{\link{WindowGPPFourier}} and
#' \code{\link{WindowGPPFourier.gts}}
#'
#' @param x Regularly sampled O2 time series. x can be a dataframe with time (POSIXt) in first column and O2 concentrations in
#' second column or a vector of concentrations sampled with time step = dt
#' @param dt Sampling time step. Not used when x is a data frame, then dt is calculated #' from the time spacing between
#' the first two samples. Either a difftime object or a numerical value. When dt is given as a numerical value, \code{units}
#' have to be provided.
#' @param units Unit of sampling time step
#' @param Nfilt Moving average filter width
#' @param fDL Relative fraction of light hours during the day. When \code{x} is a data frame, and no value for fDL is provided,
#' fDL is calculated with \code{\link{fDLfun}}
#' @param filtcorrect Logical controlling whether GPP estimate
#' is corrected for signal falsely removed by filtering
#' @param ... Other parameters to be passed to \code{\link{GPPFourierPreprocess}}
#' @param phi Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
#' @param lambda Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
#' @return Volume specific average gross primary production in units of oxgen as in \code{x}.
#' per time unit (unit of \code{dt}, given by \code{units}).
#' Specifically: if O2 concentrations are in mg/L, and \code{units = days}, calculated GPP has units mg/L/day.
#' If surface specific primary production (e.g. per m2) is desired, this number has to be multiplied by the depth of the water body.
#' @details no details
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier}},
#' \code{\link{WindowGPPFourier.gts}}
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#'
#' @keywords utilities
#' @examples
#' DO <- Kruibeke[Kruibeke$time>="2010-06-03"&Kruibeke$time<="2010-06-13",]
#' dt <- as.numeric(difftime(DO$time[2],DO$time[1],units="days"))
#' DL <- fDLfun(DO$time[1], phi=51.176, lambda=4.326)
#' GPP <- GPPFourier(DO$O2,dt=dt,Nfilt=100,fDL=DL)
#'
#' @export GPPFourier
GPPFourier <- function( x, # O2 time series from which GPP is to be calculated
dt=NULL, # Sampling time step
units=c("days","hours","mins"), # Time units of sampling time step
Nfilt=NULL, # Moving average filter width
fDL=NULL, # Relative fraction of light hours during the day
filtcorrect = FALSE, # Logical controlling whether GP estimate is corrected for signal removed by moving average filtering
phi=NULL, # Optional: Lattitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
lambda=NULL, # Optional: Longitude. Only used when fDL=NULL and daylength is calculated with \code{\link{fDLfun}}
...) # Other parameters to be passed to GPPFourierprocess
{
# Calculate average gross production from O2 time series
# Units: unit of oxgen as in ts per days
#
#
# Error handling
#
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (inherits(t, "POSIXt")){
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
} else {
if (any(0!=round(diff(diff(t)), digits=5))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- diff(c(t[1],t[2]))
}
if (is.null(fDL)){
if (is.null(phi)|is.null(lambda)) stop("Provide longitude and lattitude") else {
fDL <- fDLfun(t, phi=phi, lambda=lambda)
fDL <- mean(fDL) # average fDL of time series is used
}
}
O2 <- ts[,2]
}
} else {
if (is.null(dt)) stop("Provide dt")
if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
if (any(is.na(ts))) stop("NAs not allowed")
O2 <- ts
# fDL <- rep(fDL, length(O2))
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
# if (!inherits(dt,"difftime")){
# if (is.null(units)) {
# stop("Either provide dt as difftime object, or provide units of dt")
# } else {
# dt <- as.difftime(dt, units=unit)
# }}
# if (is.null(Nfilt)) stop("Provide moving average filter width (Nfilt)")
# if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
# if (as.numeric(length(x) < Nfilt)) stop("Filter longer than time series")
# angular frequency (rad per day) associated with diurnal periodicity, i.e. T = 1 day
w <- 2*pi/1
# Remove longer term trends and confines time series to integer number of periods
xpreprocessed <- GPPFourierPreprocess(x, dt=dt, Nfilt=Nfilt, ...)
filtx <- xpreprocessed$filt
resx <- xpreprocessed$res
indices <- xpreprocessed$indices
dt <- as.double(dt,units="days")
# direct calculation of A24
times <- (1:length(resx))*dt
A24 <- 2*sum(exp(-times*1i*w)*resx)/length(resx)
# correction for truncated sinusoidal form of GPP
theta1 <- pi*(1-fDL)
truncFac <- (sin(theta1) + (pi-theta1) * cos(theta1))/(pi-theta1 + 0.5 * sin(2*theta1))
A24 <- A24*truncFac
# Calculate time derivative in frequency domain
A24 <- 1i*w*A24
# Correct for frequency response of moving average filter
if (filtcorrect){
if (is.null(Nfilt)){
warning("Nfilt=NULL. No filter correction calculated")
} else {
L <- Nfilt
maxfreq <- (length(resx) - 1)/(length(resx))/dt
w24h <- 2*pi/maxfreq
A24 <- A24/(1-Hw(w24h,Nfilt))
}}
return(Mod(A24))
}
#'
#' @title Calculate GPP from O2 time series with gaps, in consecutive time-windows of prescribed
#' length
#'
#' @description Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length. High level function to handle time series with gaps. Short gaps are
#' interpolated, series is split at large gaps.
#'
#'
#' @param x Dataframe containing time (POSIXt) in first column and O2 concentrations in the second column.
#' @param dt Sampling time step. Either a difftime object or a numerical value.
#' When \code{dt} is given as a numerical value, \code{units} have to be provided (defaults to \code{days}). When omitted \code{dt} is calculated from the time spacing between the first two samples in the data frame.
#' @param units Unit of sampling time step
#' @param gapMaxN Minimum number of missing data points to be considered a gap in the time series. Gaps smaller than gapMaxN are interpolated. Defaults to the number of samples corresponding to a 4 hour gap.
#' @param Width Width [days] of the time-windows for which GPP is calculated
#' @param filtWidth Length of moving average filter [hours] to filter O2
#' @param ... Other parameters to be passed to WindowGPPFourier()
#' @return Average gross primary production in units of oxgen as in \code{x}
#' per day.
#' @importFrom Tides gapsts
#' @details This high level function analyses multiple regular O2 time series at once, each separated by a gap larger than gapMaxN (the number of . This is convenient e.g. for regular O2 data from a single location, where some data is missing due to sensor maintenance or replacement. \code{x} is assumed to contain all these data, pasted together in a single data frame (e.g. by \code{\link{cbind}}ing all individual regular time series or just reading from a single data file). Additionally, missing values within the regular series are allowed and are interpolated with \code{\link{gapfill}}. Uses \code{\link{WindowGPPFourier}}.
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier_t}}, \code{\link{SunRiseSet}},
#' \code{\link{WindowGPPFourier}}
#' @keywords utilities
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#' @examples
#' par(mfrow=c(2,1))
#' plot(Kruibeke, pch=20,xlim=as.POSIXct(c("2010-01-01", "2010-12-31")))
#' GPPAll_4 <- WindowGPPFourier.gts(Kruibeke,
#' gapMaxN = 10,
#' Width = 10,
#' filtWidth=1*24,
#' phi=51.176,
#' lambda=4.326,
#' Detrend=TRUE,
#' filter=TRUE,
#' filtcorrect=FALSE)
#' plot(GPPAll_4$time,GPPAll_4$GPP*9.3,
#' col="black",
#' pch=19,
#' type="b",
#' xlab="time",
#' ylab="GPP",
#' xlim=as.POSIXct(c("2010-01-01", "2010-12-31")))
#'
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @export WindowGPPFourier.gts
WindowGPPFourier.gts <- function(x,
dt = difftime(x[2,1],x[1,1],units="days"),
units = c("days", "hours", "mins"),
gapMaxN = 4/24/dt,
Width = 14,
filtWidth=16,
...) {
# x: data frame with first column t, and second column O2
# gapMaxN: number of missing samples that will be interpolated. More than gapMaxN missing samples will be handled as a gap
# with dt = 10 min, N = 6 amounts to 1h-gaps being linearly interpolated, 6*4 = 4h gaps interpolated
ts <- x # internal alias
ts$O2interpol <- ts[,2]
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
dt <- as.double(dt, units="days")
# Remove all gaps > dt*gapMaxN
gaps2 <- gapsts(ts[,1], dtMax=dt*gapMaxN, unit="days")
ts$N <- 0
ts$N[is.element(ts[,1],gaps2$t2)] <- 1
ts$N <- cumsum(ts$N)
# interpolate all gaps < dt*gapsMaxN
tslist <- split(ts[,1:2], ts$N)
tsfilledAll <- lapply(tslist, gapfill)
tsfilled <- do.call(rbind, tsfilledAll)
# re-add N
gaps2 <- gapsts(tsfilled[,1], dtMax=dt*gapMaxN, unit="days")
tsfilled$N <- 0
tsfilled$N[is.element(tsfilled$t,gaps2$t2)] <- 1
tsfilled$N <- cumsum(tsfilled$N)
# Windowed GP calculation on 'continuous' subsets of 10d or more, with small gaps filled
tslist <- split(tsfilled[,c("t","y")], tsfilled$N)
dtlist <- lapply(tslist, function(x) {x[1,1]- x[2,1]})
sublengths <- tapply(tsfilled$t, tsfilled$N, length)
tslist <- tslist[sublengths > Width/dt]
GPAll <- lapply(tslist, WindowGPPFourier, Width=Width, filtWidth=filtWidth, ...)
GPAll <- do.call(rbind, GPAll)
return(GPAll)
}
#'
#' @title Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length
#'
#' @description Calculate GPP from O2 time series in consecutive time-windows of prescribed
#' length.
#'
#' @param x Vector or dataframe containing consecutive O2 concentrations, sampled at regular intervals (time step = dt). In case of irregular time intervals, consider \code{\link{resample}}.
#' When x is a data frame, time (POSIXt) must be in first column and O2 in second column; when x is a vector, dt must be provided.
#' @param dt Sampling time step. Either a difftime object or a numerical value. When \code{dt} is given as a numerical value,
#' \code{units} have to be provided. Can be omitted when \code{x} is a data frame, then \code{dt} is calculated from the
#' time spacing between the first two samples.
#' @param units Unit of sampling time step
#' @param Width Width [days] of the time-windows for which GPP is calculated
#' @param Nblocks Number of consecutive subsections on which to calculate GPP. By default the number of windows of Width fitting into the range of x[,1]
#' Default the largest number of blocks ow width=Width that
#' @param phi Latitude of location where O2 series was recorded. Used to
#' calculate relative fraction of light hours during the day, during which
#' production takes place
#' @param lambda Longitude of location where O2 series was recorded
#' @param filtWidth [hours] Lenght of moving average filter to filter O2
#' series. See GPPFourierPreprocess()
#' @param fDL Optional. Relative fraction of daylight. When \code{x} is a data frame, fDL is calculated with \code{\link{fDLfun}}
#' @param ... Other parameters to be passed to \code{GPPFourier()}
#' @return Average gross primary production in units of oxgen as in \code{x}
#' per day.
#' @details no details
#' @author Tom Cox <tom.cox@uantwerp.be>
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier}}, \code{\link{GPPFourier_t}}, \code{\link{WindowGPPFourier.gts}}, \code{\link{SunRiseSet}}
#' @keywords utilities
#'
#' @references Cox T.J.S. et al. (2015) 'Estimating primary production from oxygen time series: a novel approach in the frequency domain', Limnology And Oceanography:Methods 13, 529-552. DOI: 10.1002/lom3.10046
#'
#'
#' @export WindowGPPFourier
WindowGPPFourier <- function(x,
dt=NULL,
units=c("days","hours","mins"),
Width=14,
Nblocks=floor(unclass(difftime(range(ts[,1])[2], range(ts[,1])[1],units="days"))[1]/Width),
fDL=NULL,
phi,
lambda,
filtWidth=16,
...){
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
O2 <- ts[,2]
if (is.null(fDL)){
fDL <- fDLfun(t, phi=phi, lambda=lambda)
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
}
} else {
if (is.null(dt)) stop("Provide dt")
if (any(is.na(ts))) stop("NAs not allowed")
if (is.null(fDL)) {
stop("Provide daylight fraction of day (fDL)")
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
O2 <- ts
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
dtnumeric <- as.double(dt, units="days")
Nfilt <- filtWidth/dtnumeric/24
blockst1 <- t[1] + (0:(Nblocks-1))*86400*Width
blockmids <- t[1] + Width/2*86400 + (0:(Nblocks-1))*86400*Width
blocknum <- approx(x=blockst1,y=0:(Nblocks-1),xout=t,method="constant", rule=2)$y
blocklength <- tapply(t, blocknum, length)
# calculate GPP on all consecutive subsets of Width
GPAll <- NULL
for (i in unique(blocknum))
{
O2i <- O2[blocknum==i]
fDLi <- mean(fDL[blocknum==i])
ti <- t[blocknum==i]
if (length(O2i) < 1.5*Nfilt){
GPAll <- c(GPAll,NA)
} else {
# RiseSet <- SunRiseSet(ti,phi=phi,lambda=lambda) + 2 #summertime = GMT+2
# DL <- diff(RiseSet)
# fDL <- DL/24
GPO2 <- GPPFourier(x=O2i,dt=dt,Nfilt= Nfilt, fDL=fDLi, ...) # Unit: mmol O2/m3/d
GPAll <- c(GPAll,GPO2)
}
}
return(data.frame(time=blockmids[unique(blocknum+1)],GPP=GPAll, t1=blockst1[unique(blocknum+1)], t2=blockst1[unique(blocknum+1)] + Width*86400))
}
#'
#' @title Complex demodulation
#' @description Retreive the time varying amplitude of a signal with fixed carrier frequency by complex demodulation
#' @param x Vector containing time series observations at constant sampling interval
#' @param dt Time step of time series. Arbitrary units. Class \code{Numeric} or \code{difftime}
#' @param w1 Carrier frequency. Inverse unit of \code{dt}.
#' @param filttype Filter type to be applied for low pass filtering. "MA" for moving average; "bw" for butterworth filter.
#' @param Nf Moving average filter width
#' @param nf Number of iterative applications of moving average filter. To remove filter bandwidth.
#' @param freq Cut-off frequency of Buttersworth filter
#' @details no details
#' @importFrom mFilter bwfilter
#' @export demod
demod <- function( x, # Time series to be demodulated
dt, # Time series time step
w1 = 2*pi/1, # Carrier frequency
filttype=c("MA","bw"), # Filter type
Nf, # Moving average filter width
nf = 1, # Number of moving averaging filter passes
freq) # Cut-off frequency of Buttersworth filter
{
filttype <- match.arg(filttype)
if (Nf>length(x)) {
warning("Demodulation filter longer than time series")
return(rep(NA, length(x)))
}
if (inherits(dt, "difftime")) dt <- unclass(dt)
times <- (1:length(x))*dt
demodfac <- exp(-1i*w1*times)
xdemod <- (x - mean(x))*demodfac
if (filttype=="MA"){
demodfiltRe <- Re(xdemod)
demodfiltIm <- Im(xdemod)
for (i in 1:nf){
demodfiltRe <- filter(demodfiltRe,rep(1/Nf,Nf))
demodfiltIm <- filter(demodfiltIm,rep(1/Nf,Nf))
}
demodfilt <- demodfiltRe + demodfiltIm*1i
} else {
demodfilt <- bwfilter(xdemod,freq=freq)
demodfilt <- demodfilt$trend
}
return(2*demodfilt)
}
#'
#' @title Calculate GPP(t) from by complex demodulation of O2 time series
#' @description Calculate GPP(t) from by complex demodulation of O2 time series
#' @param x Regularly sampled O2 time series. x can be a dataframe with time (POSIXt) in first column and O2 concentrations in second column
#' or a vector of concentrations sampled with time step = dt
#' @param dt Sampling time step. Either a difftime object or a numerical value.
#' When \code{dt} is given as a numerical value, the unit is assumed to be
#' days, unless \code{units} is provided. Can be omitted when \code{x} is a data frame, then \code{dt} is calculated from the time spacing between the first two samples.
#' @param units Unit of sampling time step.
#' @param Nfilt Moving average filter width for detrending time series. See \code{\link{GPPFourierPreprocess}}
#' @param NLowPass Moving average filter width for low pass filtering demodulated O2 series. See \code{\link{demod}}
#' @param fDL Optional. Relative fraction of light hours during the day. When \code{x} is a data frame, and no value for fDL is
#' provided, fDL is calculated with \code{\link{fDLfun}}
#' @param MAcorrect correction factor for moving average filter
#' @param nf Number of iterative applications of moving average filter, to reduce filter width
#' @param phi Latitude
#' @param lambda Longitude
#' @param ... Other parameters to be passed to GPPFourierPreprocess
#' @seealso \code{\link{GPPFourierPreprocess}}, \code{\link{GPPFourier}}, \code{\link{WindowGPPFourier.gts}}, \code{\link{SunRiseSet}}
#' @return Data frame with time and GPPt column.
# #' @details no details
#' @examples
#' # Calculate GPP(t) by complex demodulation of simulated water column time series
#' # More examples in vignette("GPPFourier").
#'
#' dt <- as.numeric(difftime(watercolumn$time[2], watercolumn$time[1] , units="days"))
#' Nfilt <- 1/dt
#' GPPt <- GPPFourier_t(watercolumn[,c("time", "O2")],
#' dt=dt,
#' units="days",
#' Detrend=TRUE,
#' filter=TRUE,
#' Nfilt=Nfilt,
#' NLowPass=Nfilt,
#' fDL=NULL,
#' circular=FALSE,
#' sides=2,
#' nf=1)
#'
#'
#' par(mfrow=c(2,1), cex=1.2)
#' plot(watercolumn[,c("time", "O2")], type="l", xlab="", ylab=expression(paste(O[2], " [", mu, "M]")))
#' title(main="Water column")
#' plot(GPPt, type="l", lwd=3, ylim=c(0,30), xlab="", ylab="GPP")
#' lines(watercolumn[,c("time","GPP")], col="red")
#' legend("topleft",
#' lty=1,
#' col=c("red", "black"),
#' legend=c( "Simulated GPP", expression(paste("Complex demodulated ", O[2], " series"))),
#' bty="n")
#'
#'
#' @references Cox et al (2017). Tune in on 11.57 muHz and listen to primary production. Biogeosciences Discussions. doi:10.5194/bg-2017-81
#' @export GPPFourier_t
GPPFourier_t <- function( x,
dt=NULL, # Sample time step
units=c("days", "hours", "mins", "secs"), # Time units of sample time step
Nfilt=NULL, # Moving average filter width for detrending time series
nf=1, # Number of moving average filter passes
NLowPass = NULL, # Moving average filter width for low pass filtering demodulated signal
fDL=NULL, # Relative fraction of light hours during the day
MAcorrect = FALSE, # Correction factor for moving average filter
phi= 51.176,
lambda=4.326,
...) # Other parameters to be passed to O22GPpreprocess
{
ts <- x # internal alias
if (inherits(ts,"data.frame")){
if (any(is.na(ts[,c(1,2)]))) stop("NAs not allowed")
if (dim(ts)[2]>=2){
t <- ts[,1]
if (inherits(t, "POSIXt")){
if (any(0!=diff(diff(round.POSIXt(t, units="secs"))))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- difftime(t[2],t[1],units="days")
} else {
if (any(0!=round(diff(diff(t)), digits=5))) stop("Samples are not equidistant and/or gaps are present in the series. Re-sample, fill gaps and/or split series in continuous subsets")
dt <- diff(c(t[1],t[2]))
}
if (is.null(fDL)){
if (is.null(phi)|is.null(lambda)) {stop("Provide longitude and lattitude")} else {
fDL <- fDLfun(t, phi=phi, lambda=lambda)
fDLfilt <- filter(fDL, rep(1/Nfilt, Nfilt))
}
} else if (length(fDL)==1){
fDL <- rep(fDL, length(ts[,2]))
} else if (length(fDL)!= length(ts[,2])){
stop("When fDL is provided, it should either be single valued or a vector of same length as the time series")
}
O2 <- ts[,2]
}
} else {
if (is.null(dt)) stop("Provide dt")
if (is.null(fDL)) stop("Provide daylight fraction of day (fDL)")
if (any(is.na(ts))) stop("NAs not allowed")
O2 <- ts
fDL <- rep(fDL, length(O2))
unit <- match.arg(units)
if (!inherits(dt,"difftime")){
if (is.null(units)) {
stop("Either provide dt as difftime object, or provide units of dt")
} else {
dt <- as.difftime(dt, units=unit)
}
}
t <- 1:length(O2)*dt
}
if (inherits(dt, "difftime")){
dtnumeric <- as.double(dt, units="days")
} else {
dtnumeric <- dt
}
if (is.null(Nfilt)) stop("Provide moving average filter width (Nfilt)")
if (Nfilt>length(O2)) stop("Filter longer than time series")
# angular frequency (cycles per day) associated with diurnal periodicity, i.e. T = 1 day
w <- 2*pi/1
xpreprocessed <- GPPFourierPreprocess(O2, dt=dt, Nfilt=Nfilt,...)
xfilt <- xpreprocessed$filt
xres <- xpreprocessed$res
timessub <- t[xpreprocessed$indices]
fDLsub <- fDL[xpreprocessed$indices]
# demodulate time series
xdemod <- demod(xres,dt=dtnumeric,Nf=NLowPass, nf=nf)
# correct for truncated sinusoid form of GPP(t)
theta1 <- pi*(1-fDLsub)
truncFac <- (sin(theta1) + (pi-theta1) * cos(theta1))/(pi-theta1 + 0.5 * sin(2*theta1))
if (MAcorrect){
# correct for moving average frequency response
L <- Nfilt
maxfreq <- (length(xdemod) - 1)/(length(xdemod))/dt
w24h <- 2*pi/maxfreq
MAcorrection <- 1/(1-Hw(w24h,L))
} else{
MAcorrection <- 1
}
# collect all correction factors
xdemodcorrect <- xdemod*truncFac*Mod(MAcorrection)
# time derivative in the frequency domain
GPt <- w*Mod(xdemodcorrect)
return(data.frame(times=timessub,GPPt=GPt))
}
Try the GPPFourier package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.