## ----loca1, message=FALSE-----------------------------------------------------
if (!require("rspat")) remotes::install_github("rspatial/rspat")
library(rspat)
auck <- spat_data("auctb")


## ----locasp-------------------------------------------------------------------
w <- relate(auck, relation="rook")


## ----loca2--------------------------------------------------------------------
Gi <- autocor(auck$TB, w, "Gi")
head(Gi)


## ----loca3, fig.height=4------------------------------------------------------
grays <- rev(gray(seq(0,1,.2)))
auck$Gi <- Gi
plot(auck, "Gi", col=grays)


## ----loca4, fig.height=4------------------------------------------------------
#include "self"
diag(w) <- TRUE
auck$Gistar <- autocor(auck$TB, w, "Gi*")
plot(auck, "Gistar", main="Gi*", col=grays)


## ----loca6, fig.height=4------------------------------------------------------
auck$loc_mean <- apply(w, 1, function(i) mean(auck$TB[i]))
plot(auck, "loc_mean", main="Local mean", col=grays)


## ----loca8, message=FALSE, fig.height=4---------------------------------------
library(spdep)
sfauck <- sf::st_as_sf(auck)
wr <- poly2nb(sfauck, row.names=sfauck$Id, queen=FALSE)
lstw <- nb2listw(wr, style="B")
auck$Ii <- localmoran(auck$TB, lstw)
plot(auck, "Ii", main="Local Moran", col=grays)


## -----------------------------------------------------------------------------
auck$TBdens <- 10000 * auck$TB / expanse(auck)


## ----loca10-------------------------------------------------------------------
cts <- spat_data("counties")
p <- spat_data("precipitation")
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 <- "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=WGS84"
sp <- vect(p, c("LONG", "LAT"), crs="+proj=longlat +datum=NAD83")
sp <- terra::project(sp, alb)
spsf <- sf::st_as_sf(sp)
vctst <- terra::project(cts, alb)
ctst <- sf::st_as_sf(vctst)


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


## ----loca16-------------------------------------------------------------------
r <- rast(vctst, res=10000)
r <- rasterize(vect(ctst), r)
newpts <- geom(as.points(r))[, c("x", "y")]


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


## ----loca18, fig.width=9------------------------------------------------------
slope <- intercept <- r
slope[!is.na(slope)] <- g$SDF$ALT 
intercept[!is.na(intercept)] <- g$SDF$"(Intercept)"
s <- c(intercept, slope)
names(s) <- c("intercept", "slope")
plot(s)