8. Local statistics

8.1 Introduction

This handout accompanies Chapter 8 in O’Sullivan and Unwin (2010).

8.2 LISA

Get the data

Download the data for this lesson here: download here: auctb, counties, precipitation. Or with the script below.

dir.create('data', showWarnings = FALSE)
files <- c('auctb.rds', 'counties.rds', 'precipitation.csv')
for (filename in files) {
    localfile <- file.path('data', filename)
    if (!file.exists(localfile)) {
        download.file(paste0('https://biogeo.ucdavis.edu/data/rspatial/', filename), dest=localfile)
    }
}

Read the Auckland data.

library(raster)
auck <- readRDS('data/auctb.rds')
## Warning in readRDS("data/auctb.rds"): invalid or incomplete compressed data
## Error in readRDS("data/auctb.rds"): error reading from connection

Local Getis Gi

library(spdep)
wr <- poly2nb(auck, row.names=auck$Id, queen=FALSE)
## Error in extends(class(pl), "SpatialPolygons"): object 'auck' not found
lstw <- nb2listw(wr, style='B')
## Error in nb2listw(wr, style = "B"): object 'wr' not found
Gi <- localG(auck$TB, lstw)
## Error in localG(auck$TB, lstw): object 'lstw' not found
head(Gi)
## Error in head(Gi): object 'Gi' not found
par(mai=c(0,0,0,0))
Gcuts <- cut(Gi, 5)
## Error in cut(Gi, 5): object 'Gi' not found
Gcutsi <- as.integer(Gcuts)
## Error in eval(expr, envir, enclos): object 'Gcuts' not found
cols <- rev(gray(seq(0,1,.2)))
plot(auck, col=cols[Gcutsi])
## Error in plot(auck, col = cols[Gcutsi]): object 'auck' not found
legend('bottomleft', levels(Gcuts), fill=cols)
## Error in levels(Gcuts): object 'Gcuts' not found

To get the Gi * we need to include each polygon as its own ‘neighbor’

ws <- include.self(wr)
## Error in include.self(wr): object 'wr' not found
lstws <- nb2listw(ws, style='B')
## Error in nb2listw(ws, style = "B"): object 'ws' not found
Gis <- localG(auck$TB, lstws)
## Error in localG(auck$TB, lstws): object 'lstws' not found
Gscuts <- cut(Gis, 5)
## Error in cut(Gis, 5): object 'Gis' not found
Gscutsi <- as.integer(Gscuts)
## Error in eval(expr, envir, enclos): object 'Gscuts' not found
cols <- rev(gray(seq(0,1,.2)))
plot(auck, col=cols[Gscutsi])
## Error in plot(auck, col = cols[Gscutsi]): object 'auck' not found
legend('bottomleft', levels(Gscuts), fill=cols)
## Error in levels(Gscuts): object 'Gscuts' not found

This looks very similar to the local average.

m <- sapply(ws, function(i) mean(auck$TB[i]))
## Error in lapply(X = X, FUN = FUN, ...): object 'ws' not found
cts <- cut(m, 5)
## Error in cut.default(m, 5): 'x' must be numeric
mcts <- as.integer(cts)
## Error in eval(expr, envir, enclos): object 'cts' not found
plot(auck, col=cols[mcts])
## Error in plot(auck, col = cols[mcts]): object 'auck' not found
legend('bottomleft', levels(cts), fill=cols)
## Error in levels(cts): object 'cts' not found

The local Moran Ii shows where there are locally high levels of autocorrelation.

Ii <- localmoran(auck$TB, lstw)
## Error in is.vector(x): object 'auck' not found
Icuts <- cut(Ii, 5)
## Error in cut(Ii, 5): object 'Ii' not found
Icutsi <- as.integer(Icuts)
## Error in eval(expr, envir, enclos): object 'Icuts' not found
plot(auck, col=cols[Icutsi])
## Error in plot(auck, col = cols[Icutsi]): object 'auck' not found
legend('bottomleft', levels(Icuts), fill=cols)
## Error in levels(Icuts): object 'Icuts' not found

8.3 Geographically weighted regression

Here is an example of GWR with California precipitation data (you can download the data with the scripts or links on the top of this page).

