R/interactive.R

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

defaultPalette <- c(red="#E41A1C", blue="#377EB8", green="#4DAF4A", violet="#984EA3",
                    orange="#FF7F00", yellow="#FFFF33", brown="#A65628", pink="#F781BF",
                    black="#000000", grey3="#333333", grey6="#666666", grey9="#999999", greyC="#CCCCCC")




##' Mouse button identification
##'
##' Identify left and right mouse button clicks from the buttons
##' vector passed to the callbacks used by getGraphicsEvent.  This is
##' non-trivial because the windows and MacOS behaviours are subtly
##' different.
##'
##' @title Mouse button identification
##' @param buttons vector passed to callbacks used by getGraphicsEvent.
##' @return Returns 1 for left click, 2 for right click and 0
##' otherwise.
mouseButton <- function(buttons) {
  n <- length(buttons)
  ## This is right click on a mac
  if(n >= 2 && buttons[1]==0 && buttons[2]==1) return(2)

  ## This is both buttons on windows => treat as neither.
  if(n >= 2 && buttons[1]==0 && buttons[2]==2) return(0)

  ## This is right click on windows
  if(n >= 1 && buttons[1]==2) return(2)

  ## This is left click
  if(n >= 1 && buttons[1]==0) return(1)

  ## Other combinations => 0
  0
}




##' Image plot with twilight points and ribbons
##'
##' Draws a light image with twilight points or twilight ribbon. If
##' tagdata is \code{NULL}, no image is drawn.  Twilight points and
##' ribbons are only drawn when twilight data is provided and the
##' point or ribbon colour specifications are non \code{NULL}.
##'
##' @title Image plot with twilights
##' @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 dataframe of twilight times as generated by
##' \code{\link{findTwilights}}.
##' @param offset the starting hour for the vertical axes.
##' @param xlim the range of dates to plot.
##' @param mark twilights to mark with a cross
##' @param points.col colours of twilight points
##' @param ribbon.col colours of twilight ribbons
##' @param zlim the range of light levels to plot.
##' @param point.cex expansion factor for plot points.
##' @importFrom graphics plot points
imageDraw <- function(tagdata=NULL,twilights=NULL,offset=0,xlim=NULL,
                       mark=NULL,points.col=NULL,ribbon.col=NULL,
                       zlim=c(0,64),point.cex=0.6) {

  ## Plot background image
  if(!is.null(tagdata)) {
    if(!is.null(xlim))
      lightImage(tagdata,offset=offset,zlim=zlim,xlim=xlim)
    else
      lightImage(tagdata,offset=offset,zlim=zlim)
  }

  ## If there is twilight data
  if(!is.null(twilights)) {

    day <- twilights$Twilight
    hour <- hourOffset(as.hour(twilights$Twilight),offset)

    ## Initialize plot axes
    if(is.null(tagdata)) plot(day,hour,type="n",xlab="Date",ylab="Hour",ylim=c(offset,offset+24))

    ## Plot twilight times
    if(length(points.col)>0) points(day,hour,pch=16,cex=point.cex,col=points.col)

    ## Plot twilights segments
    if(length(ribbon.col)>0) {
      ribbon.col <- rep(ribbon.col,length.out=2)
      rise <- twilights[twilights$Rise,]
      tsimageRibbon(.POSIXct(tapply(rise$Start,rise$Twilight,min),"GMT"),
                     .POSIXct(tapply(rise$End,rise$Twilight,max),"GMT"),
                     offset=offset,border=NA,col=ribbon.col[1])
      set <- twilights[!twilights$Rise,]
      tsimageRibbon(.POSIXct(tapply(set$Start,set$Twilight,min),"GMT"),
                     .POSIXct(tapply(set$End,set$Twilight,max),"GMT"),
                     offset=offset,border=NA,col=ribbon.col[2])
    }
    ## Mark points
    if(length(mark)>0) points(day[mark],hour[mark],pch=3)
  }
}

##' Initialize a light profile plot for plotting
##'
##' Sets up axes for drawing light profile plots.
##'
##' @title Light profile plot
##' @param date a list of three sequences of sample times as POSIXct.
##' @param light a list of three sequences of observed light levels
##' @param lmax the maximum light level.
##' @param xlab the x axis label.
##' @param main the main title.
##' @param threshold threshold levels to display
##' @param lag show profiles from neighbouring days
##' @param point show individual observations as points?
##' @param point.cex expansion factor for plot points.
##' @param profile.col the colours of the three light profiles.
##' @param threshold.col the colour of the threshold markers.
##' @importFrom graphics axis.POSIXct plot plot.new plot.window
profileInit <- function(date,light,lmax=64,xlab="",main="") {
  if(length(date[[2]]) > 0) {
    ## Draw axes for light profiles
    plot(date[[2]],light[[2]],ylim=c(0,lmax),xlab=xlab,ylab="Light",type="n",xaxt="n",main=main)
    axis.POSIXct(1,x=date[[2]],format="%H:%M")
  } else {
    plot.new()
  }
}

##' @rdname profileInit
##' @importFrom graphics abline lines points
profileOverlay <- function(date,light,threshold=NULL,point=FALSE,lag=TRUE,
                            profile.col=defaultPalette[c(9,5,2)],
                            threshold.col=defaultPalette[1],
                            point.cex=0.6) {

  ## Overlay with light threshold
  if(!is.null(threshold)) abline(h=threshold,col=threshold.col)
  ## Overlay the light profile for the current and surrounding days
  if(lag) {
    lines(date[[1]]+86400,light[[1]],col=profile.col[2])
    lines(date[[3]]-86400,light[[3]],col=profile.col[3])
  }
  lines(date[[2]],light[[2]],col=profile.col[1])
  ## Overlay observations.
  if(point) points(date[[2]],light[[2]],col=profile.col[1],pch=16,cex=point.cex)
}

##' Draw rectangle above plot
##'
##' Draws a sequence of rectangles above a plot to indicate sections
##' of the data that have been selected.  Drawing the rectangle above
##' the plot means the selection can be updated without having to
##' redraw the main plot.
##'
##' @title Selection rectangle
##' @param x1 coordinates of the left boundaries of rectangles
##' @param x2 coordinates of the right boundaries of rectangles
##' @param col rectangle colours
##' @param add add rectangles to existing rectangles.
##' @importFrom graphics grconvertX grconvertY rect
selectionRectangle <- function(x1,x2,col,add=FALSE) {
  ## Determine upper and lower limits
  rx <- grconvertX(c(0,1),from="npc",to="user")
  ry <- grconvertY(c(1.01,1.03),from="npc",to="user")
  ## Over plot with white
  if(!add) rect(rx[1],ry[1],rx[2],ry[2],border=NA,col="white",xpd=NA)
  if(length(x1) >0 && length(x2)>0)
    rect(x1,ry[1],x2,ry[2],border=NA,col=col,xpd=NA)
}



##' Convert ndc coordinates
##'
##' Given user coordinates of a sequence of points and the ndc
##' coordinates (x,y), \code{ndcClosest} returns the index of the
##' closest point to (x,y).  Given the ndc coordinates of a pixel in a
##' tsimage plot, ndcTsimageDate returns the date of the pixel.
##'
##' @title Graphics coordinate conversion
##' @param x ndc x coordinate
##' @param y ndc y coordinate
##' @param xs points user x coordinates
##' @param ys points user y coordinates
##' @return index of closest point
##' @importFrom graphics grconvertX grconvertY
ndcClosest <- function(x,y,xs,ys) {
  xs <- grconvertX(xs,from="user",to="ndc")
  ys <- grconvertY(ys,from="user",to="ndc")
  which.min((x-xs)^2+(y-ys)^2)
}

