## ---- echo=FALSE, include=FALSE-----------------------------------------------
library(knitr)
opts_chunk$set(
  fig.width = 8, fig.height = 5,  fig.cap='',  collapse = TRUE
)
opts_knit$set(
  progress = FALSE,
  global.par = TRUE
)

library(spdep)
library(rspatial)


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


## ---- loca2-------------------------------------------------------------------
library(spdep)
wr <- poly2nb(auck, row.names=auck$Id, queen=FALSE)
lstw <- nb2listw(wr, style='B')
Gi <- localG(auck$TB, lstw)
head(Gi)


## ---- loca3-------------------------------------------------------------------
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)


## ---- loca4-------------------------------------------------------------------
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)


## ---- loca6-------------------------------------------------------------------
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)


## ---- loca8-------------------------------------------------------------------
Ii <- localmoran(auck$TB, lstw)
Icuts <- cut(Ii, 5)
Icutsi <- as.integer(Icuts)
plot(auck, col=cols[Icutsi])
legend('bottomleft', levels(Icuts), fill=cols)


## ---- loca10------------------------------------------------------------------
cts <- sp_data('counties.rds')
p <- sp_data('precipitation.csv')
head(p)

plot(cts)
points(p[,c('LONG', 'LAT')], col='red', pch=20)


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


## ---- loca12------------------------------------------------------------------
m <- lm(pan ~ ALT, data=p)
m


## ---- loca13------------------------------------------------------------------
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)


## ---- loca14------------------------------------------------------------------
library( spgwr )
bw <- gwr.sel(pan ~ ALT, data=spt)
bw


## ---- loca16------------------------------------------------------------------
r <- raster(ctst, res=10000)
r <- rasterize(ctst, r)
newpts <- rasterToPoints(r)


## ---- loca17------------------------------------------------------------------
g <- gwr(pan ~ ALT, data=spt, bandwidth=bw, fit.points=newpts[, 1:2])
g


## ---- loca18, fig.width=9-----------------------------------------------------
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)