R/Raster.R

In SWotherspoon/SGAT: Solar/Satellite Geolocation for Animal Tracking

##' Bin locations to form a 2D density raster
##'
##' Constructs either a 2D histogram or kernel density estimate from
##' the samples for a sequence of locations.  Optionally, a vector of
##' weights can be provided to differentially weight samples by
##' location.
##'
##' @title Location Density Raster
##' @param s a single chain or a list of parallel chains generated by
##' \code{estelleMetropolis} or \code{stellaMetropolis}.
##' @param grid raster object that defines the sampling grid
##' @param weights weights for each location
##' @param bw optional bandwidth for the kernel density estimator.
##' @param zero.is.na if \code{TRUE}, zero counts are returned as \code{NA}.
##' @return a raster representing a 2D histogram of locations
##' @importFrom raster raster isLonLat values bbox
##' @export
locationRasterize <- function(s,grid,weights=1,zero.is.na=TRUE) {
  s <- chainCollapse(s)
  weights <- rep(weights,length=dim(s)[1L])
  r <- if(is.null(grid)) raster() else raster(grid)

  ## Project coords
  if(!is.na(projection(r)) && !isLonLat(r)) {
    from <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
    to <- CRS(projection(r))
    p <- cbind(as.vector((s[,1L,]+180)%%360-180),as.vector(s[,2L,]))
    p <- coordinates(spTransform(SpatialPoints(p,from),to))
    s[,1L,] <- p[,1L]
    s[,2L,] <- p[,2L]
  }

  ## Do the binning
  bb <- bbox(r)
  nx <- ncol(r)
  ny <- nrow(r)
  xbin <- seq.int(bb[1L,1L],bb[1L,2L],length.out=nx+1L)
  ybin <- seq.int(bb[2L,1L],bb[2L,2L],length.out=ny+1L)

  A <- 0
  W <- dim(s)[3L]*prod(res(r))
  for(k in seq_len(dim(s)[1L])) {
    A <- A+(weights[k]/W)*table(
      factor((ny+1)-.bincode(s[k,2L,],ybin,TRUE,TRUE),levels=seq_len(ny)),
      factor(.bincode(s[k,1L,],xbin,TRUE,TRUE),levels=seq_len(nx)))
  }
  if(zero.is.na) A[A==0] <- NA
  values(r) <- A
  r
}



##' @rdname locationRasterize
##' @importFrom stats dnorm var
##' @importFrom raster raster xFromCol yFromRow values<-
##' @importFrom sp CRS SpatialPoints
##' @export
locationKernelize <- function(s,grid,weights=1,bw=NULL) {

  bandwidth <- function (x) {
    x <- as.vector(x)
    1.06*min(sqrt(var(x)),diff(quantile(x,c(0.25,0.75)))/1.34)*length(x)^(-1/5)
  }

  s <- chainCollapse(s)
  weights <- rep(weights,length=dim(s)[1L])
  r <- if(is.null(grid)) raster() else raster(grid)

  ## Project coords
  if(!is.na(projection(r)) && !isLonLat(r)) {
    from <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
    to <- CRS(projection(r))
    p <- cbind(as.vector((s[,1L,]+180)%%360-180),as.vector(s[,2L,]))
    p <- coordinates(spTransform(SpatialPoints(p,from),to))
    s[,1L,] <- p[,1L]
    s[,2L,] <- p[,2L]
  }

  ## Evaluate kernels
  bw <- if(is.null(bw)) c(bandwidth(s[,1L,]),bandwidth(s[,2L,])) else rep(bw/4,length.out=2L)
  xcell <- xFromCol(r)
  ycell <- yFromRow(r)

  A <- 0
  W <- dim(s)[3L]*prod(bw)
  for(k in seq_len(dim(s)[1L])) {
      A <- A+(weights[k]/W)*tcrossprod(dnorm(outer(ycell,s[k,2L,],"-")/bw[2L]),
                                       dnorm(outer(xcell,s[k,1L,],"-")/bw[1L]))
  }
  values(r) <- A
  r
}