##' @importFrom graphics grconvertX grconvertY
##' @rdname ndcClosest
ndcTsimageDate <- function(x,y) {
  day <- .POSIXct(grconvertX(x,from="ndc",to="user"),"GMT")
  hour <- grconvertY(y,from="ndc",to="user")
  .POSIXct(day+(hour%%24-as.hour(day))*60*60,"GMT")
}


##' Interactively delete observations from a time series
##'
##' Interactively delete observations from a time series.  The time
##' series is displayed. Left mouse button clicks zoom on a region of
##' the plot; right mouse button clicks toggle deletion of individual
##' points.
##'
##' \tabular{ll}{
##' 'q' \tab Quits, returning a indicator vector \cr
##' 'r' \tab Reset zoom to full plot range\cr
##' 's' \tab Toggle reseting of y axes scale on zoom\cr
##' '+'/'-' \tab Adjust the level of zoom\cr
##' 'Left/Right arrow' \tab Jump forward or backward in hte time series \cr
##' }
##'
##'
##' @title Select data
##' @param date sample times.
##' @param r recorded response.
##' @param deleted initial logical vector indicating which
##' observations are to be deleted.
##' @param extend the period (in hours) before and after the selection
##' to be shown in the zoom window.
##' @param xlab label for the x axis
##' @param ylab label for the y axis
##' @param width width of the interface windows.
##' @param height height of the interface windows.
##' @param palette a colour palette of 4 colours.
##' @param ... additional parameters to pass to plot.
##' @return a logical vector that indicates the observations to be deleted.
##' @importFrom graphics grconvertX grconvertY plot
##' @importFrom grDevices dev.cur dev.new dev.off dev.set getGraphicsEvent setGraphicsEventHandlers
##' @export
selectData <- function(date,r,deleted=NULL,extend=48,
                        xlab="Date",ylab="",width=12,height=4,
                        palette=defaultPalette[c(2,1)],...) {

  if(is.null(deleted)) deleted <- logical(length(date))
  zoom <- 3600*extend
  minDate <- min(date)
  maxDate <- max(date)
  xlim <- c(minDate,maxDate)
  ylim <- range(r,na.rm=TRUE)
  rescale <- FALSE


  ## Select device
  setDevice <- function(dev) if(dev.cur()!=dev) dev.set(dev)
  ## Focus if possible
  focus <- if(exists("bringToTop",mode="function")) grDevices::bringToTop else setDevice

  ## Draw the selection window
  winADraw <- function() {
    setDevice(winA)
    ## Plot data
    keep <- date >= xlim[1] & date <= xlim[2]
    if(rescale) ylim <<- range(r[keep],na.rm=TRUE)
    plot(date[keep],r[keep],xlab=xlab,ylab=ylab,ylim=ylim,
         col=ifelse(deleted[keep],palette[2],palette[1]),...)
  }

  ## onMouseDown callback for selection window.
  winAOnMouseDown <- function(buttons,x,y) {
    setDevice(winA)

    if(length(buttons) > 0) {
      b <- mouseButton(buttons)
      ## Add segments to include or exclude
      if(b==1) {
        ## Set zoom
        zoom <<- 3600*extend
        x <- grconvertX(x,from="ndc",to="user")
        xlim <<- .POSIXct(x+zoom*c(-1,1),"GMT")
      }
      if(b==2) {
        ## Toggle deletion of closest point
        k <- ndcClosest(x,y,date,r)
        deleted[k] <<- !deleted[k]
      }
      winADraw()
    }
    NULL
  }

  ## onKeybd callback for both windows
  onKeybd <- function(key) {
    ## q quits
    if(key=="q") return(-1)

    ## r : reset plot
    if(key=="r") {
      xlim <<- c(minDate,maxDate)
      zoom <<- 3600*extend
    }

    if(key=="s") {
      rescale <<- !rescale
      if(!rescale) ylim <<- range(r,na.rm=TRUE)
    }

    ## +/- : expand/shrink zoom
    if(key=="+") {
      zoom <<- zoom/2
      xlim <<- .POSIXct(mean(xlim)+zoom*c(-1,1),"GMT")
    }

    if(key=="-") {
      zoom <<- 2*zoom
      if(zoom < as.numeric(maxDate)-as.numeric(minDate)) {
        xlim <<- .POSIXct(mean(xlim)+zoom*c(-1,1),"GMT")
      } else {
        xlim <<- c(minDate,maxDate)
        zoom <<- 3600*extend
      }
    }

    ## Left/Right : jump to neighbouring chunk
    if(key=="Left") {
      xlim <<- .POSIXct(xlim-2*zoom,"GMT")
      if(xlim[1] < minDate) xlim <<- .POSIXct(c(minDate,minDate+2*zoom),"GMT")
    }
    if(key=="Right") {
      xlim <<- .POSIXct(xlim+2*zoom,"GMT")
      if(xlim[2] > maxDate) xlim <<- .POSIXct(c(maxDate-2*zoom,maxDate),"GMT")
    }
    ## Redraw
    winADraw()
    NULL
  }

  ## Set up master window
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winA <- dev.cur()
  winADraw()
  focus(winA)
  setGraphicsEventHandlers(
    which=winA,
    prompt="Select Data",
    onMouseDown=winAOnMouseDown,
    onKeybd=onKeybd)
  tryCatch({
    getGraphicsEvent()
    dev.off(winA)
    deleted
  }, finally=deleted)
}

