#' Light-level geolocation data processing
#'
#' Basic data processing for data from light-level geolocation archival tags.
#' Provides facilities for importing and plotting data recorded by the
#' BAS tags, and in particular detecting times of twilight from the
#' recorded light data.
#'
#' @name TwGeos-package
#' @docType package
#' @author Simon Wotherspoon, Micheal Sumner and Simeon Lisovski.
NULL
## Data import
##' Read light data from a BAS tag
##'
##' The \code{readLig} function imports the light record from a ".lig"
##' file generated by a BAS tag.
##' @title Read Light Data
##' @param file the light file to import.
##' @param skip number of initial lines to skip
##' @return Returns a dataframe with columns
##' \item{\code{Valid}}{an indicator of whether the record is valid}
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Julian}}{date and time as Julian date}
##' \item{\code{Light}}{the recorded light level}
##' @importFrom utils read.csv
##' @export
readLig <- function(file,skip=0) {
## Read csv file and add column names
d <- read.csv(file,header=FALSE,skip=skip,
col.names=c("Valid","Date","Julian","Light"),
colClasses=c("character","character","numeric","integer"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%y %H:%M:%S",tz="GMT"))
d
}
##' Read the activity data from a BAS tag
##'
##' Functions for importing activity data from BAS tags. The
##' \code{readAct} function imports activity data from an early
##' generation BAS tag while \code{readAct2} imports activity data
##' from later generation BAS tags that run length encode activity.
##' If \code{readAct2} is passed a dataframe constructed with
##' \code{readLig}, it interpolates activity to the times at which
##' light is recorded and returns a merged dataframe with both light
##' and activity data.
##'
##' @title Read Activity Data
##' @param file the activity file to import.
##' @param d.lig dataframe of light data created with \code{readLig}
##' @param skip number of initial lines to skip
##' @return \code{readAct} returns a dataframe with columns
##' \item{\code{Valid}}{an indicator of whether the record is valid}
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Julian}}{date and time as Julian date}
##' \item{\code{Activity}}{the recorded activity level}
##' \item{\code{Temp}}{the recorded temperature}
##' If \code{readAct2} is passed a \code{d.lig} dataframe, it returns that
##' dataframe with an appended column
##' \item{\code{Activity}}{the number of seconds the tag has been wet.}
##' otherwise it returns a dataframe with columns
##' \item{\code{Valid}}{an indicator of whether the record is valid}
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Julian}}{date and time as Julian date}
##' \item{\code{Activity}}{the number of seconds of activity}
##' \item{\code{Wet}}{whether the tag was wet or dry}
##' @importFrom utils read.csv
##' @export
readAct <- function(file,skip=0) {
## Read csv file and add column names. Must explicitly specify
## column names as R cannot detect the correct number of columns to
## read
d <- read.csv(file,header=FALSE,skip=0,
col.names=c("Valid","Date","Julian","Activity","Temp"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%y %H:%M:%S",tz="GMT"))
## Add day number and time within day (as fraction of a day)
d$Day <- floor(d$Julian)
d$Time <- d$Julian-d$Day
d
}
##' @rdname readAct
##' @importFrom utils read.csv
##' @importFrom stats approx
##' @export
readAct2 <- function(file,d.lig=NULL,skip=0) {
## Read csv file and add column names. Must explicitly specify
## column names as R cannot detect the correct number of columns to
## read
d <- read.csv(file,header=FALSE,skip=0,
col.names=c("Valid","Date","Julian","Activity","Wet"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%y %H:%M:%S",tz="GMT"))
## Add day number and time within day (as fraction of a day)
d$Day <- floor(d$Julian)
d$Time <- d$Julian-d$Day
## Run length decode activity
if(!is.null(d.lig)) {
d <- d[d$Valid=="ok",]
date <- as.numeric(d$Date[1])+cumsum(c(0,d$Activity))
wet <- cumsum(c(0,ifelse(d$Wet=="wet",d$Activity,0)))
d.lig$Activity <- c(0,diff(approx(date,wet,as.numeric(d.lig$Date),rule=2)$y))
d <- d.lig
}
d
}
##' Read temperature data from a BAS tag
##'
##' The \code{readTem2} function reads temperature data from a later
##' generation BAS tag. If \code{readTem2} is passed a dataframe
##' constructed with \code{readLig}, it matches the recorded
##' temperatures to the closest recorded light reading and returns the
##' merged dataframe with both temperature and activity data.
##'
##' @title Read Temperature Data
##' @param file the temperature file to import.
##' @param d.lig dataframe of light data created with \code{readLig}
##' @param skip number of initial lines to skip
##' @return If \code{readTem2} is passed a "lig" dataframe, it returns that
##' dataframe with the appended column
##' \item{\code{Temperature}}{the recorded temperature}
##' otherwise it returns a dataframe with columns
##' \item{\code{Valid}}{an indicator of whether the record is valid}
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Julian}}{date and time as Julian date}
##' \item{\code{Temperature}}{the recorded temperature.}
##' @importFrom utils read.csv
##' @importFrom stats approx
##' @export
readTem2 <- function(file,d.lig=NULL,skip=0) {
## Read csv file and add column names. Must explicitly specify
## column names as R cannot detect the correct number of columns to
## read
d <- read.csv(file,header=FALSE,skip=0,
col.names=c("Valid","Date","Julian","Temp"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%y %H:%M:%S",tz="GMT"))
## Add day number and time within day (as fraction of a day)
d$Day <- floor(d$Julian)
d$Time <- d$Julian-d$Day
## Run length decode activity
if(!is.null(d.lig)) {
d <- d[d$Valid=="ok",]
k <- approx(d.lig$Date,1:nrow(d.lig),d$Date,method="constant",rule=2)$y
d.lig$Temp <- rep(NA,nrow(d.lig))
d.lig$Temp[k] <- d$Temp
d <- d.lig
}
d
}
##' Read light data from a Migrate Technologies tag
##'
##' The \code{readMTlux} function imports the light record from a
##' ".lux" file generated by a Migrate Technologies tag. These files
##' have an extensive header, and it may be necessary to set the
##' \code{skip} argument to skip the lines including the two column
##' headers.
##'
##' @title Read Light Data
##' @param file the lux file to import.
##' @param skip number of initial lines to skip
##' @return Returns a dataframe with columns
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Light}}{the recorded light level}
##' @importFrom utils read.table
##' @export
readMTlux <- function(file,skip=20) {
## Read csv file and add column names
d <- read.table(file,header=FALSE,skip=skip,
sep="\t",
col.names=c("Date","Light"),
colClasses=c("character","numeric"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%Y %H:%M:%S",tz="GMT"))
d
}
##' Read SST data from a Migrate Technologies tag
##'
##' The \code{readMTsst} function reads temperature data from a
##' Migrate Technologies tag. If \code{readMTsst} is passed a
##' dataframe constructed with \code{readMTlux}, it matches the
##' recorded temperatures to the closest recorded light reading and
##' returns the merged dataframe with both temperature and activity
##' data.
##'
##' @title Read Temperature Data
##' @param file the sst file to import.
##' @param d.lux dataframe of light data created with \code{readMTlux}
##' @param skip number of initial lines to skip
##' @return If \code{readMTsst} is passed a "lux" dataframe, it returns that
##' dataframe with the appended columns
##' \item{\code{MinSST}}{the minimum SST}
##' \item{\code{MaxSST}}{the maximum SST}
##' \item{\code{SST}}{the mean SST}
##' \item{\code{NSST}}{the number of SST samples}
##' otherwise it returns a dataframe with columns
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{MinSST}}{the minimum SST}
##' \item{\code{MaxSST}}{the maximum SST}
##' \item{\code{SST}}{the mean SST}
##' \item{\code{NSST}}{the number of SST samples}
##' @importFrom utils read.csv
##' @importFrom stats approx
##' @export
readMTsst <- function(file,d.lux=NULL,skip=20) {
## Read csv file and add column names
d <- read.table(file,header=FALSE,skip=skip,
sep="\t",
col.names=c("Date","MinSST","MaxSST","SST","NSST"),
colClasses=c("character","numeric","numeric","numeric","numeric"))
## Parse date
d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%Y %H:%M:%S",tz="GMT"))
## Run length decode activity
if(!is.null(d.lux)) {
k <- approx(d.lux$Date,1:nrow(d.lux),d$Date,method="constant",rule=2)$y
for(v in c("MinSST","MaxSST","SST")) {
d.lux[,v] <- rep(NA,nrow(d.lux))
d.lux[k,v] <- d[,v]
}
d.lux[,"NSST"] <- rep(0,nrow(d.lux))
d.lux[k,"NSST"] <- d[,"NSST"]
d <- d.lux
}
d
}
## Data export
##' To pre- process geolocator datasets with TAGS it may be necessary
##' to compress the data (i.e. remove repeated light levels).
##'
##' @title Export data to match TAGS requirements
##' @rdname export2TAGS
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param path if full path is specified, the output will be writen as a txt file ready to import into TAGS
##' @param file name of txt file (without suffix)
##' @param min ...
##' @param max ...
##' @return Return a data frame that matches rquirements for TAGS.
##' @export
export2TAGS <- function(tagdata, path = NULL, file = "exportTAGS", min = 0, max = 9999999999) {
light <- tagdata$Light
datetime <- tagdata$Date
light[light > max] <- max #trim high values to designated maximum
light[light < min] <- min #raise low values to designated minimum
l1 <- light[2:(length(light)-1)] #all but first and last light levels
l2 <- light[3:length(light)] #light levels subsequent to each in l1
l3 <- light[1:(length(light)-2)] #light levels previous to each in l1
keep <- c(TRUE, l1!=l2 | l1!=l3, TRUE) #True/False vector corresponding with lt. True means there is a change in light levels (keep these rows)
out <- data.frame(datetime = datetime[keep], light = light[keep]) #new dataframe with just the changing light levels
if(!is.null(path)) {
write.table(out, file = paste0(path, "/", file, ".txt"), sep = ",", quote = FALSE, row.names = FALSE, col.names = c("datetime", "light"))
}
out
}
##' Transforms twGeos twilight data frames into format reguiered by 'GeoLight'
##'
##' @title Transforms twiligth data to match GeoLight requirements
##' @rdname export2GeoLight
##' @param twilights a dataframe generated by \code{\link{findTwilights}}.
##' @return Return a data frame that matches GeoLight rquirements.
##' @export
export2GeoLight <- function(twilights) {
diff <- c(apply(cbind(twilights$Twilight[-nrow(twilights)], twilights$Twilight[-1]), 1, function(x) c(x[2]-x[1])/60/60), 0)
out <- data.frame(tFirst = as.POSIXct("1900-01-01 01:01",
"GMT"), tSecond = as.POSIXct("1900-01-01 01:01", "GMT"), type = 0, diff.max = 0)
rw <- 1
for (k in 1:(nrow(twilights) - 1)) {
if (as.numeric(difftime(twilights$Twilight[k], twilights$Twilight[k + 1])) < 24 & twilights$Twilight[k] != twilights$Twilight[k + 1]) {
out[rw, 1] <- twilights$Twilight[k]
out[rw, 2] <- twilights$Twilight[k + 1]
out[rw, 3] <- ifelse(twilights$Rise[k], 1, 2)
out[rw, 4] <- max(diff[k:(k+1)])
rw <- rw + 1
}
}
out <- subset(out, diff.max<23)
out[,-ncol(out)]
}
##' Utilities for manipulating hours
##'
##' Given a vector of POSIXct dates, \code{as.hour} extracts the time
##' of day component of the date and returns it as decimal hours.
##' Given a vector of decimal hours, \code{hourOffset} recodes the
##' decimal hour into a new 24 hour interval.
##' @title Hour Manipulation
##' @rdname hours
##' @param tm them timestamp as POSIXct.
##' @param hr the decimal hour to be wrap.
##' @param offset minimum hour of the interval to wrap into.
##' @return Return a decimal hour.
##' @examples
##' as.hour(as.POSIXct("2005-11-12 19:58:00"))
##' hourOffset(1:10,5)
##' @export
as.hour <- function(tm) {
(as.numeric(tm)-as.numeric(as.POSIXct(as.Date(tm))))/3600
}
##' @rdname hours
##' @export
hourOffset <- function(hr,offset=0) {
(hr-offset)%%24+(offset%%24)
}
##' Display time series data as an image.
##'
##' The \code{tsimage} function divides a sequence of samples into 24
##' hour periods, packs these into the columns of a matrix and
##' displays the result as an image. The data are resampled to ensure
##' the time interval between samples is regular and evenly divides
##' the 24 hour period. The time interval to which the data are
##' resampled is selected automatically, or can be specified with the
##' \code{dt} argument.
##'
##' The \code{tsimageLocator} function allows the user to select
##' pixels with the mouse pointer in the style of \code{locator}, and
##' returns the times corresponding to the selected pixels.
##'
##' These functions ignore clock drift and assume the sampling period
##' exactly divides 24 hours.
##'
##' @title Display Tag Data as an Image
##' @param date the sequence of sample times as POSIXct.
##' @param y the sequence of responses.
##' @param offset the starting hour for the vertical axis.
##' @param dt the time interval to which the data are resampled (secs).
##' @param xlab the x axis label.
##' @param ylab the y axis label.
##' @param xaxt a character that specifies the x axis type (see \code{par}).
##' @param ... additional arguments to be passed to \code{image}.
##' @param im the axes coordinates returned by \code{tsimage}.
##' @param n the number of points to select.
##' @return
##' \code{tsimage} returns the date and hour coordinates of the image
##' grid.
##' \code{tsimageLocator} returns the times corresponding to the
##' selected pixels.
##' @importFrom stats approx
##' @importFrom graphics axis axis.POSIXct box image
##' @export
tsimage <- function(date,y,offset=0,dt=NA,xlab="Date",ylab="Hour",xaxt=par("xaxt"),...) {
## Reasonable resampling intervals (secs)
dts <- c(5, 10, 15, 20, 30, 60, 90, 120, 180, 240,
300, 360, 400, 480, 540, 600, 720, 900, 960, 1200)
if(is.na(dt)) {
## Estimate resampling interval
dt <- mean(diff(as.numeric(date)))
dt <- dts[which.min(abs(dt-dts))]
}
## Get range to resample into
tmin <- .POSIXct(as.POSIXct(as.Date(date[1]))+offset*60*60,"GMT")
if(as.numeric(tmin) > as.numeric(date[1])) tmin <- tmin-24*60*60
tmax <- .POSIXct(as.POSIXct(as.Date(date[length(date)]))+offset*60*60,"GMT")
if(as.numeric(tmax) < as.numeric(date[length(date)])) tmax <- tmax+24*60*60
## Resample data
y <- approx(as.numeric(date),y,seq(as.numeric(tmin)+dt/2,as.numeric(tmax)-dt/2,dt))$y
## Restructure as mxn matrix
m <- 24*60*60/dt
n <- length(y)/m
y <- matrix(y,m,n)
## Show as image, with hours wrapped to [0,24].
hour <- seq(offset,offset+24,length=m+1)
day <- seq(tmin,tmax,length=n+1)
image(as.numeric(day),hour,t(y),axes=FALSE,xlab=xlab,ylab=ylab,...)
if(xaxt!="n") axis.POSIXct(1,day)
axis(2,at=seq(0,48,by=4),labels=seq(0,48,by=4)%%24)
box()
invisible(list(date=day,hour=hour,offset=offset))
}
##' @rdname tsimage
##' @importFrom graphics locator
##' @export
tsimageLocator <- function(im,n=512) {
## Force evaluation of im.
im <- im
## Select pixels and convert to date/time.
loc <- locator(n)
im$date[.bincode(loc$x,as.numeric(im$date),right=FALSE)]+3600*(loc$y-min(im$hour))
}
##' Plot times as hour vs date
##'
##' Analogs of \code{plot}, \code{points} and \code{lines} that plot
##' the time of day against date. The \code{tsimageRibbon} uses
##' \code{polygon} to produce a polygonal "ribbon" that represents a
##' series of daily time intervals.
##'
##' @title tsimage plot
##' @param date times as a vector of POSIXct.
##' @param date1 lower bound as a vector of POSIXct.
##' @param date2 upper bound as a vector of POSIXct.
##' @param offset the starting hour for the vertical axes.
##' @param xlab the label for the x axis.
##' @param ylab the label for the y axis.
##' @param ... additional arguments parameters to pass to \code{plot},
##' \code{points} or \code{lines}.
##' @importFrom graphics plot
##' @export
tsimagePlot <- function(date,offset,xlab="Date",ylab="Hour",...) {
plot(date,hourOffset(as.hour(date),offset%%24),xlab=xlab,ylab=ylab,...)
}
##' @rdname tsimagePlot
##' @importFrom graphics points
##' @export
tsimagePoints <- function(date,offset,...) {
points(date,hourOffset(as.hour(date),offset%%24),...)
}
##' @rdname tsimagePlot
##' @importFrom graphics lines
##' @export
tsimageLines <- function(date,offset,...) {
## Unfold so that lines do not cross from top to bottom
hour <- hourOffset(as.hour(date), offset)
hr <- hour[!is.na(hour)]
hour[!is.na(hour)] <- cumsum(c(hr[1], (diff(hr) + 12)%%24 - 12))
## Plot 24 hour translations of the unfolded line
lines(date,hour-24,...)
lines(date,hour+24,...)
lines(date,hour,...)
}
##' @rdname tsimageDeployment
##' @importFrom graphics lines
##' @export
tsimageDeploymentLines <- function(date, lon, lat, offset, zenith = 96, ...) {
tm <- seq(min(date), max(date), by = "day")
rise <- rep(c(TRUE, FALSE), length(tm))
c.dat <- data.frame(Twilight = twilight(rep(tm, each = 2), lon = lon, lat = lat,
rise = rise, zenith = zenith), Rise = rise)
tsimageLines(c.dat$Twilight[c.dat$Rise], offset = offset, ...)
tsimageLines(c.dat$Twilight[!c.dat$Rise], offset = offset, ...)
} #test for warning: Undocumented code objects; all functions need to have at least on documentation entry
##' @rdname tsimagePlot
##' @importFrom graphics polygon
##' @export
tsimageRibbon <- function(date1,date2,offset,...) {
hour1 <- hourOffset(as.hour(date1),offset)
hour2 <- hourOffset(as.hour(date2),offset)
hour2 <- ifelse(hour1 > hour2,hour2+24,hour2)
xs <- c(date1,rev(date2))
ys <- c(hour1,rev(hour2))
polygon(xs,ys-24,...)
polygon(xs,ys,...)
polygon(xs,ys+24,...)
}
##' Display light series as an image.
##'
##' The \code{lightImage} function displays sequence of light
##' samples into 24 hour periods, packs these into the columns of a
##' matrix and displays the result as an image.
##'
##' This function is essentially a wrapper for \code{tsimage}.
##' @title Display Light Data as an Image
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param offset the offset for the vertical axis in hours.
##' @param zlim the range of light levels to plot.
##' @param dt the time interval to which the data are resampled (secs).
##' @param xlab the x axis label.
##' @param ylab the y axis label.
##' @param col the colour palette
##' @param ... additional arguments to pass to image.
##' @return Returns the date and hour coordinates of the image grid.
##' @importFrom grDevices grey
##' @export
lightImage <- function(tagdata,offset=0,zlim=c(0,64),dt=NA,
xlab="Date",ylab="Hour",col=grey(seq(0,1,length=256)),...) {
if(is.na(zlim[1])) zlim[1] <- min(tagdata$Light,na.rm=T)
if(is.na(zlim[2])) zlim[2] <- max(tagdata$Light,na.rm=T)
tsimage(tagdata$Date,pmax(zlim[1],pmin(zlim[2],tagdata$Light)),
offset=offset,dt=dt,zlim=zlim,col=col,xlab=xlab,ylab=ylab,...)
}
##' Search for pairs of twilights spanning night.
##'
##' Search for sunset, sunrise pairs that correspond to a given light
##' threshold.
##'
##' Given a set of times (\code{include}) known to fall in the night,
##' \code{findTwilights} determines the twilights that span these times, and
##' computes the corresponding midnights. It then searches for periods
##' of darkness that lie approximately 24 hours from these midnights,
##' repeating the process until no new twilight pairs are found.
##'
##' If \code{interleave=TRUE}, the sunrise and sunset times are
##' interleaved andreturned as a single sequence of twilights,
##' otherwise sunset and sunrise times are returned separately. The
##' function \code{interleave.twilights} takes a dataframe of separate
##' sunset and sunrise times and interleaves them to form a sequence
##' of twilight times.
##'
##' @title Search for twilight times
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param threshold the light threshold that defines twilight.
##' @param include a vector of times as POSIXct. Nights that span these
##' times are included in the search.
##' @param exclude a vector of POSIXct times. Nights that span these
##' times are excluded from the search.
##' @param extend a time in minutes. The function seeks periods of
##' darkness that differ from one another by 24 hours plus or minus
##' this interval.
##' @param dark.min a time in minutes. Periods of darkness shorter
##' than this interval will be excluded.
##' @return A dataframe with columns
##' \item{\code{Twilight}}{times of twilight}
##' \item{\code{Rise}}{logical indicating sunrise}
##' where each row corresponds to a single twilight.
##' @export
findTwilights <- function(tagdata,threshold,include,
exclude=NULL,extend=0,dark.min=0) {
## Extract date and light data
date <- tagdata$Date
light <- tagdata$Light
## Is any x in each [a,b]
containsAny <- function(a,b,x) {
f <- logical(length(a))
for(k in seq_along(x))
f <- f | (a <= x[k] & b >= x[k])
f
}
## Convert to minutes
extend <- 60*extend
dark.min <- 60*dark.min
## Calculate intervals [a,b] of darkness
## a is first points before light drops below threshold
## b is first points before light rises to or above threshold
l <- (light>=threshold)
f <- diff(l)
a <- which(f==-1)
b <- which(f==1)
## Keep only fall-rise pairs
if(b[1] < a[1]) b <- b[-1]
a <- a[1:length(b)]
## Only keep intervals that do not include excluded data and are
## less than 24 hours in length
keep <- (!containsAny(date[a],date[b+1],exclude) &
(as.numeric(date[b+1])-as.numeric(date[a]) < 86400) &
(as.numeric(date[b+1])-as.numeric(date[a]) > dark.min))
a <- a[keep]
b <- b[keep]
## Compute bounding dates and midpoint
aDate <- date[a]
bDate <- date[b+1]
mDate <- aDate+(bDate-aDate)/2
## Iteratively expand set of twilights by searching for additional
## twilights +/- 24 hrs from midpoints of existing set.
keep <- logical(length(a))
add <- containsAny(aDate,bDate,include) & !keep
while(any(add)) {
keep <- keep | add
mid <- c(mDate[add]-86400,mDate[add]+86400)
add <- containsAny(aDate-extend,bDate+extend,mid) & !keep
}
a <- a[keep]
b <- b[keep]
## Interpolate times to get exact twilights.
ss <- date[a]+(threshold-light[a])/(light[a+1]-light[a])*(date[a+1]-date[a])
sr <- date[b]+(threshold-light[b])/(light[b+1]-light[b])*(date[b+1]-date[b])
data.frame(Twilight=.POSIXct(as.vector(t(cbind(ss,sr))),"GMT"),
Rise=rep(c(F,T),length(ss)))
}
##' Calibration functon.
##'
##' ...
##'
##' @title Search for twilight times
##' @param twilight ...
##' @param rise ...
##' @param lon ...
##' @param lat ...
##' @param method ...
##' @param plot ...
##' @return a vector with the _zero_ and the _median_ zenith angle as well as the parameters of the error distribution.
##' @export
thresholdCalibration <- function(twilight, rise, lon, lat, method = "log-normal", plot= TRUE) {
if (!(method %in% c("gamma", "log-norm")))
stop("Method can only be `gamma` or `log-norm`.")
tab <- data.frame(Twilight = twilight, Rise = rise)
sun <- solar(tab[, 1])
z <- refracted(zenith(sun, lon, lat))
inc = 0
repeat {
twl_t <- SGAT::twilight(tab[, 1], lon, lat, rise = tab[, 2], zenith = max(z) + inc, iters = 5, closest = TRUE)
twl_dev <- ifelse(tab$Rise, as.numeric(difftime(tab[, 1], twl_t, units = "mins")),
as.numeric(difftime(twl_t, tab[, 1], units = "mins")))
if (all(twl_dev >= 0) | inc>1.5) {
break
}
else {
inc <- inc + 0.01
}
}
OutS <- rep(TRUE, length(twl_dev))
if(any(abs(twl_dev)>12*60)) OutS[abs(twl_dev)>12*60] <- FALSE
twl_dev[!OutS] <- NA
twl_t[!OutS] <- NA
z0 <- max(z) + inc
seq <- seq(0, max(twl_dev, na.rm = T), length = 100)
if (method == "log-norm") {
fitml_ng <- suppressWarnings(fitdistr(twl_dev[!is.na(twl_dev)], "log-normal"))
lns <- dlnorm(seq, fitml_ng$estimate[1], fitml_ng$estimate[2])
}
if (method == "gamma") {
fitml_ng <- suppressWarnings(fitdistr(twl_dev[!is.na(twl_dev)], "gamma"))
lns <- dgamma(seq, fitml_ng$estimate[1], fitml_ng$estimate[2])
}
diffz <- as.data.frame(cbind(min = apply(cbind(tab[, 1], twl_t), 1,
function(x) abs(x[1] - x[2]))/60, z = z))
mod <- lm(z ~ min, data = diffz)
mod2 <- lm(min ~ z, data = diffz)
z1 <- median(z)
if (plot) {
opar <- par(mar = c(10, 4, 1, 1))
if (method == "log-norm")
hist(twl_dev, freq = F, breaks = 26, main = "Twilight Model (log-norm)",
xlab = "twilight error (min)")
if (method == "gamma")
hist(twl_dev, freq = F, breaks = 26, main = "Twilight Model (gamma)",
xlab = "twilight error (min)")
lines(seq, lns, col = "firebrick", lwd = 3, lty = 2)
points(predict(mod2, newdata = data.frame(z = z0)), 0,
pch = 21, cex = 5, bg = "white", lwd = 2)
text(predict(mod2, newdata = data.frame(z = z0)), 0,
"0")
points(predict(mod2, data.frame(z = z1)), 0, pch = 21, cex = 5, bg = "white", lwd = 2)
text(predict(mod2, data.frame(z = z1)), 0, "1")
axis(1, at = seq(0, max(twl_dev, na.rm = T), 6), labels = round(90 -
predict(mod, newdata = data.frame(min = seq(0, max(twl_dev, na.rm = T),
6))), 1), line = 5)
mtext("sun elevation angle (degrees)", 1, line = 8)
if (method == "log-norm")
legend("topright", paste(c("0. Zenith angle (zero)",
"1. Zenith angle (median)", "log-mean", "log-sd"),
round(c(z0, z1, fitml_ng$estimate[1], fitml_ng$estimate[2]),
3)), bty = "n")
if (method == "gamma")
legend("topright", paste(c("0. Zenith angle (zero)",
"1. Zenith angle (median)", "shape", "scale"),
round(c(z0, z1, fitml_ng$estimate[1], fitml_ng$estimate[2]),
3)), bty = "n")
par(opar)
}
c(z1 = z1, z0 = z0, log.mean = fitml_ng$estimate[1], log.sd = fitml_ng$estimate[2])
}
##' Search for crepuscular intervals
##'
##' Using the twilights calculated by \code{findTwilights},
##' \code{findCrepuscular} searchs for the intervals of reduced light
##' intensity around twilight and returns the dataframe constructed by
##' \code{findTwilights} augmented with the start and end times of
##' the crepuscular segment. The sugset of data corresponding to these
##' segments can be extracted with \code{extractCrepuscular}.
##'
##' @title Search for crepuscular light segments
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param twilights a dataframe generated by
##' \code{\link{findTwilights}}.
##' @param min.threshold the light level that defines the dark end of
##' the crepuscular period.
##' @param max.threshold the light level that defines the light end of
##' the crepuscular period.
##' @param extend.dark the number of hours to search towards midnight.
##' @param extend.light the number of hours to search towards noon.
##' @return A dataframe with columns.
##' \item{\code{Twilight}}{times of twilight}
##' \item{\code{Rise}}{logical indicating sunrise}
##' \item{\code{Start}}{the start of the crepuscular period}
##' \item{\code{End}}{the end of the crepuscular period}
##' where each row corresponds to a single twilight.
##' @export
findCrepuscular <- function(tagdata,twilights,
min.threshold=0,max.threshold=60,
extend.dark=1,extend.light=2) {
runlength <- function(x) {
y <- cumsum(x)
y-cummax(ifelse(!x,y,0))
}
## Extract date and light
date <- tagdata$Date
light <- tagdata$Light
twilight <- twilights$Twilight
rise <- twilights$Rise
start <- twilight
end <- twilight
ks <- findInterval(twilight,date)
## Loop over twilights
for(i in seq_along(twilight)) {
if(rise[i]) {
## Search forward to noon
k <- ks[i]
ls <- light[date > date[k] & date <= date[k]+extend.light*60*60]
end[i] <- date[k+min(which(ls >= max.threshold),length(ls))]
## Search backward to midnight
k <- ks[i]+1
ls <- light[date <= date[k] & date >= date[k]-extend.dark*60*60]
rs <- rev(runlength(ls<=min.threshold))
start[i] <- date[k+1-min(which(rs>=min(max(rs),4)))]
} else {
## Search backward to noon
k <- ks[i]+1
ls <- light[date < date[k] & date >= date[k]-extend.light*60*60]
start[i] <- date[k-min(which(rev(ls) >= max.threshold),length(ls))]
## Search forward to midnight
k <- ks[i]
ls <- light[date >= date[k] & date <= date[k]+extend.dark*60*60]
rs <- rev(runlength(rev(ls)<=min.threshold))
end[i] <- date[k-1+min(which(rs>=min(max(rs),4)))]
}
}
## Return
cbind(twilights,
Start=start,
End=end)
}
##' Extract crepuscular light segments
##'
##' Extract the crepuscular subsets of the light data, given a
##' dataframe of start and end points generated by
##' \code{\link{findCrepuscular}}.
##'
##' The \code{adjust.interval} argument can be used to specify a
##' timing adjustment for tags that report the maximum light level
##' observed in the preceeding sampling interval. If this argument is
##' zero, no adjustment is made, otherwise the timestamps of sunset
##' intervals will be adjusted by this interval to compensate for the
##' recording behaviour of the tag.
##'
##' If \code{filter.threshold} is non \code{NULL}, the light levels
##' are filtered to remove shaded observations. The light levels at
##' sunrise and sunset are assumed to be monotonic increasing, and the
##' running maximum in the direction of day is calculated. Light
##' levels that are \code{filter.threshold} units below the running
##' maximum are removed.
##'
##' @title Extract crepuscular light segments
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param twilights a dataframe generated by
##' \code{\link{findCrepuscular}}.
##' @param adjust.interval timing adjustment for sunset intervals.
##' @param filter.threshold \code{NULL} or a positive threshold.
##' @return Returns a dataframe with columns
##' \item{\code{Date}}{the date and time of the observation}
##' \item{\code{Light}}{the recorded light level}
##' \item{\code{Segment}}{an integer identifying the segment}
##' @export
extractCrepuscular <- function(tagdata,twilights,
adjust.interval=0,filter.threshold=NULL) {
date <- vector(mode="list",nrow(twilights))
lght <- vector(mode="list",nrow(twilights))
sgmt <- vector(mode="list",nrow(twilights))
indices <- as.numeric(factor(as.numeric(twilights$Twilight)))
## For each twilight
for(k in seq_len(nrow(twilights))) {
is <- which(tagdata$Date >= twilights$Start[k] &
tagdata$Date <= twilights$End[k])
## Remove non-increasing segments
if(!is.null(filter.threshold)) {
light <- tagdata$Light[is]
lower <- if(twilights$Rise[k]) cummax(light) else rev(cummax(rev(light)))
is <- is[light >= lower-filter.threshold]
}
## Timing offset
off <- if(!twilights$Rise[k]) adjust.interval else 0
date[[k]] <- tagdata$Date[is]-off
lght[[k]] <- tagdata$Light[is]
sgmt[[k]] <- rep(indices[k],length.out=length(is))
}
data.frame(Date=.POSIXct(unlist(date),"GMT"),
Light=unlist(lght),
Segment=unlist(sgmt))
}
##' Find intervals of darkness spanning more than 24 hours.
##'
##' The \code{\link{findTwilights}} function find twilight times by
##' searching for adjacent nights, but is unable to search beyond a
##' period of light or dark spaaning greater than 24 hours. This function
##' finds extended periods of light or dark that span more than 24 hours.
##' @title Find extended periods of light or dark
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param threshold the light threshold that defines twilight.
##' @return A dataframe with one row for each intervals and columns
##' \item{\code{Start}}{the start of the interval as POSIXct}
##' \item{\code{Start}}{the end of the interval as POSIXct}
##' \item{\code{Midpoint}}{the midpoint of the interval as POSIXct}
##' @export
findChunk24 <- function(tagdata,threshold) {
## Extract date and light data
date <- tagdata$Date
light <- tagdata$Light
## Calculate light/dark transitions
l <- (light>=threshold)
f <- diff(l)
ab <- which(abs(f)==1)
ab <- ab+f[ab]==1
## Find blocks longer than 24 hrs and return
tm <- as.numeric(date[ab])
keep <- which(diff(tm) >= 86400)
a <- ab[keep]
b <- ab[keep+1]
data.frame(Start=date[a],
End=date[b],
Midpoint=date[a] + (date[b+1]-date[a])/2)
}
##' Adjust twilight estimates for max light recording.
##'
##' Some tags record the maximum light level observd in the preceding
##' sampling interval, and this behaviour biases naive estimates of
##' the times of sunset. This function offsets sunset times by the
##' sampling interval to correct for this bias.
##' @title Twilight Adjustment
##' @param twilights dataframe of twilight times.
##' @param interval the sampling interval.
##' @param fixed logical vector indicating twilights to be held fixed.
##' @return a dataframe of bias correct twilight times.
##' @export
twilightAdjust <- function(twilights,interval,fixed=FALSE) {
adj <- !twilights$Rise & !fixed
twilights$Twilight[adj] <- twilights$Twilight[adj]-interval
twilights
}
##' Automated (Objective) editing and deletion of twilight times (experimental)
##'
##' This function searches for outliers and either deletes or adjusts them based on three defined thresholds.
##'
##' \tabular{ll}{
##' 'window' \tab ... \cr
##' 'outlier.mins' \tab ... \cr
##' 'stationary.mins' \tab ... \cr
##' }
##' @title Search and edit twilight outliers
##' @param twilights a dataframe with columns \code{Twilight} and \code{Rise}
##' @param offset the starting hour for the vertical axes.
##' @param window the number of neighboring twilights
##' @param outlier.mins threshold for outliers (in minutes)
##' @param stationary.mins threshold for variation of twiligths at stationary site (in mins)
##' @return the dataframe of edited twilights, with columns
##' \item{\code{Twilight}}{times of (edited) twilight}
##' \item{\code{Rise}}{logical indicating sunrise}
##' \item{\code{Deleted}}{logical indentifying deleted twilights}
##' \item{\code{Edited}}{logical indentifying edited twilights}
##' \item{\code{Twilight0}}{Original twilight time}
##' @importFrom zoo rollapply
##' @importFrom graphics plot
##' @export
twilightEdit <- function(twilights, offset = 17, window = 4, outlier.mins = 45, stationary.mins = 15, zlim=c(0,64), plot = T){
day <- twilights$Twilight
hour <- hourOffset(as.hour(twilights$Twilight),offset)
sunr <- which(twilights$Rise)
suns <- which(!twilights$Rise)
fnc <- function(x) {
ind0 <- abs(x[(window/2)+1,1] - x[-((window/2)+1),1]) > median(diff(as.numeric(day[suns])))*((window/2)+1)
if(any(ind0)) x <- x[-((window/2)+1),][-which(ind0),]
if(nrow(x)<window/2) {
out <- cbind(x[(window/2)+1,1], FALSE, FALSE, x[(window/2)+1,1])
} else {
diffr <- abs(x[(window/2)+1,2]-median(x[-((window/2)+1),2]))*60
if(diffr>=outlier.mins & all(dist(x[-((window/2)+1),2])*60<=stationary.mins)) {
out <- cbind(x[(window/2)+1,1] + (median(x[c(window/2,(window/2)+2),2]) - x[(window/2)+1,2])*60*60, FALSE, TRUE, x[(window/2)+1,1])
}
if(diffr>=outlier.mins & !all(dist(x[-((window/2)+1),2])*60<=stationary.mins)) {
out <- cbind(x[(window/2)+1,1], TRUE, FALSE, x[(window/2)+1,1])
}
if(diffr<outlier.mins) out <- cbind(x[(window/2)+1,1], FALSE, FALSE, x[(window/2)+1,1])
}
return(out)
}
sunrT <- rollapply(cbind(day,hour)[sunr,], width = window+1, FUN = fnc, fill = FALSE, by.column = F)
sunsT <- rollapply(cbind(day,hour)[suns,], width = window+1, FUN = fnc, fill = FALSE, by.column = F)
sunrT[which(sunrT[,1]==0), c(1,4)] <- day[sunr][which(sunrT[,1]==0)]
sunsT[which(sunsT[,1]==0), c(1,4)] <- day[suns][which(sunsT[,1]==0)]
out <- data.frame(Twilight = as.POSIXct(c(sunrT[,1], sunsT[,1]), origin = "1970-01-01", tz = "GMT"),
Rise = c(rep(TRUE, nrow(sunsT)), rep(FALSE, nrow(sunsT))),
Deleted = ifelse(c(sunrT[,2], sunsT[,2])==1, TRUE, FALSE),
Edited = ifelse(c(sunrT[,3], sunsT[,3])==1, TRUE, FALSE),
Twilight0 = as.POSIXct(c(sunrT[,4], sunsT[,4]), origin = "1970-01-01", tz = "GMT"))
out <- out[order(out[,1]),]
rownames(out) <- 1:nrow(out)
if(plot) {
day0 <- out$Twilight0
hour0 <- hourOffset(as.hour(out$Twilight0),offset)
day <- out$Twilight
hour <- hourOffset(as.hour(out$Twilight),offset)
plot(day0,hour0,type="n",xlab="Date",ylab="Hour")
points(day[!out$Deleted],hour[!out$Deleted],pch=16,cex=0.5,col=ifelse(out$Rise[!out$Deleted], "firebrick", "cornflowerblue"))
arrows(day0[out$Edited], hour0[out$Edited], day[out$Edited], hour[out$Edited], length = 0.1)
points(day0[out$Deleted | out$Edited],hour0[out$Deleted | out$Edited],pch=16,col="grey50")
points(day[out$Edited],hour[out$Edited],pch=16,col=ifelse(out$Rise[out$Edited], "firebrick", "cornflowerblue"))
points(day0[out$Deleted],hour0[out$Deleted], pch = "X")
}
out
}
##' Calibrate threshold light level to solar zenith angle
##'
##' Given data for a known location, this plots the recorded light
##' against the solar zenith angle to allow observed light levels to
##' be calibrated to solar zenith angle.
##' @title Light Threshold Calibration
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' and light levels recorded by the tag.
##' @param lon calibration longitude.
##' @param lat calibration latitude.
##' @param max.adjust adjust twilights for tags that report the
##' maximum light interval observed in the preceeding sampling
##' interval.
##' @param col colors for increasing and decreasing zenith angles
##' @param ... additional arguments to pass to \code{plot}.
##' @importFrom SGAT zenith solar
##' @importFrom graphics plot
##' @export
thresholdCalibrate <- function(tagdata,lon,lat,max.adjust=TRUE,
col=c("dodgerblue","firebrick"),...) {
Light <- tagdata$Light
Zenith <- zenith(solar(tagdata$Date),lon,lat)
if(max.adjust) {
n <- length(Light)
Light <- Light[-1]
Zenith <- pmin(Zenith[-1],Zenith[-n])
}
n <- length(Zenith)
dz <- Zenith[c(2:n,n)] - Zenith[c(1,1:(n-1))]
plot(Zenith,Light,col=ifelse(dz > 0,col[1],col[2]),...)
}
##' Adjust time for clock drift
##'
##' Linearly rescale a sequence of dates to a new start and end time
##' to correct for clock drift in the tag.
##' @title Clock Drift Adjustment
##' @param time a vector of POSIXct times.
##' @param start new start time as POSIXct.
##' @param end new end time as POSIXct.
##' @importFrom stats approx
##' @export
driftAdjust <- function(time,start,end) {
n <- length(time)
.POSIXct(approx(time[c(1,n)],c(start,end),time)$y)
}
##' Summarise auxiliary data in specified intervals
##'
##' Given a time series of responses and a sequence of start and end
##' times defining time intervals spanned by the time series,
##' \code{intervalData} computes a summary statistic for each
##' interval. The \code{twilightData} is similar, but the time
##' intervals are defined as a number of hours before and after a
##' given sequence of twilight times.
##'
##' @title Extract auxiliary tag data
##' @param date the sequence of sample times as POSIXct.
##' @param y the sequence of responses.
##' @param start the start of the time interval to summarize.
##' @param end the end of the time interval to summarize.
##' @param twl the sequence of estimated twilight times as POSIXct.
##' @param before the interval (hours) before twilight to include in the interval.
##' @param after the interval (hours) after twilight to include in the interval.
##' @param FUN summary function
##' @return the data summary for each twilight.
##' @export
intervalData <- function(date,y,start,end,
FUN=function(y) if(all(is.na(y))) NA else mean(y,na.rm=TRUE)) {
mapply(
function(start,end) FUN(y[start <= date & date <= end]),
start,end)
}
##' @rdname intervalData
##' @export
twilightData <- function(date,y,twl,before=4,after=4,
FUN=function(y) if(all(is.na(y))) NA else mean(y,na.rm=TRUE)) {
mapply(
function(start,end) FUN(y[start <= date & date <= end]),
twl-3600*before,twl+3600*after)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.