R/utilities.R

In objornstad/ncf: Spatial Covariance Functions

#' @title Simple wrapper around symbols to visualize spatial data
#' @description \code{spatial.plot} is a quick function to visualize spatial data using bubble plots
#' @param x vector of length n representing the x coordinates.
#' @param y vector of length n representing the y coordinates.
#' @param z vector of n representing the observation at each location.
#' @param ctr If TRUE, observations will be centered before plotting (zero-sized symbols represents average observations); if FALSE, the original observations are used.
#' @param add If TRUE, a lisa-plot will be added to a pre-existing plot.
#' @param inches scales the size of the symbols
#' @param \dots other arguments
#' @return A bubble-plot of the spatial data is produced.
#' @details This is a simple function to visualize spatial data. Positive (or above average) observations are shown by red circles, Negative (or below average) observations are shown as black squares. For hot/coldspot analysis using Local indicators of spatial association use \code{\link{lisa}}.
#' @references Ripley, B.D. (1987). Stochastic Simulation. Wiley.
#' @author Ottar N. Bjornstad \email{onb1@psu.edu}
#' @seealso \code{\link{lisa}}
#' @examples
#' # first generate some sample data
#' x <- expand.grid(1:20, 1:5)[, 1]
#' y <- expand.grid(1:20, 1:5)[, 2]
#' 
#' # z data from an exponential random field
#' z <- rmvn.spa(x = x, y = y, p = 2, method = "gaus")
#' 
#' # plot data
#' \dontrun{spatial.plot(x = x, y = y, z = z, ctr = FALSE)}
#' @keywords spatial
#' @export
################################################################################
spatial.plot <- function(x, y, z, ctr = TRUE, add = FALSE, inches = 0.2, ...) {
  ##############################################################################
  if (ctr) {
    z <- z - mean(z, na.rm = TRUE)
  }
  
  if (add == FALSE) {
    plot(x, y, type = "n", ...)
  }
  sel <- is.finite(z)
  x <- split(x, z > 0)
  y <- split(y, z > 0)
  sel <- split(sel, z > 0)
  z2 <- split(z, z > 0)
  
  if (!is.null(length(z2[[1]][sel[[1]]]))) {
    symbols(x[[1]][sel[[1]]], y[[1]][sel[[1]]], squares = -z2[[1]][sel[[1]]], 
            inches = inches, add = TRUE, fg = 1, bg = 1)
  }
  
  if (!is.null(length(z2[[1]][sel[[2]]]))) {
    symbols(x[[2]][sel[[2]]], y[[2]][sel[[2]]], circles = z2[[2]][sel[[2]]], 
            inches = inches, add = TRUE, fg = 2, bg = 2)}
}

#' @title Simulate spatially correlated data
#' @description Function to generate spatially autocorrelated random normal variates using the eigendecomposition method. Spatial covariance can follow either and exponential or Gaussian model.
#' @param x vector of length n representing the x coordinates (or latitude; see latlon).
#' @param y vector of length n representing the y coordinates (or longitude).
#' @param p the range of the spatial models.
#' @param method correlation function "exp" (exponential) or "gaus" (gaussian). Exponential is the default.
#' @param nugget correlation at the origin (defaults to one)
#' @return A vector of spatially correlated random normal variates with zero mean and unit variance is returned
#' @details A target covariance matrix A between the n units is generated by calculating the distances between the locations and thereafter evaluating the covariance function in each pairwise distance. A vector, Z, of spatially correlated normal data with the target covariance is subsequently generated using the eigendecomposition method (Ripley, 1987).
#' @references Ripley, B.D. (1987). Stochastic Simulation. Wiley.
#' @author Ottar N. Bjornstad \email{onb1@psu.edu}
#' @seealso \code{\link{mSynch}}
#' @keywords smooth regression
#' @export
################################################################################
rmvn.spa <- function(x, y, p, method = "exp", nugget = 1) {
  ##############################################################################
  # Function to generate spatially autocorrelated random normal variates using the 
  # eigendecomposition method. Spatial covariance can follow either and exponential 
  # or Gaussian model. 
  ##############################################################################
  imeth <- charmatch(method, c("exp", "gaus"), nomatch = NA)
  if (is.na(imeth)) stop("method should be \"exp\", or \"gaus\"")
  
  ch.fn <- function(n, mu, vmat, tol = 1e-007) {
    p <- ncol(vmat)
    vs <- svd(vmat)
    vsqrt <- t(vs$v %*% (t(vs$u) * sqrt(vs$d)))
    ans <- matrix(rnorm(n * p), nrow = n) %*% vsqrt
    ans <- sweep(ans, 2, mu, "+")
    dimnames(ans) <- list(NULL, dimnames(vmat)[[2]])
    ans
  }
  
  n <- length(x)
  z <- matrix(0, ncol = 1, nrow = n)
  tmpmat <- cbind(x + 0, y + 0)
  xy <- tmpmat[, 1:2] + 0
  dmat <- sqrt(outer(x, x, "-")^2 + outer(y, y, "-")^2)
  if (imeth == 1) {
    covmat <- nugget * exp(-dmat/p)
  } else {
    if (imeth == 2) {
      covmat <- nugget * exp(-(dmat/p)^2)
    }
  }
  z <- ch.fn(1, mu = rep(0, n), vmat = covmat)
  t(z)[, 1]
}