##' Generate density estimates for time slices of a trip
##'
##' These functions divide a trip into time slices, and generate
##' location density estimates for each time slice.
##'
##' The \code{slices} function defines the slices into which the trip
##' is divided, and whether density estimates are generated for the
##' primary or intermediate locations. If \code{breaks} is NULL, each
##' location forms a separate time slice, otherwise \code{breaks}
##' divides the trip into time slices in the same style as
##' \code{\link{cut.POSIXt}}.  A default set of samples and a raster
##' defining the spatial grid may also be specified.
##'
##' The \code{slice} function generates a density estimate for a time
##' slice, \code{sliceInterval} returns the corresponding time
##' interval, and \code{sliceIndices} returns the indices that will
##' yield non null raster.
##'
##' @title Location Density Estimates
##' @param type construct density rasters for \code{primary} (x) or
##' \code{intermediate} (z) locations
##' @param slices an object generated by \code{slices}
##' @param k the index of the time slice to bin.
##' @param breaks NULL, or a specification suitable for
##' \code{cut.POSIXt} to define the time slices to bin.
##' @param mcmc object generated by generated by
##' \code{estelleMetropolis} or \code{stellaMetropolis}.
##' @param grid raster object that defines the sampling grid.
##' @param weights weights for each location.
##' @param include.lowest parameter to \code{cut.POSIXt}.
##' @param right parameter to \code{cut.POSIXt}.
##' @param method return a kernel density estimate (\code{"kde"})
##' or simple histogram (\code{"hist"})?
##' @param bandwidth bandwidth of the kernel density estimator.
##' @param chains NULL or the subset of chains to bin.
##' @return \code{slice} returns the locations for time slice of the
##' track binned into a raster, \code{slices} returns a slices
##' object that defines the time slices into which to bin,
##' \code{sliceInterval} returns the time interval spanned by a
##' slice, and \code{sliceIndices} returns the valid set of indices
##' that will yield a raster.
##' @importFrom raster raster
##' @export
slice <- function(slices,k,
                  mcmc=slices$mcmc,grid=slices$grid,weights=slices$weights,
                  method=slices$method,bandwidth=slices$bandwidth,chains=NULL) {
  ## Split times
  time <- mcmc$model$time
  if(!is.null(slices$breaks)) {
    ks <- unclass(cut(if(slices$type=="intermediate") time[-length(time)] else time,
                      breaks=slices$breaks,
                      include.lowest=slices$include.lowest,
                      right=slices$right))
    k <- which(ks %in% k)
  }
  if(length(k)>0) {
    ## Select x or z
    if(slices$type=="intermediate") {
      w <- if(!is.null(weights)) weights else diff(as.numeric(time)/(60*60))
      s <- chainCollapse(mcmc$z,chains=chains)
    } else {
      w <- if(!is.null(weights)) weights else rep(1,length(time))
      s <- chainCollapse(mcmc$x,chains=chains)
    }

    ## Subset
    w <- w[k]
    s <- s[k,,,drop=FALSE]

    ## Convert
    if(method=="hist")
      locationRasterize(s,grid,w)
    else
      locationKernelize(s,grid,w,bandwidth)
  }
}


##' @rdname slice
##' @export
slices <- function(type=c("primary","intermediate"),
                   breaks=NULL,mcmc=NULL,grid=raster(),weights=NULL,
                   method=c("hist","kde"),bandwidth=0,
                   include.lowest=TRUE,right=FALSE) {
  r <- list(type=match.arg(type),
            breaks=breaks,
            mcmc=mcmc,
            grid=grid,
            weights=weights,
            method=match.arg(method),
            bandwidth=bandwidth,
            include.lowest=include.lowest,
            right=right)
  class(r) <- "slices"
  r
}


##' @rdname slice
##' @export
sliceInterval <- function(slices,k,mcmc=slices$mcmc) {
  time <- mcmc$model$time
  if(!is.null(slices$breaks)) {
    ks <- unclass(cut(if(slices$type=="intermediate") time[-length(time)] else time,
                      breaks=slices$breaks,
                      include.lowest=slices$include.lowest,
                      right=slices$right))
    k <- which(ks %in% k)
  }
  if(length(k)>0) range(time[if(slices$type=="intermediate") c(k,k+1L) else k])
}

##' @rdname slice
##' @export
sliceIndices <- function(slices,mcmc=slices$mcmc) {
  time <- mcmc$model$time
  if(slices$type=="intermediate") time <- time[-length(time)]
  if(!is.null(slices$breaks))
    unique(unclass(cut(time,
                       breaks=slices$breaks,
                       include.lowest=slices$include.lowest,
                       right=slices$right)))
  else
    seq_along(time)
}