##' Interactively edit crepuscular intervals
##'
##' Interactively edit the crepuscular intervals based on the light
##' profile.  A plot of the estimated sunrise and sunset intervals is
##' displayed, and the user can select the twilight to be edited with
##' a left mouse click.
##'
##' The light profile for the selected twilight is shown in a separate
##' window, and the selected segments of the light profile are
##' highlighted.  The corresponding light profiles from the preceeding
##' and following days are also shown for reference.
##'
##' The user may select a new candidate interval by clicking and
##' dragging with the left mouse button.  Individual points may be
##' selected or deselected with the right mouse button.
##'
##' Note however, no actual change to the selection is made until the
##' candidate edits are accepted depressing the 'a' key.
##'
##' Each twilight may be marked with the an integer 0 to 9 with the
##' numeric keys. By default each day is given the mark 0.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits, returning the dataframe of edited twilight segments \cr
##' 'a' \tab Accepts the candidate edit \cr
##' 'r' \tab Resets the selection \cr
##' 'i' \tab Toggles the display of the light image \cr
##' 'p' \tab Toggles the display of individual points \cr
##' 'l' \tab Toggles the display of surrounding profiles \cr
##' '+'/'-' \tab Zoom in or out \cr
##' 'Left arrow' \tab Jump to previous twilight \cr
##' 'Right arrow' \tab Jump to next twilight \cr
##' '0'-'9' \tab Mark this twilight \cr
##' }
##'
##' @title Edit crepuscular 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 dataframe of twilight times as generated by
##' \code{\link{findCrepuscular}}.
##' @param offset the starting hour for the vertical axes.
##' @param extend the period (in hours) before and after twilight for
##' which the light profile should be plotted.
##' @param threshold the light threshold that defines twilight.
##' @param lmax the maximum light level to plot.
##' @param zlim the range of light levels to plot in images.
##' @param point.cex expansion factor for plot points.
##' @param width width of the interface windows.
##' @param height height of the interface windows.
##' @param palette a colour palette of 8 colours.
##' @seealso \code{\link{preprocessLight}}
##' @return the dataframe of edited twilights, with columns
##' \item{\code{Twilight}}{edited times of twilight}
##' \item{\code{Rise}}{logical indicating sunrise}
##' \item{\code{Start}}{date of first observation in the crepuscular segment}
##' \item{\code{End}}{date of last observation in the crepuscular segment}
##' @importFrom graphics grconvertX grconvertY lines points
##' @importFrom grDevices dev.cur dev.new dev.off dev.set getGraphicsEvent setGraphicsEventHandlers
##' @export
crepuscularEdit <- function(tagdata,twilights,offset=0,extend=6,threshold=NULL,lmax=64,zlim=c(0,lmax),
                             point.cex=0.5,width=12,height=4,
                             palette=defaultPalette[c(5,2,9,3,4,1,1)]) {

  ## Order twilights
  if(is.null(twilights$Marker)) twilights$Marker <- integer(nrow(twilights))
  twilights <- twilights[order(twilights$Twilight,twilights$Start),
                         c("Twilight","Rise","Start","End","Marker")]
  ## Extract date and hour of twilight
  day <- twilights$Twilight
  hour <- hourOffset(as.hour(twilights$Twilight),offset)

  ## Cached data subsets
  index <- 1
  first <- 1
  indices <- as.numeric(factor(as.numeric(twilights$Twilight)))
  date <- vector(3,mode="list")
  lght <- vector(3,mode="list")
  selected <- NULL
  start <- end <- 0
  changed <- FALSE
  showobs <- FALSE
  showimg <- FALSE
  showlag <- TRUE

  ## Set cached values
  cache <- function(k) {
    index <<- k
    changed <<- FALSE
    first <<- which(index==indices)[1]
    ## Get profiles
    for(k in 1:3) {
      mid <- twilights$Twilight[index]+(k-2)*86400
      keep <- (tagdata$Date >= mid-3600*extend) & (tagdata$Date <= mid+3600*extend)
      date[[k]] <<- tagdata$Date[keep]
      lght[[k]] <<- tagdata$Light[keep]
    }
    ## Determined selected range
    selected <<- logical(length(date[[2]]))
    for(i in which(index==indices))
      selected[date[[2]] >= twilights$Start[i] & date[[2]] <= twilights$End[i]] <<- TRUE
  }

  ## Select device
  setDevice <- function(dev) if(dev.cur()!=dev) dev.set(dev)
  ## Focus if possible
  focus <- if(exists("bringToTop",mode="function")) grDevices::bringToTop else setDevice

  ## Draw the twilights window
  winADraw <- function() {
    setDevice(winA)
    imageDraw(if(showimg) tagdata,
               twilights,offset=offset,mark=index,
               ribbon.col=palette[1:2],
               zlim=zlim,point.cex=point.cex)
  }


  ## Draw axes for light profiles
  winBInit <- function() {
    setDevice(winB)
    marker <- twilights$Marker[first]
    profileInit(date,lght,lmax=lmax,
                 xlab=if(marker>0) paste("Marker: ",marker) else "",
                 main=as.character(twilights$Twilight[first]))
  }

  ## Draw light profiles
  winBDraw <- function() {
    setDevice(winB)
    ## Overlay with light profiles
    profileOverlay(date,lght,threshold,showobs,showlag,
                    profile.col=palette[3:5],threshold.col=palette[6],
                    point.cex=point.cex)

    ## Show selection
    col <- palette[if(changed) 7 else (if(twilights$Rise[first]) 1 else 2)]
    ## Hightlight selected segments
    x <- ifelse(selected,date[[2]],NA)
    y <- ifelse(selected,lght[[2]],NA)
    lines(x,y,col=col)
    if(showobs) points(x,y,pch=16,cex=point.cex,col=col)

    ## Selection rectangle
    x1 <- x2 <- NULL
    if(any(selected)) {
      x1 <- date[[2]][diff(c(FALSE,selected))==1]
      x2 <- date[[2]][diff(c(selected,FALSE))==-1]
    }
    selectionRectangle(x1,x2,col)
  }


  ## onMouseDown callback for twilights window.
  winAOnMouseDown <- function(buttons,x,y) {
    setDevice(winA)
    if(length(buttons) > 0 && mouseButton(buttons)==1) {
      ## Determine selected profile.
      setDevice(winA)
      xs <- grconvertX(c(day,day,day),from="user",to="ndc")
      ys <- grconvertY(c(hour-24,hour,hour+24),from="user",to="ndc")
      k <- (which.min((x-xs)^2+(y-ys)^2)-1)%%length(day)+1
      ## Redraw
      cache(k)
      winADraw()
      winBInit()
      winBDraw()
      focus(winB)
    }
    NULL
  }

  ## onKeybd callback for both windows
  onKeybd <- function(key) {
    ## q quits
    if(key=="q") return(-1)
    ## +/- : zoom time window around threshold crossing
    if(key=="+") {
      extend <<- max(extend-1,1)
      cache(index)
    }
    if(key=="-") {
      extend <<- min(extend+1,24)
      cache(index)
    }
    ## x : reset selection
    if(key=="r") {
      cache(index)
    }
    ## i : toggle display light image
    if(key=="i") {
      showimg <<- !showimg
    }
    ## p : toggle display of points in the profile window
    if(key=="p") {
      showobs <<- !showobs
    }
    ## l : toggle display of lagged twilights in the profile window
    if(key=="l") {
      showlag <<- !showlag
    }
    ## Left/Right : jump to neighbouring twilight
    if(key=="Left") {
      cache(max(index-1,1))
    }
    if(key=="Right") {
      cache(min(index+1,max(indices)))
    }
    ## a : accept current edit
    if(key=="a") {
      d <- split(twilights,twilights$Twilight)
      d[[index]] <- cbind(Twilight=d[[index]]$Twilight[1],
                          Rise=d[[index]]$Rise[1],
                          data.frame(Start=date[[2]][diff(c(FALSE,selected))==1],
                                     End=date[[2]][diff(c(selected,FALSE))==-1]),
                          Marker=d[[index]]$Marker[1])
      d <- do.call("rbind",d)
      d$Twilight <- .POSIXct(d$Twilight,"GMT")
      d$Start <- .POSIXct(d$Start,"GMT")
      d$End <- .POSIXct(d$End,"GMT")
      twilights <<- d
      indices <<- as.numeric(factor(as.numeric(twilights$Twilight)))
      changed <<- FALSE
    }
    if(key >= "0" && key <= "9") {
      twilights$Marker[which(index==indices)] <<- as.numeric(key)
    }

    ## Redraw
    winADraw()
    winBInit()
    winBDraw()
    NULL
  }

  ## onMouseDown callback for profile window
  winBOnMouseDown <- function(buttons,x,y) {
    if(length(buttons) > 0) {
      b <- mouseButton(buttons)
      ## Button 1 -> record location and do complete draw
      if(b==1) {
        changed <<- TRUE
        start <<- grconvertX(x,from="ndc",to="user")
        selected <<- logical(length(date[[2]]))
        end <<- start
        winBInit()
      }
      ## Button 2 -> toggle selected points
      if(b==2) {
        changed <<- TRUE
        xs <- grconvertX(date[[2]],from="user",to="ndc")
        ys <- grconvertY(lght[[2]],from="user",to="ndc")
        k <- which.min((x-xs)^2+(y-ys)^2)
        selected[k] <<- (selected[k]==FALSE)
      }
    }
    winBDraw()
    NULL
  }

  ## onMouseMove callback for profile window
  winBOnMouseMove <- function(buttons,x,y) {
    ## Button 1 drag to select crepuscular period
    if(length(buttons) > 0 && mouseButton(buttons)==1) {
      end <<- grconvertX(x,from="ndc",to="user")
      sel <- (date[[2]] >= min(start,end) & date[[2]] <= max(start,end))
      if(any(sel!=selected)) {
        selected <<- sel
        winBDraw()
      }
    }
    NULL
  }

  ## Set up twilights window
  index <- 1
  cache(index)
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winA <- dev.cur()
  winADraw()
  setGraphicsEventHandlers(
    which=winA,
    prompt="Select Twilight",
    onMouseDown=winAOnMouseDown,
    onKeybd=onKeybd)
  ## Set up profile window
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winB <- dev.cur()
  winBInit()
  winBDraw()
  setGraphicsEventHandlers(
    which=winB,
    prompt="Light Profile",
    onMouseDown=winBOnMouseDown,
    onMouseMove=winBOnMouseMove,
    onKeybd=onKeybd)
  focus(winA)
  ## Monitor for events
  tryCatch({
      getGraphicsEvent()
      dev.off(winB)
      dev.off(winA)
      twilights
  }, finally=twilights)
}