#' @title Utility function
#' @description Called by various functions to calculate various intercepts.
#' @param u a vector.
#' @param v a vector.
#' @param w a vector.
#' @param moran a matrix.
#' @param df a scalar.
#' @param xpoints a vector.
#' @param filter a logical.
#' @param fw a scalar
#' @return A list is returned.
#' @details An auxiliary function to ease  maintenance.
#' @keywords misc
#' @export
gather <- function(u, v, w, moran, df, xpoints, filter, fw) {
  cbar <- mean(v, na.rm = TRUE)
  sobj <- smooth.spline(u, v, df = df)
  x <- xpoints
  y <- predict(sobj, x = x)$y
  if (filter == TRUE) {
    if (fw > 0) {
      y[x > fw] <- 0
    }
    y <- ff.filter(y)
  }
  #if(is.null(w)){
  yint <- y[1]
  #	}
  #else {
  #	 yint <- mean(diag(moran), na.rm= TRUE)
  #	# yint <- ifelse(is.finite(mean(diag(moran), na.rm= TRUE)), mean(diag(moran), na.rm= TRUE),  y[1])
  #	}
  x2 <- x[x >= 0]
  y2 <- y[x >= 0]
  konst <- ifelse(y2[1] > 0, 1, -1)
  xint <- eint <- cint <- NA
  wh <- (1:length(x2))[y2 < 0][1]
  int <- c(x2[wh], x2[wh - 1])
  if (length(int) == 2 & all(is.finite(int))) {
    #xi=nleqslv(0, fn=function(x) predict(sobj, x)$y)$x
    xi <- uniroot(f = function(x) predict(sobj, x)$y, interval = int)$root
    xint <- c(0, xi*konst, NA)[findInterval(xi, range(xpoints)) + 1]
  }
  wh <- (1:length(x2))[(y2 - 1/exp(1)) < 0][1]
  int <- c(x2[wh], x2[wh - 1])
  if (length(int) == 2 & all(is.finite(int))) {
    ei <- uniroot(f = function(x) { predict(sobj, x)$y - 1/exp(1) }, interval = int)$root
    eint <- c(0, ei, NA)[findInterval(ei, range(xpoints)) + 1]
  }
  wh <- (1:length(x2))[(y2 - cbar) < 0][1]
  int <- c(x2[wh], x2[wh - 1])
  if (length(int) == 2 & all(is.finite(int))) {
    ci <- uniroot(f = function(x) { predict(sobj, x)$y - cbar }, interval = int)$root
    cint <- c(0, ci, NA)[findInterval(ci, range(xpoints)) + 1]
  }
  list(x = x, y = y, yint = yint, cbar = cbar, xint = xint, eint = eint, cint = cint)
}

#' @title Great-circle distance
#' @description Great-circle distance function to calculate spatial distance from lat-long data.  Called by various functions.
#' @param x vector of longitudes.
#' @param y vector of latitudes.
#' @return The distance in km is returned
#' @keywords misc
#' @export
################################################################################
gcdist <- function(x, y) {
  # vecotorized gcdist function
  ##############################################################################
  r <- 360/(2 * pi)
  lon <- x/r
  lat <- y/r
  dlon <- outer(lon, lon, "-")
  dlat <- outer(lat,lat,"-")
  a <- (sin(dlat/2))^2 + outer(cos(lat), cos(lat), "*") * (sin(dlon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1 - a))
  return(6370 * c)
}

