R/BAStag.R

In SWotherspoon/BAStag: Basic data processing for light based geolocation archival tags

#' BAS tag data processing
#'
#' Basic data processing for British Antarctic Survey 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 BAStag-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
}


##' Read Temperature data from a Migrate Technologies tag
##'
##' The \code{readMTdeg} function reads temperature data from a
##' Migrate Technologies tag.
##'
##' @title Read Temperature Data
##' @param file the sst file to import.
##' @param skip number of initial lines to skip
##' @return Return 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{Wet}}{the number of wet samples}
##' \item{\code{Conductivity}}{the condicivity}
##' @importFrom utils read.csv
##' @export
readMTdeg <- 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","MinSST","MaxSST","SST","Wet","Conductivity"),
                  colClasses=c("character","numeric","numeric","numeric","numeric","numeric"))
  ## Parse date
  d$Date <- as.POSIXct(strptime(d$Date,"%d/%m/%Y %H:%M:%S",tz="GMT"))
  ## Return data frame.
  d
}



##' 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 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)))
}





##' 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 BAS light recording.
##'
##' BAS 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
}



##' 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)
}