R/xytract.R

In mdsumner/xytract: What the Package Does (One Line, Title Case)

## copied from https://github.com/AustralianAntarcticDivision/raadtools/blob/822b4340cc5bc63b3a03a01c6832bd3da9fb932b/R/extract.R
## 2018-06-09

setOldClass("trip")
.read.generic  <- function(x, y, ...) {
  ## read function "x", takes "y" as Date, POSIXct, character
  x(y, ...)
}

.standard.assumeXYT.TimeError <- function() {
  stop("invalid times in data, ensure that y is a data.frame with values of longitude, latitude, and date-times")
}


.determine.time.resolution <- function(x, ...) {
  rng <- range(difftime(x[-1L], x[-length(x)], units = "days"))
  a <- round(min(rng))
  if (a == 0) {
    a <- round(24 * as.numeric(min(rng)))
    return(sprintf("%shourly", a))
  }
  if (a == 1) {
    return("daily")
  }
  if (a %in% 5:9) {
    val = "weekly"
  } else {
    val = "monthly"
  }
  val

}



.interp <- function(x1, x2, proportion) {
  x1 * (1 - proportion) + x2 * proportion
}

.calcProportion <- function(xmin, xmax, x) {
  (unclass(x) - unclass(xmin) ) / (unclass(xmax) - unclass(xmin))
}

#' @importFrom raster projection xmin xmax ymin ymax extent
.big.extract <-  function (x, y,
                           ctstime = FALSE,
                           inputfiles = NULL, ...) {
  result <- rep(as.numeric(NA), nrow(y))
  ## progress
  pb <- progress::progress_bar$new(
    format = "getting ready                  [:bar] :percent in :elapsed",
    total = 10, clear = FALSE, width= 80)
  pb$tick(0) ## ---------------------------------------------

  files <- inputfiles
  pb$tick(0) ## ---------------------------------------------
  ## data.frame input has  assumed structure
  ## we assume y is lon,lat,time
  y1 <- sp::SpatialPoints(as.matrix(y[,1:2]),
                      sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
  pb$tick(0) ## ---------------------------------------------
  times <- try(as.POSIXct(y[[3L]], tz = "UTC"))
  y <- y1

  pb$tick(0) ## ---------------------------------------------

  dummy <- x(inputfiles = files)
  yp <- sp::spTransform(y1, raster::projection(dummy))
  pb$tick(0) ## ---------------------------------------------
  xylim <- raster::extent(yp)
  ## expand out a bit for single-location queries
  if (xmax(xylim) == xmin(xylim) | ymax(xylim) == ymin(xylim)) {
    xylim <- xylim + res(dummy)
  }
  dx <- xmax(xylim)-xmin(xylim)
  dy <- ymax(xylim)-ymin(xylim)
  xylim <- xylim + c(dx, dy) / 10
  pb$tick(0) ## ---------------------------------------------

  ## unique indexes
  pb$tick(0) ## ---------------------------------------------

  windex <- findInterval(times, files$date)
  findex <- unique(windex)
  pb$tick(0) ## ---------------------------------------------
  ## this won't always work, need to zap anything out of range . . .
  if (max(times) == max(files$date[findex])) findex <- c(findex, max(findex) + 1)
  findex <- findex[findex <= nrow(files)]
  date <- files$date[findex]
  l <- list(...)
  if ("inputfiles" %in% names(l)) warning("using inputfiles explicitly is deprecated, please don't do it")
  mess1 <- ""
  pb$tick(0) ## ---------------------------------------------
  ## progress

  pb <- progress::progress_bar$new(
    format = "extracting :what file :ith of :nn [:bar] :percent in :elapsed",
    total = length(date), clear = FALSE, width= 80)
  pb$tick(0, tokens = list(ith = 1, nn = length(date)))


  ## interpolate in time?
  if (ctstime) {
    ## we need to store start and end values
    resm <- cbind(result, result)
    thisx1 <- x(date[1L], verbose = FALSE, inputfiles = files, xylim = xylim, ...)  ## inputfiles direct
    for (i in seq_along(date)[-1]) {

      thisx2 <- x(date[i], verbose = FALSE, inputfiles = files, xylim = xylim,  ...)
    ## findInterval is too hard to use reliably
      ## asub <- findInterval(times, date) == (i - 1)
      asub <- windex == findex[i]

      ## interpolation in time, controlled in space by "method" for xy
      if (any(asub)) {resm[asub, ] <- suppressWarnings(extract(stack(thisx1, thisx2), y[asub, ], ...))}
      ## use date since agggregate smashes the Z
      result[asub] <- .interp(resm[asub,1], resm[asub,2], .calcProportion(date[i-1L], date[i], times[asub]))
      ## setup to do the next loop
      thisx1 <- thisx2
      ## report happy times
      if (interactive() & verbose) {
        pb$tick(tokens = list(what = time.resolution, ith = i, nn = length(date)))


      }
    }
    ##                      message("", appendLF = TRUE)
    # if (interactive() & verbose) cat("\n")
  } else {

    ## TODO, fix up the if/else here with an exception for the first/last for ctstime
    for (i in seq_along(date)) {
      thisx <- x(date[i], inputfiles = files, xylim = xylim,  ...)

      asub <- windex == findex[i]
      ## no interpolation in time, controlled by "method" for xy
      if (any(asub)) {result[asub] <- suppressWarnings(extract(thisx, y[asub, ], ...))}

       pb$tick(tokens = list(ith = i, nn = length(date)))


    }
    ##message("", appendLF = TRUE)
    ##if (interactive() & verbose)   cat("\n")

  }
  result

}


##' Extract methods for raster read functions
##'
##' Extract data from read functions in various ways.
##' @title extract
##' @param x A raster read function, with arguments date/inputfiles.
##' @param y Object to use for querying from the rasterread functions, such as a vector of character, Date, or POSIXt values,  data.frame, trip, etc.
##' @param ctstime specify whether to find the nearest value in time (\code{FALSE}), or interpolate between slices (\code{TRUE})
##' @param fact integer. Aggregation factor expressed as number of cells in each direction (horizontally and vertically). Or two integers (horizontal and vertical aggregation factor). See Details in \code{\link[raster]{aggregate}}
##' @param ... Additional arguments passed to the read function.
##' @return data values extracted by the read functions
##' @seealso  \code{\link{extract}}
##' @importFrom raster extract
##' @export
##' @aliases extract,function,data.frame-method
xytract <- function(x, y, ctstime = FALSE, ...) {
  stopifnot(is.function(x))
   if (inherits(x, "trip")) {
     out <- .trip.extract(x, y, ctstime = ctstime, ...)
   } else {
    out <- .big.extract(x, y, ctstime = ctstime, ...)
   }
  out
}

#setMethod("extract", signature(x = 'function', y = 'data.frame'), .big.extract)

longlat_coords <- function(x) {
  x <- as(x, "SpatialPoints")
  if (!raster::couldBeLonLat(x)) {
    x <- sp::spTransform(x, sp::CRS("+init=epsg:4326"))
  }
  as.data.frame(coordinates(x))
}
.trip.extract <- function(x, y, ...) {
  xyt <- longlat_coords(y)
  xyt[["time"]] <- y[[y@TOR.columns[1L]]]
  xytract(x, xyt, ...)
}