#' @title Circular correlation
#' @description A vectorized function to calculate a correlation matrix for panels of data.
#' @param x a matrix.
#' @param y an optional second matrix.
#' @return A correlation matrix is returned.
#' @details Missing values are not allowed.
#' @references Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 8.2, World Scientific Press, Singapore.
#' @keywords misc
#' @export
################################################################################
circ.cor2 <- function(x, y = NULL) {
  ##############################################################################
  # Fast vectorized circular correlation 
  ##############################################################################
  circ.mean <- function(x) {
    x <- na.omit(x)
    sinr <- sum(sin(x))
    cosr <- sum(cos(x))
    circmean <- atan2(sinr, cosr)
    circmean
  }
  
  if (is.vector(x) & is.vector(y)) {
    x <- as.matrix(x)
    y <- as.matrix(y)
  }
  x.bar <- apply(x, 2, circ.mean)
  y.bar <- apply(y, 2, circ.mean)
  num <- tcrossprod(sin(t(x) - x.bar), sin(t(y) - y.bar))
  den <- sqrt(outer(
    apply(t(sin(t(x) - x.bar))^2, 2, sum),
    apply(t(sin(t(y) - y.bar))^2, 2, sum), "*"))
  r <- num/den
  return(r)
}

#' @title Utility function
#' @description Called by various functions to calculate Pearson or angular correlation matrices.
#' @param x a matrix.
#' @param y an optional second matrix.
#' @param circ If TRUE, the observations are assumed to be angular (in radians), and circular correlation is used. If FALSE, Pearson product moment correlations is returned.
#' @return A correlation matrix is returned.
#' @details An auxilliary function to ease the maintenance.
#' @references Jammalamadaka, S. Rao and SenGupta, A. (2001). Topics in Circular Statistics, Section 8.2, World Scientific Press, Singapore.
#' @keywords misc
#' @export
################################################################################
cor2 <- function(x, y = NULL, circ = FALSE) {
  ##############################################################################
  circ.cor <- function(x, y) {
    circ.mean <- function(x) {
      sinr <- sum(sin(x))
      cosr <- sum(cos(x))
      circmean <- atan2(sinr, cosr)
      circmean
    }
    ok <- (is.finite(x) & is.finite(y))
    x <- x[ok]
    y <- y[ok]
    n <- length(x)
    r <- NA
    if (n >= 2) {
      x.bar <- circ.mean(x)
      y.bar <- circ.mean(y)
      num <- sum(sin(x - x.bar) * sin(y - y.bar))
      den <- sqrt(sum(sin(x - x.bar)^2) * sum(sin(y - y.bar)^2))
      r <- num/den
    }
    return(r)
  }
  
  if (is.data.frame(x))
    x <- as.matrix(x)
  if (is.data.frame(y))
    y <- as.matrix(y)
  if (!is.matrix(x) && is.null(y))
    stop("supply both x and y or a matrix-like x")
  if (!circ) {
    cor <- cor(x = x, y = y, use = "pairwise.complete.obs", method = "pearson")
  }
  
  if (circ) {
    if (is.null(y)) {
      y <- x
    }
    
    if (all(is.finite(x) & is.finite(y))) {
      cor <- circ.cor2(x, y)
    }
    
    if (!all(is.finite(x) & is.finite(y))) {	
      cat("Missing values with circular data: calculations may take some time.")
      m <- dim(x)[2]
      # this part (for circular correlations is likely to be SLOW (it's a double loop)
      # I should use distm() from geosphere instead
      cor <- matrix(NA, nrow = m, ncol = m)
      for (i in 1:m) {
        for (j in i:m) {
          tmp <- circ.cor(as.vector(x[, i]), as.vector(y[, j]))
          cor[j, i] <- tmp
          cor[i, j] <- tmp
        }
      }
    }
  }
  return(cor)
}

#' @title Fourier filter for correlation functions.
#' @description Fourier filter to ensure positive semi-definite correlation functions. Called by various functions.
#' @param x a vector.
#' @return A vector is returned whose Fourier-transform has no non-negative coefficients.
#' @seealso \code{\link{Sncf}}
#' @keywords misc
#' @export
################################################################################
ff.filter <- function(x) {
  ##############################################################################
  # Fourier filter function
  ##############################################################################
  n <- length(x)
  x2 <- c(x, rev(x[-c(1, n)]))
  fo <- Re(fft(x2))
  fo[fo < 0] <- 0
  ny <- Re(fft(fo, inverse = TRUE)/length(x2))
  return(ny[1:n])
}