Spatial Analysis in R

Lancaster University

Spatial Statistics

Intro

  • Its not always on a map
  • Lots of packages
  • Lots of techniques
  • Lots of subtleties

Point Patterns

When the coordinate locations are the data.

  • Q: Are these points just 'random'?
  • Q: Is that a cluster?
  • Q: Are these points like those points?
  • Q: Do we get those points where its like this or that?

Geostatistics

When measurements at locations are the data.

  • Q: Where's the gold?
  • Q: How much oil is left?
  • Q: Who has all the rare earth metals?
  • Q: How do I draw a pollution map?

Areal Spatial Regression

Like linear regression models, on areas

  • Q: How are these measures related?
  • Q: Have we missed something in our model?

All the same...

  • A Point Pattern is a Poisson process from a Gaussian field.
  • Geostatistical data are samples from a Gaussian field.
  • Areal measurements are integrals over a Gaussian field...

Point Pattern Analysis

The spatstat package

Point-pattern analysis package

  • Doesn't use sp classes
  • Point randomness tests
  • Point clustering tests
  • Point pattern simulation
  • Point model fitting

spatstat examples

Quadrat counts

data(humberside)
plot(humberside)

##    case control 
##       1       2

qH <- quadratcount(humberside, 8, 8)
plot(qH, add = TRUE, col = "blue", cex = 1.5, lwd = 2)
plot of chunk spatstatex

Point pattern simulation and test

pPois = rpoispp(200)
plot(pPois)
plot of chunk simsim
pSSI <- rSSI(0.05, 200)
plot(pSSI)
plot of chunk simsim
plot(Kest(pPois, correction = "best"), . - theo ~ r)
plot of chunk simsim
##      lty col  key                       label
## iso    1   1  iso hat(K)[iso](r) - K[pois](r)
## theo   2   2 theo     K[pois](r) - K[pois](r)
##                                           meaning
## iso  Ripley isotropic correction estimate of K(r)
## theo                     theoretical Poisson K(r)

plot(Kest(pSSI, correction = "best"), . - theo ~ r)
plot of chunk simsim
##      lty col  key                       label
## iso    1   1  iso hat(K)[iso](r) - K[pois](r)
## theo   2   2 theo     K[pois](r) - K[pois](r)
##                                           meaning
## iso  Ripley isotropic correction estimate of K(r)
## theo                     theoretical Poisson K(r)

The splancs package

  • Very old. Almost prehistoric
  • Mostly superceded by spatstat and other packages
  • One or two things not implemented elsewhere

Splancs space-time analysis

plot of chunk stanalysis

Space-time packages

  • spacetime - data structures for space-time data
  • stpp - space-time point patterns

GeoStatistics

Packages

Automatic interpolation

require(automap)
data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y
kriging_result = autoKrige(zinc ~ soil + ffreq + dist, meuse, meuse.grid)

## [using universal kriging]

plot(kriging_result)
plot of chunk unnamed-chunk-1

Kriging with geoR

geoR uses its own class of objects for its data

require(geoR)
plot(s100)
plot of chunk kriggeor
grid = expand.grid(seq(0, 1, l = 101), seq(0, 1, l = 101))
vario100 <- variog(s100, max.dist = 1)

## variog: computing omnidirectional variogram

ini.vals <- expand.grid(seq(0, 1, l = 5), seq(0, 1, l = 5))
vfit <- variofit(vario100, ini = ini.vals, fix.nug = TRUE, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "1"     "0.25" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 0.167887026209413

plot(vario100)
lines(vfit)
plot of chunk kriggeor
kc <- krige.conv(s100, loc = grid, krige = krige.control(cov.pars = vfit$cov.pars))

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

image(kc)
points(s100, add = TRUE)
plot of chunk kriggeor

Kriging with gstat

gstat was written by one of the sp developers so it uses sp objects everywhere.

data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y
m <- vgm(0.59, "Sph", 874, 0.04)
# ordinary kriging:
x <- krige(log(zinc) ~ 1, meuse, meuse.grid, model = m)

## [using ordinary kriging]

spplot(x["var1.pred"], main = "ordinary kriging predictions")
plot of chunk gstatkrige
spplot(x["var1.var"], main = "ordinary kriging variance")
plot of chunk gstatkrige

Spatial Regression

Packages

  • spdep - spatial dependency
  • spgwr - GWR methodology

spdep functions

Global Moran Test

require(spdep)
require(rgdal)
NY8 = readOGR("./Data/", "NY8_utm18")

## OGR data source with driver: ESRI Shapefile 
## Source: "./Data/", layer: "NY8_utm18"
## with 281 features and 17 fields
## Feature type: wkbPolygon with 2 dimensions

nb = poly2nb(NY8)
moran.test(NY8$Cases, listw = nb2listw(nb))

## 
## 	Moran's I test under randomisation
## 
## data:  NY8$Cases  
## weights: nb2listw(nb)  
##  
## Moran I statistic standard deviate = 4.393, p-value = 5.594e-06
## alternative hypothesis: greater 
## sample estimates:
## Moran I statistic       Expectation          Variance 
##          0.156020         -0.003571          0.001320

moran.plot(NY8$Cases, listw = nb2listw(nb, style = "C"))
plot of chunk nymoranspc

Local Moran Statistic

localM = localmoran(NY8$Cases, listw = nb2listw(nb))
NY8$moranZ = localM[, "Z.Ii"]
spplot(NY8, "moranZ")
plot of chunk locmoran
as.character(NY8[NY8$moranZ > 8, ]$AREANAME)

## [1] "Vestal town" "Auburn city" "Auburn city"

spautolm

m <- glm(PROPCAS ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8)
sm <- spautolm(PROPCAS ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, listw = nb2listw(nb), 
    family = "SAR", method = "full")
summary(m)$coef

##               Estimate Std. Error t value  Pr(>|t|)
## (Intercept)  4.129e-04  1.524e-04   2.709 0.0071671
## PEXPOSURE    7.112e-05  3.371e-05   2.110 0.0357572
## PCTAGE65P    2.127e-03  5.821e-04   3.654 0.0003084
## PCTOWNHOME  -4.062e-04  1.637e-04  -2.481 0.0136952

summary(sm)$Coef

##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  4.036e-04  0.0001371   2.944 3.241e-03
## PEXPOSURE    6.632e-05  0.0000285   2.327 1.997e-02
## PCTAGE65P    2.445e-03  0.0005466   4.473 7.712e-06
## PCTOWNHOME  -4.500e-04  0.0001454  -3.096 1.962e-03

Geographically Weighted Regression

spgwr

require(spgwr)

## Loading required package: spgwr

## NOTE: This package does not constitute approval of GWR as a method of
## spatial analysis