R/gwrgrid.R

In fionahevans/agric: Functions for accessing & handling agricultural data.

#' Gaussian kernel weights
#' 
#' Gaussian weights 2d kernel of given radius
#'
#' @param s size.
#'
#' @author Fiona Evans
#' 
#' @return Returns an array with dim=c(s, s).
#' @export
gaussianKernel <- function(radius, sigma=1) {
  x <- c(-radius:radius)
  y <- c(-radius:radius)
  outer(X = x, Y = y, FUN = function(X, Y) return(exp(-(X^2 + Y^2)/(2 * sigma^2))))
}



# Required by gwrgrid.cpp (which is used by gwrgrid)

#' Linear regression 
#' 
#' Linear regression with extra return data
#'
#' @param y vector; Predictor.
#' @param x matrix; Covariate.
#' @param weights vector; Weights.
#'
#' @author Fiona Evans
#' 
#' @return Returns a raster stack.
#' @export
mylm <- function(y, x, weights){
  d <- na.omit(data.frame(x, y, weights))
  ret <- rep(NA, 8)
  if (nrow(d) >= 2 && length(unique(x)) >= 2 && length(unique(y)) >= 2){
    m <- lm(y ~ x, weights=d$weights, data=d)
    ret <- as.vector(summary(m)$coef)
  }
  rn <- c("Intercept", "Slope")
  cn <- c("Estimate", "Std. Error",   "t value",    "Pr(>|t|)")
  names(ret) <- c(outer(rn, cn, FUN=paste, sep=" "))
  ret
}

#' Perform geographically weighted regression on gridded data
#' 
#' Perform geographically weighted regression on gridded data
#'
#' @param y raster; Raster containing response variable.
#' @param x raster; Raster containing covariate.
#' @param radius integer; Radius of GWR window.
#' @param weights matrix; Matrix of weights to use (nrow = 2*radius+1, ncol=2*radius+1).
#'
#' @author Fiona Evans
#' 
#' @return Returns a raster stack.
#' @export
gwrgrid <- function(y, x, radius, weights){

  y_array <- raster::as.array(y)
  x_array <- raster::as.array(x)
  
  # Run to get names for return array
  #z <- lm_R(as.numeric(y_array), as.matrix(as.numeric(x_array)), rep(1, length(y_array)))
  z <- mylm(as.numeric(y_array), as.matrix(as.numeric(x_array)), rep(1, length(y_array)))
  
  
  # Run gwr in C++, result is an array
  res <- gwrgrid_cpp(y_array[,,1], x_array[,,1], radius=radius, weights=weights)
  
  e <- extent(x)
  
  ls <- list()
  for (i in 1:dim(res)[3]) {
    tmp <- res[,,i]
    tmp[tmp==-99999] <- NA
    ls[[i]] = raster(tmp, xmn=e@xmin, xmx=e@xmax, ymn=e@ymin, ymx=e@ymax)
  }
  names(ls) <- names(z)
  
  ret <- raster::stack(ls)
  ret
}
 
subst <- function(Command, ...) do.call("substitute", list(Command, list(...)))


fitCommand <- function(covariates, Var, data = "data", weights=NULL) {
  rhs <- parse(text="1")[[1]]
  if (length(covariates) > 0) rhs <- parse(text=paste(covariates, collapse = "+"))[[1]]
  fitCommand <- Quote(lm(Var ~ rhs, data = d, na.action = na.exclude))
  if (!is.null(weights))  fitCommand <- Quote(lm(Var ~ rhs, data = d, weights=w, na.action = na.exclude))
  fitCommand <- subst(fitCommand, Var = as.name(Var))
  fitCommand <- subst(fitCommand, d = as.name(data))
  if (!is.null(weights)) fitCommand <- subst(fitCommand, w = as.name(weights))
  fitCommand <- subst(fitCommand, rhs = rhs)
  fitCommand
}



#' Perform geographically weighted regression on gridded data
#' 
#' Perform geographically weighted regression on gridded data
#'
#' @param y string; Name of response variable.
#' @param variables string vector; names of covariates/predictors.
#' @param stack raster stack; Data
#' @param radius integer; Radius of GWR window
#'
#' @author Fiona Evans
#' 
#' @return Returns a dataframe.
#' @export
gwrgrid.old <- function(y, variables, stack, radius) {
  # y is a string - name of response variable
  # variables is a vector of names of covariates
  
  n <- radius
  
  nx <- nrow(stack[[y]])
  ny <- ncol(stack[[y]]) 
  
  fw <- gaussianKernel(radius)
  
  out.n <- 4 * (length(variables) + 1)
  
  # buffer by n rows and columns
  inpy <- matrix(nrow=nx+2*n, ncol=ny+2*n)
  inpy[(n+1):(nx+n), (n+1):(ny+n)] <- as.matrix(stack[[y]])
  
  inp <- array(dim=c(nrow=nx+2*n, ncol=ny+2*n, length(variables)))
  for (k in 1:length(variables)){
    inp[(n+1):(nx+n), (n+1):(ny+n), k] <- as.matrix(stack[[variables[k]]])
  }
  
  out <- array(NA, dim=c(nrow(inpx), ncol(inpy), out.n))
  
  for (i in (n+1):(nx+n)){
    for (j in (n+1):(ny+n)) {
      if (!is.na(inpx[i, j])) {
        d <- data.frame(y=as.vector(inpy[(i-n):(i+n),(j-n):(j+n)]))
        
        for (k in 1:length(variables)){
          d[, k+1] <- as.vector(inp[(i-n):(i+n),(j-n):(j+n), k])
        }
        names(d) <- c(y, variables)
        
        dd <<- d
        
        # Make the fitCommand
        fc <- fitCommand(variables, y, data = "d")
        
        # Check the data isn't all NAs
        dofit <- T
        if (all(is.na(dd[,y]))) dofit <- F
        for (k in 1:length(variables)) {
          if (all(is.na(dd[,variables[k]]))) dofit <- F
        }
        
        # Fit the model 
        if (dofit) m <- eval(fc)
        out[i, j, ] <- as.vector(summary(m)$coef)
      }
    }
  }
  
  out <- out[(n+1):(nx+n), (n+1):(ny+n), ]
  
  s <- stack()
  for (k in 1:out.n) {
    r <- stack[[y]]
    r <- setValues(r, out[,, k])
    s <- stack(s, r)
  }
  
  rn <- c("Intercept", variables)
  cn <- c("Estimate", "Std. Error",   "t value",    "Pr(>|t|)")
  
  names(s) <- c(outer(rn, cn, FUN=paste, sep=" "))
  
  s
}