##' Extract longitude and latitude of raster cells.
##'
##' Extract the longitude and latitude of the center of the requested
##' cells of a Raster* object, similar to \code{xyFromCell} and \code{cellFromXY}.
##' @title Raster cell longitude and latitudes
##' @param raster a raster object
##' @param cells the cell numbers
##' @param pts the lon, lat locations as a 2 column matrix
##' @param spatial return locations as SpatialPoints object instead of a matrix.
##' @return \code{lonlatFromCell} returns the lon,lat locations for
##' the requested cells as a 2 column matrix, and
##' \code{cellFromLonLat} returns the cells corresponding to the 2
##' column matrix of lon,lat positions.
##' @importFrom raster xyFromCell cellFromXY
##' @importFrom sp CRS spTransform coordinates
##' @export
lonlatFromCell <- function(raster,cells,spatial=FALSE) {
  if(is.na(projection(raster)) || isLonLat(raster)) {
    xyFromCell(raster,cells,spatial=spatial)
  } else {
    p <- spTransform(xyFromCell(raster,cells,spatial=TRUE),
                     CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
    if(spatial) p else coordinates(p)
  }
}


##' @rdname lonlatFromCell
##' @export
cellFromLonLat <- function(raster,pts) {
  if(!(is.na(projection(raster)) || isLonLat(raster))) {
    pts <- SpatialPoints(pts,proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
    pts <- coordinates(spTransform(pts,projection(raster,FALSE)))
  }
  cellFromXY(raster,pts)
}



##' Spatial maps of twilight residuals
##'
##' This function calculates the twilight residuals corresponding to
##' a single observed twilight across a grid of locations.
##'
##' @title Twilight Residuals
##' @param twilight an observed time of twilight
##' @param rise \code{TRUE} if twilight is a sunrise
##' @param grid raster object that define the sampling grid.
##' @param zenith the solar zenith angle that defines twilight.
##' @return a raster of twilight residuals (in minutes).
##' @seealso \code{\link{twilightResiduals}}
##' @importFrom raster raster
##' @export
twilightResidualsMap <- function(twilight,rise,grid,zenith=96) {
  p <- lonlatFromCell(grid,seq_len(ncell(grid)))
  legal <- !(is.na(p[,1]) | is.na(p[,2]))
  sgn <- if(rise) 1 else -1
  s <- solar(twilight)
  r <- raster(grid)
  r[legal] <- 4*sgn*(s$solarTime-twilightSolartime(s,p[legal,1L],p[legal,2L],rise,zenith))
  r
}



## ##' Simple Land Mask Example
## ##'
## ##' This function provides a low resolution land/sea mask based on
## ##' wrld_simpl from \pkg{maptools}.  While this provides a good
## ##' general starting point for many analyses, users are encouraged to
## ##' develop higher resolution masks specific to their analyses.
## ##'
## ##' This implementation discretizes wrld_simpl onto a simple grid in
## ##' longitude and latitude. The bounding box is set by \code{xlim} and
## ##' \code{ylim}.  In general \code{xlim} must be contained in the
## ##' interval [-540,540]. If \code{xlim} is \code{NULL} or spans more
## ##' than 360 degrees, then \code{xlim} is taken to be [-180,180] and
## ##' longitudes outside this range will be mapped back into this
## ##' region.  Similarly, if \code{ylim} is \code{NULL}, it is taken to
## ##' be [-90,90].
## ##'
## ##' For points outside the bounding box, the mask function returns
## ##' \code{NA}. For points inside the bounding box, the function
## ##' returns \code{land} for points on land (within the resolution of
## ##' the grid) and \code{!land} for points at sea.
## ##'
## ##' @title Example Land Mask
## ##' @param xlim the longitude range.
## ##' @param ylim the latitude range.
## ##' @param n number of grid cells per degree.
## ##' @param land mask for land or sea?
## ##' @return A vectorized function that determines (approximately)
## ##' whether a location is land or sea.
## ##' @examples
## ##'
## ##'
## ##' ## Define mask for sea
## ##' is.sea <- land.mask(xlim=c(80,170),ylim=c(-70,-10),n=4,land=FALSE)
## ##' ## Land
## ##' is.sea(cbind(120,-30))
## ##' ## Sea
## ##' is.sea(cbind(120,-60))
## ##' ## Out of bounding box
## ##' is.sea(cbind(0,0))
## ##' @importFrom maptools elide
## ##' @export
## ## Function for constructing a land/sea mask from wrld_simpl
## land.mask <- function(xlim=NULL,ylim=NULL,n=4,land=TRUE) {
##   data("wrld_simpl",package="maptools",envir=environment())
##   ## Set limits
##   wrap360 <- is.null(xlim) || diff(xlim)>=360
##   if(wrap360) xlim <- c(-180,180)
##   if(is.null(ylim)) ylim <- c(-90,90)
##
##   ## Construct lookup raster
##   r <- raster(nrows=n*diff(ylim),ncols=n*diff(xlim),
##               xmn=xlim[1],xmx=xlim[2],ymn=ylim[1],ymx=ylim[2],
##               crs=proj4string(wrld_simpl))
##   if(xlim[1] >= -180 && xlim[2] <= 180) {
##     r <- rasterize(wrld_simpl,r,1,silent=TRUE)
##   } else if(xlim[1] < -180) {
##     r <- cover(rasterize(maptools::elide(wrld_simpl,shift=c(-360,0)),r,1,silent=TRUE),
##                rasterize(wrld_simpl,r,1,silent=TRUE))
##   } else {
##     r <- cover(rasterize(wrld_simpl,r,1,silent=TRUE),
##                rasterize(maptools::elide(wrld_simpl,shift=c(360,0)),r,1,silent=TRUE))
##   }
##   r <- as.matrix(is.na(r))[nrow(r):1,]
##   if(land) r <- !r
##
##   ## Lookup bins
##   xbin <- seq(xlim[1],xlim[2],length=ncol(r)+1)
##   ybin <- seq(ylim[1],ylim[2],length=nrow(r)+1)
##
##   if(wrap360) {
##     function(p) {
##       r[cbind(.bincode(p[,2],ybin),.bincode((p[,1]-180)%%360+180,xbin))]
##     }
##   } else {
##     function(p) {
##       r[cbind(.bincode(p[,2],ybin),.bincode(p[,1],xbin))]
##     }
##   }
## }