##' Overlay a data when editing an initial path.
##'
##' The \code{pathEdit} function allows the user to specify a
##' function to overlay data on the map to assist is selecting an
##' initial path.
##'
##' This function overlays a contour map of the difference between the
##' observed twilight and true time of twilight for that location and
##' day, where positive and negative contours are shown in different
##' colours.  This can be useful is selecting an initial path for
##' which the residuals are positive.
##'
##' @title Overlay data for path editing
##' @param twilights the full twilight data
##' @param index the index of the current twilight
##' @param mode the overlay mode requested by \code{pathEdit}
##' @param xlim the horizontal limits of the map
##' @param ylim the vertical limits of the map
##' @param twilight.contours the contour levels to display
##' @param zenith the solar zenith angle that defines twilight
##' @seealso \code{\link{pathEdit}}
##' @importFrom raster raster
##' @importFrom SGAT twilightResidualsMap
##' @importFrom graphics contour
##' @importFrom grDevices col2rgb rgb
##' @export
overlayTwilightResiduals <- function(twilights,index,mode,xlim,ylim,
                                       twilight.contours=c(10,20,50),
                                       zenith=NULL) {
  addAlpha <- function(col,alpha) {
    col <- col2rgb(col)/255
    rgb(red=col[1],green=col[2],blue=col[3],alpha=alpha)
  }

  if(mode!=0 && !is.null(zenith)) {
    grid <- raster(nrows=30,ncols=30,xmn=xlim[1],xmx=xlim[2],ymn=ylim[1],ymx=ylim[2])
    grid <- twilightResidualsMap(twilights$Twilight[index],twilights$Rise[index],grid,zenith=zenith)
    contour(grid,add=TRUE,levels=c(0,twilight.contours),col=addAlpha(defaultPalette[3],0.5))
    contour(grid,add=TRUE,levels=-twilight.contours,col=addAlpha(defaultPalette[4],0.5))
  }
}


##' Interactively edit a path of twilight locations.
##'
##' Interactively edit a path of twilight locations.  A plot of the
##' estimated sunrise and sunset times is displayed, and the user can
##' select the location corresponding to a particular twilight with a
##' left mouse click.
##'
##' The path is dislayed in another window, with the editable location
##' highlighted.  The user can select the point to edit with a left
##' mouse click, or move the editable location with a right mouse
##' click.  In auto advance mode, when a location is edited, the
##' editable location advances to the next location in the sequence.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits, returning the dataframe of edited twilight segments \cr
##' 'a' \tab Toggle auto advance mode \cr
##' 'r' \tab Resets the zoom to the encompass the entire track \cr
##' 'c' \tab Centres the zoomed window on the current point \cr
##' 'z' \tab Zooms to the locations surrounding the current location \cr
##' 'u' \tab Resets any edits to the current point \cr
##' '+'/'-' \tab Zoom in or out \cr
##' 'Left arrow' \tab Jump to previous location \cr
##' 'Right arrow' \tab Jump to next location \cr
##' '0'-'9' \tab Set the mode of the overlay plot \cr
##' }
##'
##' The user may supply a function \code{plotMap} that plots a
##' background map.  This must be a function of two arguments
##' \code{xlim} and \code{ylim} the determine the extent of the plot.
##'
##' The user may also supply a function \code{plot.overlay} that
##' overlays twilight specific data on the map.  This must be a
##' function of six arguments, \code{mode} is an integer flag used to
##' indicate the style of the overlay plot, \code{twilights} is the
##' dataframe of twilight data, \code{index} is the index of the
##' current twilight, \code{xlim} and \code{ylim} determine the extent
##' of the plot, and \code{...} is used to pass any additional
##' arguments.  The mode of the overlay plot is set with the keys 0 to
##' 9. We adopt the convention that mode 0 suppresses the overlay, but
##' this is not enforced.
##'
##' The user may also supply a function \code{is.invalid} that accepts
##' the current path as an argument and returns a logical vector
##' indicating which locations along the path are in some way invalid.
##' The twilights for these locations will be highlighted.
##'
##' @title  Adjust a path
##' @param path a two column matrix of the (lon,lat) locations at the
##' twilight times.
##' @param twilights dataframe of twilight times as generated by
##' \code{\link{findTwilights}}.
##' @param offset the starting hour for the vertical axes.
##' @param fixed logical vector indicating which locations to hold
##' fixed.
##' @param aspect aspect ratio of the map.
##' @param extend the number of locations before and after the current
##' location to highlight.
##' @param auto.advance advance to next point afet edit.
##' @param plotMap A function to plot the background map.
##' @param plot.overlay A function to overlay twilight specific data on the map.
##' @param is.invalid A function that indicates if a location is not valid.
##' @param point.cex expansion factor for plot points.
##' @param width width of the selection window.
##' @param height height of the selection window.
##' @param map.width width of the map window.
##' @param map.height height of the map window.
##' @param palette a colour palette of 6 colours.
##' @param ... additional arguments passed to plot.overlay
##' @return a two column matrix of (lon,lat) locations.
##' @importFrom graphics grconvertX grconvertY lines par plot.new plot.window points
##' @importFrom grDevices dev.cur dev.new dev.off dev.set getGraphicsEvent setGraphicsEventHandlers recordPlot replayPlot
##' @export
pathEdit <- function(path,twilights,offset=0,fixed=FALSE,
                      aspect=1,extend=3,auto.advance=FALSE,
                      plotMap=NULL,plot.overlay=NULL,
                      is.invalid=function(path) logical(nrow(path)),
                      point.cex=0.5,width=12,height=4,
                      map.width=8,map.height=8,
                      palette=defaultPalette[c(5,2,1,12,12,1)],
                      ...) {


  ## Store original path
  path0 <- path
  ## Order twilights and check deleted
  twilights <- twilights[order(twilights$Twilight),]
  ## Extract date and hour of twilight
  day <- twilights$Twilight
  hour <- hourOffset(as.hour(twilights$Twilight),offset)
  ## Set fixed points
  fixed <- rep(fixed,length.out=nrow(path))
  invalid <- !fixed & is.invalid(path)
  ## Overlay mode
  mode <- 0


  ## Map scale parameters
  xlim <- ylim <- centre <- xyscl <- NULL
  setZoom <- function(w) {
    xyscl <<- w
    xlim <<- centre[1]+xyscl*c(-0.5,0.5)
    ylim <<- pmax(pmin(centre[2]+aspect*xyscl*c(-0.5,0.5),90),-90)
    map <<- NULL
  }

  setWindow <- function(path) {
    xl <- range(path[,1])
    yl <- range(path[,2])
    centre <<- c(mean(xl),mean(yl))
    setZoom(max(diff(xl),diff(yl)/aspect))
  }
  setWindow(path)
  ## Cached map
  map <- NULL

  if(is.null(plotMap))
    plotMap <- function(xlim,ylim) { plot.new(); plot.window(xlim,ylim) }

  ## Select device
  setDevice <- function(dev) if(dev.cur()!=dev) dev.set(dev)
  ## Focus if possible
  focus <- if(exists("bringToTop",mode="function")) grDevices::bringToTop else setDevice

  ## Draw the twilights window
  winADraw <- function() {
    setDevice(winA)
    col <- palette[ifelse(fixed,4,
                          ifelse(is.na(invalid) | invalid,3,
                                 ifelse(twilights$Rise,1,2)))]
    cex <- point.cex*ifelse(is.na(invalid) | invalid,2,1)
    imageDraw(NULL,twilights,offset=offset,mark=index,
               points.col=col,point.cex=cex)
  }

  ## Draw path
  winBDraw <- function() {
    setDevice(winB)
    ## Create underlying map
    if(is.null(map)) {
      ## User defined map function
      plotMap(xlim,ylim)
      map <<- recordPlot()
      p <- par()$usr
      xlim <<- p[1:2]
      ylim <<- p[3:4]
    } else {
      ## Replot stored plot
      replayPlot(map)
    }

    ## Overlays
    if(!is.null(plot.overlay))
      plot.overlay(twilights,index,mode,xlim,ylim,...)
    ## Show full path
    lines(path[,1],path[,2],col=palette[5])
    points(path[,1],path[,2],col=palette[5],pch=16,cex=0.4)
    ## Highlight current point
    ks <- max(1,index-extend):min(nrow(path),index+extend)
    lines(path[ks,1],path[ks,2],col=palette[6])
    points(path[index,1],path[index,2],col=palette[if(fixed[index]) 4 else 6],pch=16,cex=1)
  }

  ## onMouseDown callback for twilights window.
  winAOnMouseDown <- function(buttons,x,y) {
    setDevice(winA)
    ## Determine selected profile.
    index <<- (ndcClosest(x,y,c(day,day,day),c(hour-24,hour,hour+24))-1)%%length(day)+1
    ## Redraw
    winADraw()
    winBDraw()
    focus(winB)
    NULL
  }

  ## onKeybd callback for both windows
  onKeybd <- function(key) {
    ## q quits
    if(key=="q") return(-1)
    if(key=="a") {
      auto.advance <<- !auto.advance
    }
    ## Reset zoom to the full path
    if(key=="r") {
      setWindow(path)
    }
    ## Centre on current point
    if(key=="c") {
      centre <<- path[index,]
      setZoom(xyscl)
    }

    ## Zoom to locations surrounding current location
    if(key=="z") {
      setWindow(path[max(1,index-extend):min(nrow(path),index+extend),])
    }
    ## +/- : zoom map window around current centre
    if(key=="+") {
      setZoom(5/8*xyscl)
    }
    if(key=="-") {
      setZoom(8/5*xyscl)
    }
    ## Left/Right : jump to neighbouring twilight
    if(key=="Left") {
      index <<- max(index-1,1)
    }
    if(key=="Right") {
      index <<- min(index+1,nrow(path))
    }
    ## Undo edit - reset to original location
    if(key=="u") {
      path[index,] <<- path0[index,]
      invalid <<- !fixed & is.invalid(path)
    }
    if(key >= "0" && key <= "9") {
      mode <<- key
    }

    ## Redraw
    winADraw()
    winBDraw()
    NULL
  }

  ## onMouseDown callback for path window
  winBOnMouseDown <- function(buttons,x,y) {
    setDevice(winB)
    if(length(buttons) > 0) {
      b <- mouseButton(buttons)
      ## Button 1 -> select point
      if(b==1) {
        ## Select nearest point
        index <<- ndcClosest(x,y,path[,1],path[,2])
      }
      ## Button 2 -> move location
      if(b==2) {
        if(!fixed[index]) {
          path[index,] <<- c(grconvertX(x,from="ndc",to="user"),
                             grconvertY(y,from="ndc",to="user"))
          invalid <<- !fixed & is.invalid(path)
        }
        if(auto.advance) index <<- min(index+1,nrow(path))
      }
    }
    winADraw()
    winBDraw()
    NULL
  }

  ## Set up twilights window
  index <- 1
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winA <- dev.cur()
  winADraw()
  setGraphicsEventHandlers(
    which=winA,
    prompt="Select Twilight",
    onMouseDown=winAOnMouseDown,
    onKeybd=onKeybd)
  ## Set up path window
  dev.new(width=map.width,height=map.height,noRStudioGD=TRUE)
  winB <- dev.cur()
  winBDraw()
  setGraphicsEventHandlers(
    which=winB,
    prompt="Path",
    onMouseDown=winBOnMouseDown,
    onKeybd=onKeybd)
  focus(winA)
  ## Monitor for events
  tryCatch({
    getGraphicsEvent()
    dev.off(winB)
    dev.off(winA)
    path
  }, finally=path)
}




