R codes used in Spatial Relationships Between Two Georeferenced Variables

Below you will find links to R scripts and codes used in the book; please inform us if you notice anything missing.

Codispersion map for non-rectangular grids

The codispersion coefficient quantifies the association between two spatial processes for a particular direction (spatial lag) on the two-dimensional space. Following R code reflects the computational procedure to estimate the values of the codispersion and the plotting procedure to yield the codispersion map. This code also is available as an R source file: codismap.R 

codisp.map <- function(x, y, coords, nclass = 13, ncell = 40)
{
  require("akima")
  require("fields")
  require("geoR")

  rho <- function(x, y, uvec, max.dist, angle) {
    z  <- as.geodata(cbind(x$coords, x$data + y$data))
    nz <- variog(z, uvec = uvec, max.dist = max.dist, direction = angle, messages = FALSE)
    dx <- variog(x, uvec = uvec, max.dist = max.dist, direction = angle, messages = FALSE)
    dy <- variog(y, uvec = uvec, max.dist = max.dist, direction = angle, messages = FALSE)
    cf <- .5 * (nz$v - dx$v - dy$v) / sqrt(dx$v * dy$v)
    cf
  }

  x <- as.geodata(cbind(coords, x))
  y <- as.geodata(cbind(coords, y))
  dmax <- .5 * max(dist(coords))
  angles <- seq(from = 0, to = pi, by = 0.01)
  nangles <- length(angles)
  uvec <- seq(from = 0, to = dmax, length = nclass + 1)[-1]

  xcirc <- ycirc <- 0

  for (i in seq_len(nclass)) {
    xcirc[(nangles*(i-1)+1):(nangles*i)] <- seq(-uvec[i], uvec[i], length = nangles)
    ycirc[(nangles*(i-1)+1):(nangles*i)] <- sqrt(uvec[i]^2 - xcirc[(nangles*(i-1)+1):(nangles*i)]^2)
  }

  z <- matrix(0, nrow = nangles, ncol = nclass)
  for (i in seq_len(nangles))
    z[i,] <- rho(x, y, uvec = uvec, max.dist = dmax, angle = angles[i])
  z <- as.vector(z)

  par(pty = "s")
  xl <- seq(min(xcirc), max(ycirc), length = ncell)
  yl <- seq(min(ycirc), max(ycirc), length = ncell)
  image.plot(interp(xcirc, ycirc, z, xo = xl,yo = yl), col = topo.colors(256), xlab = "x", ylab = "y")
  title(main = "Codispersion map")
}

Source:
Vallejos, R., Osorio, F., Mancilla, D. (2015). The codispersion map: a graphical tool to visualize the association between two spatial variables. Statistica Neerlandica 69, 298-314.

Codispersion map for rectangular grids

The R function hcodisp.map provides a graphical tool called codispersion map to visualize the spatial correlation between two sequences measured for regular lattices on the plane. Sources for this codispersion map are available at the next link:

To use this funtion, first download the source as a tarball (.tar.gz) file. Unpack the tarball (thereby creating a directory named, hcodisp.map) and to create the Dynamically Loaded (DL) library, enter at the console prompt: 

R CMD SHLIB -o hcodisp.so *.c

To use this funtion, start R and enter: 

source("hcodisp.map.R")
dyn.load("hcodisp.so")
z <- hcodisp.map(x, y)

Source:
Vallejos, R., Osorio, F., Mancilla, D. (2015). The codispersion map: a graphical tool to visualize the association between two spatial variables. Statistica Neerlandica 69, 298-314.

Codispersion map using kernel estimation

Codispersion map between two spatial processes using a nonparametric approach. An R script and the source for the computatioon of the codispersion coefficient are available at: codisp.kern_script.R and codisp.kern.R , respectively. 

## Script to compute the codispersion map using kernel estimation

library(fields) # required to see the output image

# Warning: this operation could take a long time
lenx <- length(x)
leny <- length(y)
Mcodisp <- matrix(0, ncol = leny, nrow = lenx)

for (i in 1:lenx) {
  for (j in 1:leny) {
    Mcodisp[i,j] <- codisp.kern(x, y, c(x[i], y[j]), h)[4]
  }
}

image.plot(x, y, Mcodisp, col = topo.colors(128), ylab = "", xlab = "", main = "Codispersion map")

Function to compute the codispersion coefficient between two spatial processes. The routine uses geodata objects (from geoR package): 

codisp.kern <- function(x, y, k, h, gamma = 1)
{
  kernel <- function(u, gamma) {
    v = 0
    v = ifelse(abs(u) <= 1, (1 / beta(0.5, gamma + 1)) * (1 - u^2)^gamma, 0)
  }

  ifelse(x$coords == y$coords, 1, stop("The coordinates of x and y are different"))

  n <- length(x$data)
  xx <- matrix(0, ncol = n, nrow = n)
  yy <- matrix(0, ncol = n, nrow = n)
  xy <- matrix(0, ncol = n, nrow = n)
  zx <- matrix(0, ncol = n, nrow = n)
  zy <- matrix(0, ncol = n, nrow = n)

  for (i in 1:n) {
    for (j in i:n) {
      xx[i,j] <- (x$data[i] - x$data[j])^2
      yy[i,j] <- (y$data[i] - y$data[j])^2
      xy[i,j] <- (x$data[i] - x$data[j]) * (y$data[i] - y$data[j])
      zx[i,j] <- x$coords[i,1] - x$coords[j,1]
      zy[i,j] <- x$coords[i,2] - x$coords[j,2]
    }
  }

  zx <- zx - t(zx)
  zy <- zy - t(zy)
  xx <- xx + t(xx)
  yy <- yy + t(yy)
  xy <- xy + t(xy)

  kern.xx <- kernel((k - zx) / h[1], gamma) * kernel((k - zy) / h[1], gamma)
  kern.yy <- kernel((k - zx) / h[2], gamma) * kernel((k - zy) / h[2], gamma)
  kern.xy <- kernel((k - zx) / h[3], gamma) * kernel((k - zy) / h[3], gamma)

  num  <- sum(kern.xy * xy) / (2 * sum(kern.xy))
  den1 <- sum(kern.xx * xx) / (2 * sum(kern.xx))
  den2 <- sum(kern.yy * yy) / (2 * sum(kern.yy))

  o1 <- den1
  o2 <- den2
  o3 <- num
  o4 <- num / sqrt(den1 * den2)
  c(o1, o2, o3, o4)
}

Source:
Cuevas, F., Porcu, E., Vallejos, R. (2013). Study of spatial relationships between two sets of variables: A nonparametric approach. Journal of Nonparametric Statistics 25, 695-714.