cts <- readRDS('data/counties.rds')
## Warning in readRDS("data/counties.rds"): invalid or incomplete compressed
## data
## Error in readRDS("data/counties.rds"): error reading from connection
p <- read.csv('data/precipitation.csv')
head(p)
##      ID                 NAME   LAT    LONG ALT  JAN FEB MAR APR MAY JUN
## 1 ID741         DEATH VALLEY 36.47 -116.87 -59  7.4 9.5 7.5 3.4 1.7 1.0
## 2 ID743  THERMAL/FAA AIRPORT 33.63 -116.17 -34  9.2 6.9 7.9 1.8 1.6 0.4
## 3 ID744          BRAWLEY 2SW 32.96 -115.55 -31 11.3 8.3 7.6 2.0 0.8 0.1
## 4 ID753 IMPERIAL/FAA AIRPORT 32.83 -115.57 -18 10.6 7.0 6.1 2.5 0.2 0.0
## 5 ID754               NILAND 33.28 -115.51 -18  9.0 8.0 9.0 3.0 0.0 1.0
## 6 ID758        EL CENTRO/NAF 32.82 -115.67 -13  9.8 1.6 3.7 3.0 0.4 0.0
##   JUL  AUG SEP OCT NOV DEC
## 1 3.7  2.8 4.3 2.2 4.7 3.9
## 2 1.9  3.4 5.3 2.0 6.3 5.5
## 3 1.9  9.2 6.5 5.0 4.8 9.7
## 4 2.4  2.6 8.3 5.4 7.7 7.3
## 5 8.0  9.0 7.0 8.0 7.0 9.0
## 6 3.0 10.8 0.2 0.0 3.3 1.4

plot(cts)
## Error in plot(cts): object 'cts' not found
points(p[,c('LONG', 'LAT')], col='red', pch=20)
## Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

Compute annual average precipitation

p$pan <- rowSums(p[,6:17])

Global regression model

m <- lm(pan ~ ALT, data=p)
m
##
## Call:
## lm(formula = pan ~ ALT, data = p)
##
## Coefficients:
## (Intercept)          ALT
##      523.60         0.17

Create Spatial* objects with a planar crs.

alb <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
sp <- p
coordinates(sp) = ~ LONG + LAT
crs(sp) <- "+proj=longlat +datum=NAD83"
spt <- spTransform(sp, alb)
ctst <- spTransform(cts, alb)
## Error in spTransform(cts, alb): object 'cts' not found

Get the optimal bandwidth

library( spgwr )
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
bw <- gwr.sel(pan ~ ALT, data=spt)
## Bandwidth: 526221.1 CV score: 64886883
## Bandwidth: 850593.6 CV score: 74209073
## Bandwidth: 325747.9 CV score: 54001118
## Bandwidth: 201848.6 CV score: 44611213
## Bandwidth: 125274.7 CV score: 35746320
## Bandwidth: 77949.39 CV score: 29181737
## Bandwidth: 48700.74 CV score: 22737197
## Bandwidth: 30624.09 CV score: 17457161
## Bandwidth: 19452.1 CV score: 15163436
## Bandwidth: 12547.43 CV score: 19452191
## Bandwidth: 22792.75 CV score: 15512988
## Bandwidth: 17052.67 CV score: 15709960
## Bandwidth: 20218.99 CV score: 15167438
## Bandwidth: 19767.99 CV score: 15156913
## Bandwidth: 19790.05 CV score: 15156906
## Bandwidth: 19781.39 CV score: 15156902
## Bandwidth: 19781.48 CV score: 15156902
## Bandwidth: 19781.47 CV score: 15156902
## Bandwidth: 19781.47 CV score: 15156902
## Bandwidth: 19781.47 CV score: 15156902
## Bandwidth: 19781.47 CV score: 15156902
bw
## [1] 19781.47

Create a regular set of points to estimate parameters for.

r <- raster(ctst, res=10000)
## Error in raster(ctst, res = 10000): object 'ctst' not found
r <- rasterize(ctst, r)
## Error in rasterize(ctst, r): object 'ctst' not found
newpts <- rasterToPoints(r)
## Error in nlayers(x): object 'r' not found

Run the gwr function

g <- gwr(pan ~ ALT, data=spt, bandwidth=bw, fit.points=newpts[, 1:2])
## Error in is(fit.points, "Spatial"): object 'newpts' not found
g
## Error in eval(expr, envir, enclos): object 'g' not found

Link the results back to the raster

slope <- r
## Error in eval(expr, envir, enclos): object 'r' not found
intercept <- r
## Error in eval(expr, envir, enclos): object 'r' not found
slope[!is.na(slope)] <- g$SDF$ALT
## Error in eval(expr, envir, enclos): object 'g' not found
intercept[!is.na(intercept)] <- g$SDF$'(Intercept)'
## Error in eval(expr, envir, enclos): object 'g' not found
s <- stack(intercept, slope)
## Error in stack(intercept, slope): object 'intercept' not found
names(s) <- c('intercept', 'slope')
## Error in names(s) <- c("intercept", "slope"): object 's' not found
plot(s)
## Error in plot(s): object 's' not found

See this page for a more detailed example of Geographically weighted regression.