R/peExtract.R

In ssaxe-usgs/parallelExtract:

library(raster)
library(rgdal)
library(dplyr)
# polygon
y <- "C:/Users/ssaxe/Documents/Scripts/R Scripts/Model Evaluation Tool/shapefiles/arcSimplify/HUC10/HUC_resolved2.shp" %>%
  rgdal::readOGR(.,
                 stringsAsFactors = F,
                 verbose = F)
# raster and reproject
x <- "C:/Users/ssaxe/Documents/Scripts/R Scripts/Model Evaluation Tool/raster/Storage/Grace-JPL-Unprojected/GRACE-JPL.tif" %>%
  raster::raster(.) %>%
  raster::rotate(.) %>%
  raster::projectRaster(from = .,
                        crs = raster::crs(raster::projection(y))) %>%
  raster::crop(., raster::extent(y), snap = 'out')




peExtract <- function(x,
                      y,
                      fun=NULL,
                      na.rm=FALSE,
                      weights=FALSE,
                      normalizeWeights=TRUE,
                      cellnumbers=FALSE,
                      small=TRUE,
                      df=FALSE,
                      layer){
  # libraries
  require(raster)

  # Compare projections and reproject if necessary
  px <- raster::projection(x, asText=FALSE)
  comp <- raster::compareCRS(px, raster::projection(y), unknown=TRUE)
  if (!comp) {
    stop('Raster and SpatialPolygons have conflicting projections.')
  }

  # get bounding boxes for each input
  spbb <- bbox(y)
  rsbb <- bbox(x)

  # get res of raster
  addres <- max(res(x))

  # get number of features in shapefile (npol = 15955)
  npol <- length(y@polygons)

  # dealing with user-supplied function
  if (!is.null(fun)) {
    cellnumbers <- FALSE
    if (weights) {
      if (!is.null(fun)) {
        test <- try(methods::slot(fun, 'generic') == 'mean', silent=TRUE)
        if (!isTRUE(test)) {
          warning('"fun" was changed to "mean"; other functions cannot be used when "weights=TRUE"' )
        }
      }
      fun <- function(x, ...) {
        # some complexity here because different layers could
        # have different NA cells
        if ( is.null(x) ) {
          return(rep(NA, nl))
        }
        w <- x[,nl+1]
        x <- x[,-(nl+1), drop=FALSE]
        x <- x * w
        w <- matrix(rep(w, nl), ncol=nl)
        w[is.na(x)] <- NA
        w <- colSums(w, na.rm=TRUE)
        x <- apply(x, 1, function(X) { X / w } )
        if (!is.null(dim(x))) {
          rowSums(x, na.rm=na.rm)
        } else {
          sum(x, na.rm=na.rm)
        }
      }
    }

    if (sp) {
      df <- TRUE
    }

    doFun <- TRUE

  } else {
    if (sp) {
      sp <- FALSE
      df <- FALSE
      warning('argument sp=TRUE is ignored if fun=NULL')
      #} else if (df) {
      #	df <- FALSE
      #	warning('argument df=TRUE is ignored if fun=NULL')
    }

    doFun <- FALSE
  }

  # identify correct layer to use
  if (missing(layer)) {
    layer <- 1
  } else {
    layer <- max(min(nlayers(x), layer), 1)
  }

  # identify number of layers to extract from (hardcode, one layer only.  If there is more than one layer in
  # the raster, we will extract cell numbers)
  nl <- 1


  if (spbb[1,1] >= rsbb[1,2] | spbb[1,2] <= rsbb[1,1] | spbb[2,1] >= rsbb[2,2] | spbb[2,2] <= rsbb[2,1]) {
    return(NULL)
  }


  # make sure raster is a single layer
  rr <- raster(x)
  rr2 <- x

  # Extract



}