##' Interactively derive and edit twilight times.
##'
##' This function allows the user to interactively search for and edit
##' twilight times corresponding to a given light threshold. The
##' process consists of four stages,
##' \enumerate{
##' \item Subset -- selection of a subset of the data for processing
##' \item Search -- semi-automated search for the twilights
##' \item Insert -- optionally, twilights are inserted where the light record is incomplete
##' \item Edit -- optioanlly, individual twilights are manually adjusted based on the light profiles.
##' }
##'
##' At each stage, user is presented with two windows.  The first
##' shows the data in its entirety and allows the user to select a
##' subset of the data to be manipulated in the second window.
##'
##' In the first stage the user can restrict the data to a subset for
##' processing.  The first window shows a light image and a coloured
##' rectangle above the image shows the selected subset of data.  In
##' the first window, clicking with either mouse button produces a
##' zoomed image of the clicked region in the second window.  In
##' either window, left mouse clicks define the start of the selected
##' subset and right mouse clicks defines the end of the selected
##' subset. The "a" key accepts the current selection an proceeds to
##' the next stage.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits\cr
##' 'a' \tab Accepts any changes and advances to the next stage \cr
##' '+'/'-' \tab Zoom in or out \cr
##' }
##'
##' In the second stage the user guides an semi-automated search for
##' the twilight times.  Again, a right mouse click in the first
##' window produces a zoomed image of the clicked region in the
##' second.  In the second window, the user must click in periods of
##' nights.  Left mouse clicks select regions to be searched for
##' twilights, right mouse clicks select regions that should not be
##' searched.  The "u" key provides an undo facility, allowing mouse
##' clicks to be deleted, and the "a" key terminates the search and
##' proceeds to the next stage.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits, returning the dataframe of edited twilights\cr
##' 'a' \tab Accepts any changes and advances to the next stage \cr
##' 'b' \tab Returns back to the previous stage \cr
##' 'u' \tab Removes the most recent search point\cr
##' '+'/'-' \tab Zoom in or out \cr
##' }
##'
##' In the third stage the user may insert additional twilights where
##' the light record is incomplete. As for the first two stages, right
##' mouse clicks in the first window produce a zoomed image of the
##' clicked region in the second window.  In the second window, the
##' user may ad times of sunrise or sunset.  A left mouse click adds a
##' sunset, and a right mouse click adds a sunrise. The numeric keys
##' on the keyboard mark the most recent inerted point with a "mark" 0
##' to 9. If the user provides both the \code{zenith} argument and a
##' two column matrix of locations for the \code{fixed} argument,
##' marking a point with the digit n replaces the marked twilight with
##' the computed twilight time for the n-th fixed location.  Again "u"
##' provides an undo facility and the "a" key proceeds to the next
##' stage.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits, returning the dataframe of edited twilight segments \cr
##' 'a' \tab Accepts any changes and advances to the next stage \cr
##' 'b' \tab Returns back to the previous stage \cr
##' 'u' \tab Removes the most recent insertion \cr
##' '+'/'-' \tab Zoom in or out \cr
##' '0'-'9' \tab Mark the most recent addition \cr
##' }
##'
##' In the fourth stage, the user may adjust individual twilights
##' based on the observed light profile.  The first window now shows
##' the sequence of twilights.  A right mouse click on any twilight
##' shows the light profile for that twilight in the second window,
##' together with the profiles for the preceeding and following days.
##' A left click in the second window proposes a new location for the
##' current twilight, but no change is made until the edit is accepted
##' with the "a" key. A right click in the second window marks the
##' twilight as deleted, and the numeric keys mark the twilight with a
##' "mark" 0 to 9. If the user provides both the \code{zenith}
##' argument and a two column matrix of locations for the \code{fixed}
##' argument, marking a point with the digit n replaces the marked
##' twilight with the computed twilight time for the n-th fixed
##' location.  The "u" key reverts the edits to the current twilight.
##'
##' In either window
##' \tabular{ll}{
##' 'q' \tab Quits, returning the dataframe of edited twilight segments \cr
##' 'a' \tab Accepts the candidate edit \cr
##' 'b' \tab Returns back to the previous stage \cr
##' 'd' \tab Toggle deletion of this twilight \cr
##' 'i' \tab Toggles the display of the light image \cr
##' 'p' \tab Toggles the display of individual points \cr
##' 'p' \tab Toggles the display of surrounding profiles \cr
##' 'r' \tab Resets the selection \cr
##' 'u' \tab Revert changes to this twilight \cr
##' '+'/'-' \tab Zoom in or out \cr
##' 'Left arrow' \tab Jump to previous twilight \cr
##' 'Right arrow' \tab Jump to ext twilight \cr
##' '0'-'9' \tab Mark this twilight \cr
##' }
##'
##' If \code{map=TRUE} an additional window is constructed that shows
##' the approximate track based on hte current estimated twilights. A
##' function to draw a background map for this window can be supplied
##' with the \code{plotMap} argument.  This must be ia function of two
##' arguments \code{xlim} and \code{ylim} the determine the extent of
##' the plot. A left mouse click in this window selects the twilight
##' associated with the nearest location, a right mouse click in this
##' window displays the time of twilight at the selected location for
##' the currently selected date as an edit in the light profile
##' window, and this edit can be accepted with the "a" key.
##'
##' The results of a previous call of \code{preprocessLight} can be
##' updated by providing the dataframe generated by the previous call
##' as the \code{twilights} argument, and using the \code{stage}
##' argument to determine to point from which to restart the process.
##'
##' @title Interactively derive twilight
##' @param tagdata a dataframe with columns \code{Date} and
##' \code{Light} that are the sequence of sample times (as POSIXct)
##' @param threshold the light threshold that defines twilight.
##' @param offset the starting hour for the vertical axes.
##' @param lmax the maximum light level to plot.
##' @param zlim the range of light levels to plot in images.
##' @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.
##' @param zenith the solar zenith angle that defines twilight.
##' @param fixed a two column array of fixed locations
##' @param map whether to plot an approximate path
##' @param plotMap a function to plot a background map when \code{map=TRUE}.
##' @param twilights a result of a previous run to re-edit
##' @param stage when re-editing, the stage to commence at
##' @param point.cex expansion factor for plot points.
##' @param width width of the interface windows.
##' @param height height of the interface windows.
##' @param palette a colour palette of 7 colours.
##' @return A dataframe with columns
##' \item{\code{Twilight}}{times of twilight}
##' \item{\code{Rise}}{logical indicating sunrise}
##' \item{\code{Deleted}}{logical indentifying deleted twilights}
##' \item{\code{Marker}}{integer vector of marks}
##' \item{\code{Inserted}}{logical vector of marks}
##'\item{\code{Twilight3}}{Twilight from stage 3}
##' \item{\code{Marker3}}{Marker from stage 3}
##' where each row corresponds to a single twilight.
##' @importFrom SGAT twilight thresholdPath
##' @importFrom stats median
##' @importFrom graphics abline grconvertX grconvertY lines plot.new points title
##' @importFrom grDevices dev.cur dev.new dev.off dev.set getGraphicsEvent setGraphicsEventHandlers recordPlot replayPlot
##' @export
preprocessLight <- function(tagdata,threshold,offset=0,lmax=64,zlim=c(0,lmax),
                            extend=0,dark.min=0,
                            zenith=96,fixed=NULL,
                            map=FALSE,plotMap=NULL,
                            twilights=NULL,stage=1,
                            point.cex=0.8,width=12,height=4,
                            palette=defaultPalette[c(5,2,9,3,4,1,13)]) {
  
  ## Default is no map
  if(is.null(plotMap))
    plotMap <- function(xlim,ylim) { plot.new(); plot.window(xlim,ylim) }

  ## Round down/up to nearest offset
  floorDate <- function(date) date - ((as.hour(date)-offset)%%24)*60*60
  ceilingDate <- function(date) date + ((offset-as.hour(date))%%24)*60*60

  ## Set initial selection
  minDate <- if(is.null(attr(twilights,"minDate"))) floorDate(min(tagdata$Date)) else attr(twilights,"minDate")
  maxDate <- if(is.null(attr(twilights,"maxDate"))) ceilingDate(max(tagdata$Date)) else attr(twilights,"maxDate")

  ## Set initial seed points
  seed <- attr(twilights,"seed")
  include <- attr(twilights,"include")

  Date1 <- minDate
  Date2 <- maxDate
  DateZ <- NULL

  zoom <- 6

  ## Cached data subsets
  index <- 1
  editpt <- NULL
  changed <- FALSE
  twls <- NULL
  date <- vector(3,mode="list")
  lght <- vector(3,mode="list")
  showobs <- FALSE
  showimg <- TRUE
  showlag <- TRUE
  rmap <- NULL
  path <- NULL

  ## Set cached values
  cache <- function(k) {
    index <<- k
    editpt <<- NULL
    changed <<- FALSE
    ## Get twilight times
    twl <- twilights$Twilight[index]
    keep <- (twilights$Twilight >= twl-3600*zoom) & (twilights$Twilight <= twl+3600*zoom)
    keep[index] <- FALSE
    twls <<- twilights$Twilight[keep]
    ## Get profiles
    for(k in 1:3) {
      mid <- twilights$Twilight[index]+(k-2)*86400
      keep <- (tagdata$Date >= mid-3600*zoom) & (tagdata$Date <= mid+3600*zoom)
      date[[k]] <<- tagdata$Date[keep]
      lght[[k]] <<- tagdata$Light[keep]
    }
  }



  ## Select device
  setDevice <- function(dev) if(dev.cur()!=dev) dev.set(dev)
  ## Focus if possible
  focus <- if(exists("bringToTop",mode="function")) grDevices::bringToTop else setDevice

  ## Show vertical lines in gaps (stage 3)
  gapline <- function(ts,col) {
    dt <- diff(ts)
    ks <- which(dt > 36*60*60)
    if(length(ks) > 0)
      abline(v=.POSIXct(ts[ks]+dt[ks]/2,"GMT"),lwd=2,col=col)
  }

  ## Draw the twilights window
  winADraw <- function() {
    setDevice(winA)
    imageDraw(if(showimg) tagdata,
               twilights,offset=offset,mark=index,
               points.col=palette[ifelse(twilights$Deleted,7,ifelse(twilights$Rise,1,2))],
               zlim=zlim,point.cex=point.cex)
    title(main=switch(stage,
            "Select subset",
            "Find twilights",
            "Insert twilights",
            "Edit twilights"))

    if(stage==1)
      selectionRectangle(Date1,Date2,col=palette[6])

    if(!is.null(DateZ) && (stage==2 || stage==3)) {
      start <- floorDate(max(DateZ - zoom*24*60*60,minDate))
      end <- ceilingDate(min(start+2*zoom*24*60*60,maxDate))
      start <- floorDate(max(end - 2*zoom*24*60*60,minDate))
      selectionRectangle(start,end,col=palette[6])
    }


    if(stage==2 && length(include)>0)
      tsimagePoints(seed,offset=offset,pch=16,col=palette[ifelse(include,4,5)])

    if(stage==3) {
      gapline(sort(as.numeric(twilights$Twilight[twilights$Rise])),col=palette[1])
      gapline(sort(as.numeric(twilights$Twilight[!twilights$Rise])),col=palette[2])
    }

  }


  winBDraw <- function() {
    setDevice(winB)
    if(stage <= 3) {
      ## Display zoomed image
      xlim <- if(!is.null(DateZ)) {
        start <- floorDate(max(DateZ - zoom*24*60*60,minDate))
        end <- ceilingDate(min(start+2*zoom*24*60*60,maxDate))
        start <- floorDate(max(end - 2*zoom*24*60*60,minDate))
        c(start,end)
      }
      imageDraw(tagdata,twilights,offset=offset,
                 points.col=palette[ifelse(twilights$Rise,1,2)],
                 zlim=zlim,xlim=xlim)
      title(main=as.character(DateZ))
    } else {
      ## Plot light profiles
      marker <- twilights$Marker[index]
      profileInit(date,lght,lmax=lmax,
                  xlab=if(marker>0) paste("Marker: ",marker) else "",
                  main=as.character(twilights$Twilight[index]))
      ## Overlay with light profiles
      profileOverlay(date,lght,threshold,showobs,showlag,
                     profile.col=palette[3:5],threshold.col=palette[6],
                     point.cex=point.cex)
      abline(v=twls,col=palette[7])
      abline(v=twilights$Twilight[index],col=if(!twilights$Deleted[index]) palette[6] else palette[7])
      if(changed) points(editpt[1],editpt[2],pch=16,col=palette[6])
    }

    if(stage==1)
      tsimagePoints(c(Date1,Date2),offset=offset,pch=16,col=palette[6])

    if(stage==2 && length(include)>0)
      tsimagePoints(seed,offset=offset,pch=16,col=palette[ifelse(include,4,5)])

    if(stage==3) {
      gapline(sort(as.numeric(twilights$Twilight[twilights$Rise])),col=palette[1])
      gapline(sort(as.numeric(twilights$Twilight[!twilights$Rise])),col=palette[2])
    }
  }


  ## Draw path
  winCDraw <- function() {
    if(!is.null(path)) {
      setDevice(winC)
      ## Create underlying map
      if(is.null(rmap)) {
        ## User defined map function
        plotMap(range(path[,1]),range(path[,2]))
        rmap <<- recordPlot()
      } else {
        ## Replot stored plot
        replayPlot(rmap)
      }
      ## Show full path
      lines(path[,1],path[,2],col=palette[5])
      points(path[,1],path[,2],col=palette[5],pch=16,cex=0.4)
      ## Highlight current point
      marker <- twilights$Marker[index]
      points(path[index,1],path[index,2],col=palette[if(marker>0 & marker <= NROW(fixed)) 4 else 6],pch=16,cex=1)
    }
  }

  ## onMouseDown callback for twilights window.
  winAOnMouseDown <- function(buttons,x,y) {
    setDevice(winA)
    if(length(buttons) > 0) {
      b <- mouseButton(buttons)

      if(stage==1) {
        ## Set endpoints and zoom window B
        if(b==1) {
          Date1 <<- ndcTsimageDate(x,y)
          if(Date1 > Date2) Date2 <<- maxDate
          DateZ <<- Date1
        }
        if(b==2) {
          Date2 <<- ndcTsimageDate(x,y)
          if(Date2 < Date1) Date1 <<- minDate
          DateZ <<- Date2
        }
      }

      if(stage==2 || stage==3) {
        ## Zoom window B.
        if(b==1) {
          DateZ <<- ndcTsimageDate(x,y)
        }
      }

      if(stage==4) {
        ## Determine selected profile.
        if(b==1) {
          day <- twilights$Twilight
          hour <- hourOffset(as.hour(twilights$Twilight),offset)
          k <- (ndcClosest(x,y,c(day,day,day),c(hour-24,hour,hour+24))-1)%%length(day)+1
          ## Redraw
          cache(k)
        }
      }
    }

    winADraw()
    winBDraw()
    if(map && stage==4) winCDraw()
    NULL
  }


  ## onMouseDown callback for profile window
  winBOnMouseDown <- function(buttons,x,y) {
    setDevice(winB)

    if(length(buttons) > 0) {
      b <- mouseButton(buttons)
      if(stage==1) {
        ## Set endpoint and zoom
        if(b==1) {
          Date1 <<- ndcTsimageDate(x,y)
          if(Date1 > Date2) Date2 <<- maxDate
        }
        if(b==2) {
          Date2 <<- ndcTsimageDate(x,y)
          if(Date2 < Date1) Date1 <<- minDate
        }
      }

      if(stage==2) {
        ## Add segments to include or exclude
        if(b==1 || b==2) {
          seed <<- c(ndcTsimageDate(x,y),seed)
          include <<- c(buttons[1]==0,include)

          ## Recompute twilights
          twilights <<- findTwilights(tagdata,threshold=threshold,
                                       include=seed[include],exclude=seed[!include],
                                       extend=extend,dark.min=dark.min)
          twilights$Deleted <<- logical(nrow(twilights))
          twilights$Marker <<- integer(nrow(twilights))
          twilights$Inserted <<- logical(nrow(twilights))

        }
      }

      if(stage==3) {
        ## Add twilights
        if(b==1) {
          twilights <<- rbind(twilights,
                              data.frame(Twilight=ndcTsimageDate(x,y),
                                         Rise=FALSE,
                                         Deleted=FALSE,
                                         Marker=0,
                                         Inserted=TRUE))
        }
        if(b==2) {
          twilights <<- rbind(twilights,
                              data.frame(Twilight=ndcTsimageDate(x,y),
                                         Rise=TRUE,
                                         Deleted=FALSE,
                                         Marker=0,
                                         Inserted=TRUE))
        }
      }

      if(stage==4) {
        ## Button 1 -> record location
        if(b==1) {
          changed <<- TRUE
          editpt <<- c(grconvertX(x,from="ndc",to="user"),
                       grconvertY(y,from="ndc",to="user"))
        }
        ## Button 2 -> toggle deletion
        if(b==2) {
          twilights$Deleted[index] <<- !twilights$Deleted[index]
          cache(index)
        }
      }
    }
    winADraw()
    winBDraw()
    if(map && stage==4) winCDraw()
    NULL
  }

  ## onMouseDown callback for map window.
  winCOnMouseDown <- function(buttons,x,y) {
    if(!is.null(path)) {
      setDevice(winC)
      if(length(buttons) > 0) {
        b <- mouseButton(buttons)
        ## Button 1 -> select point
        if(b==1) {
          ## Select nearest point
          k <- ndcClosest(x,y,path[,1],path[,2])
          cache(k)
        }
        if(b==2) {
          lon <- grconvertX(x,from="ndc",to="user")%%360
          lat <- pmin(90,pmax(-90,grconvertY(y,from="ndc",to="user")))
          twl <- twilight(twilights$Twilight[index],lon,lat,rise=twilights$Rise[index],zenith=zenith)
          changed <<- TRUE
          editpt <<- c(twl,threshold)
        }

      }

      winADraw()
      winBDraw()
      winCDraw()
    }
    NULL
  }

  ## onKeybd callback for both windows
  onKeybd <- function(key) {
    ## Common keybindings
    ## q : quit
    if(key=="q") {
      if(!is.null(twilights)) {
        twilights <<- twilights[order(twilights$Twilight),]
        attr(twilights,"interval") <- c(minDate,maxDate)
      }
      return(-1)
    }
    ## +/- : zoom time window around threshold crossing
    if(key=="+") {
      zoom <<- max(zoom-2,1)
      if(stage==4) cache(index)
    }
    if(key=="-") {
      zoom <<- min(zoom+2,36)
      if(stage==4) cache(index)
    }

    ## Left/Right movement of zoomed image
    if(stage < 4 && !is.null(DateZ)) {
      if(key=="Left") {
        DateZ <<- .POSIXct(DateZ-zoom*24*60*60,"GMT")
        if(DateZ < minDate) DateZ <<- minDate
      }
      if(key=="Right") {
        DateZ <<- .POSIXct(DateZ+zoom*24*60*60,"GMT")
        if(DateZ > maxDate) DateZ <<- maxDate
      }
    }

    ## i : toggle image display
    if(stage > 1 && key=="i") {
      showimg <<- !showimg
    }

    ## b : go back a stage
    if(stage > 1 && key=="b") {
      stage <<- stage-1
    }




    if(stage==1) {
      ## Stage 1 keybindings
      ## a : accept and advance to next stage
      if(key=="a") {
        minDate <<- Date1
        maxDate <<- Date2
        tagdata <<- tagdata[tagdata$Date >= minDate & tagdata$Date <= maxDate,]
        DateZ <<- NULL
        stage <<- 2
      }
    } else if(stage==2) {
      ## Stage 2 keybindings
      ## a : accept and advance to next stage
      if(key=="a") {
        index <<- 1
        DateZ <<- NULL
        stage <<- 3
      }
      if(key=="u") {
        if(length(seed)>0) {
          seed <<- seed[-1]
          include <<- include[-1]

          ## Recompute twilights
          twilights <<- findTwilights(tagdata,threshold=threshold,
                                       include=seed[include],exclude=seed[!include],
                                       extend=extend,dark.min=dark.min)
          twilights$Deleted <<- logical(nrow(twilights))
          twilights$Marker <<- integer(nrow(twilights))
          twilights$Inserted <<- logical(nrow(twilights))
        }
      }
    } else if(stage==3) {
      ## Stage 3 keybindings
      ## a : accept and advance to next stage
      if(key=="a") {
        twilights$Twilight3 <<- twilights$Twilight
        twilights$Marker3 <<- twilights$Marker
        ks <- order(twilights$Twilight)
        twilights <<- twilights[ks,]
        cache(ks[index])
        if(map) path <<- thresholdPath(twilights$Twilight,twilights$Rise,zenith=zenith)$x
        zoom <<- 6
        stage <<- 4
      }
      ## u : undo twilight insertions
      if(key=="u") {
        n <- nrow(twilights)
        if(twilights$Inserted[n]) twilights <<- twilights[-n,]
        index <<- nrow(twilights)
      }
      ## Markers
      if(key >= "0" && key <= "9") {
        marker <- as.numeric(key)
        index <<- nrow(twilights)
        twilights$Marker[index] <<- marker
        if(marker > 0 && marker <= NROW(fixed))
          twilights$Twilight[index] <<- twilight(twilights$Twilight[index],
                                                 fixed[marker,1],fixed[marker,2],
                                                 rise=twilights$Rise[index],
                                                 zenith=zenith)
      }
    } else if(stage==4) {
      ## Stage 4 keybindings
      ## a : accept changes
      if(key=="a") {
        if(!is.null(editpt)) {
          twilights$Twilight[index] <<- .POSIXct(editpt[1],"GMT")
          if(map) path <<- thresholdPath(twilights$Twilight,twilights$Rise,zenith=zenith)$x
          editpt <<- NULL
          changed <<- FALSE
        }
      }
      if(key=="m") {
        ks <- which((twilights$Twilight <= twilights$Twilight[index]+60*60*60) &
                      (twilights$Twilight >= twilights$Twilight[index]-60*60*60) &
                        (twilights$Rise==twilights$Rise[index]))
        mdn <- median((as.numeric(twilights$Twilight[ks])-
                         as.numeric(twilights$Twilight[index])+12*60*60)%%(24*60*60)-12*60*60)
        twilights$Twilight[index] <<- .POSIXct(twilights$Twilight[index]+mdn,"GMT")
        if(map) path <<- thresholdPath(twilights$Twilight,twilights$Rise,zenith=zenith)$x
      }
      ## x : reset selection
      if(key=="r") {
        cache(index)
      }
      ## d : toggle deletion
      if(key=="d") {
        twilights$Deleted[index] <<- !twilights$Deleted[index]
        if(map) path <<- thresholdPath(twilights$Twilight,twilights$Rise,zenith=zenith)$x
      }
      ## p : toggle display of points in the profile window
      if(key=="p") {
        showobs <<- !showobs
      }
      ## l : toggle display of lagged twilights in the profile window
      if(key=="l") {
        showlag <<- !showlag
      }
      ## u : undo edits
      if(key=="u") {
        twilights$Twilight[index] <<- twilights$Twilight3[index]
        twilights$Marker[index] <<- twilights$Marker3[index]
        if(map) path <<- thresholdPath(twilights$Twilight,twilights$Rise,zenith=zenith)$x
      }
      ## Left/Right : jump to neighbouring twilight
      if(key=="Left") {
        cache(max(index-1,1))
      }
      if(key=="Right") {
        cache(min(index+1,nrow(twilights)))
      }
      ## Markers
      if(key >= "0" && key <= "9") {
        marker <- as.numeric(key)
        twilights$Marker[index] <<- marker
        if(marker > 0 && marker <= NROW(fixed))
          twilights$Twilight[index] <<- twilight(twilights$Twilight[index],
                                                 fixed[marker,1],fixed[marker,2],
                                                 rise=twilights$Rise[index],
                                                 zenith=zenith)
      }
    }

    ## Redraw
    winADraw()
    winBDraw()
    if(map && stage==4) winCDraw()
    NULL
  }


  ## Set up the stage
  if(is.null(twilights)) stage <- 1
  if(stage > 1)
    tagdata <- tagdata[tagdata$Date >= minDate & tagdata$Date <= maxDate,]
  if(stage==4) cache(1)

  ## Set up twilights window
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winA <- dev.cur()
  winADraw()
  setGraphicsEventHandlers(
    which=winA,
    prompt="Selection Window",
    onMouseDown=winAOnMouseDown,
    onKeybd=onKeybd)
  ## Set up profile window
  dev.new(width=width,height=height,noRStudioGD=TRUE)
  winB <- dev.cur()
  winBDraw()
  setGraphicsEventHandlers(
    which=winB,
    prompt="Edit Window",
    onMouseDown=winBOnMouseDown,
    onKeybd=onKeybd)
  if(map) {
    dev.new(width=8,height=8,noRStudioGD=TRUE)
    winC <- dev.cur()
    plotMap(c(0,360),c(-90,90))
    setGraphicsEventHandlers(
      which=winC,
      prompt="Edit Window",
      onMouseDown=winCOnMouseDown,
      onKeybd=onKeybd)
  }
  focus(winA)
  ## Monitor for events
  tryCatch({
    getGraphicsEvent()
    dev.off(winB)
    dev.off(winA)
    if(map) dev.off(winC)
  }, finally={
    dev.off(winB)
    dev.off(winA)
    if(map) dev.off(winC)})
  attr(twilights,"minDate") <- minDate
  attr(twilights,"maxDate") <- maxDate
  attr(twilights,"seed") <- seed
  attr(twilights,"include") <- include
  twilights
}