# Local statistics¶

## Introduction¶

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

## LISA¶

### Get the data¶

if (!require("rspatial")) devtools::install_github('rspatial/rspatial')
library(rspatial)
auck <- sp_data('auctb.rds')


Local Getis Gi

library(spdep)
wr <- poly2nb(auck, row.names=auck$Id, queen=FALSE) lstw <- nb2listw(wr, style='B') Gi <- localG(auck$TB, lstw)
## [1]  0.4241759  0.5623307  0.4047413  0.6819762 -1.3278352 -1.4086435

par(mai=c(0,0,0,0))
Gcuts <- cut(Gi, 5)
Gcutsi <- as.integer(Gcuts)
cols <- rev(gray(seq(0,1,.2)))
plot(auck, col=cols[Gcutsi])
legend('bottomleft', levels(Gcuts), fill=cols)


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

ws <- include.self(wr)
lstws <- nb2listw(ws, style='B')
Gis <- localG(auck$TB, lstws) Gscuts <- cut(Gis, 5) Gscutsi <- as.integer(Gscuts) cols <- rev(gray(seq(0,1,.2))) plot(auck, col=cols[Gscutsi]) legend('bottomleft', levels(Gscuts), fill=cols)  This looks very similar to the local average. m <- sapply(ws, function(i) mean(auck$TB[i]))
cts <- cut(m, 5)
mcts <- as.integer(cts)
plot(auck, col=cols[mcts])
legend('bottomleft', levels(cts), fill=cols)


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

Ii <- localmoran(auck$TB, lstw) Icuts <- cut(Ii, 5) Icutsi <- as.integer(Icuts) plot(auck, col=cols[Icutsi]) legend('bottomleft', levels(Icuts), fill=cols)  ## 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 <- sp_data('counties.rds') p <- sp_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) points(p[,c('LONG', 'LAT')], col='red', pch=20)  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
spt <- spTransform(sp, alb)
ctst <- spTransform(cts, alb)


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)
r <- rasterize(ctst, r)
newpts <- rasterToPoints(r)


Run the gwr function

g <- gwr(pan ~ ALT, data=spt, bandwidth=bw, fit.points=newpts[, 1:2])
g
## Call:
## gwr(formula = pan ~ ALT, data = spt, bandwidth = bw, fit.points = newpts[,
##     1:2])
## Kernel function: gwr.Gauss
## Fixed bandwidth: 19781.47
## Fit points: 4087
## Summary of GWR coefficient estimates at fit points:
##                    Min.    1st Qu.     Median    3rd Qu.      Max.
## X.Intercept. -702.40121   79.54254  330.48807  735.42718 3468.8702
## ALT            -3.91270    0.03058    0.20461    0.41542    4.6133


Link the results back to the raster

slope <- r
intercept <- r
slope[!is.na(slope)] <- g$SDF$ALT
intercept[!is.na(intercept)] <- g$SDF$'(Intercept)'
s <- stack(intercept, slope)
names(s) <- c('intercept', 'slope')
plot(s)


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