## ----getDataLocal-------------------------------------------------------------

if (!require("rspat")) remotes::install_github('rspatial/rspat')

library(rspat)
counties <- spat_data("counties")
p <- spat_data("precipitation")
head(p)

plot(counties)
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 +units=m"
sp <- vect(p, c("LONG", "LAT"), crs="+proj=longlat +datum=WGS84")
spt <- project(sp, alb)
ctst <- project(counties, alb)


## ----loca14-------------------------------------------------------------------
library( spgwr )
bw <- gwr.sel(pan ~ ALT, data=as.data.frame(spt), coords=geom(spt)[,c("x", "y")])
bw


## ----loca16-------------------------------------------------------------------
r <- rast(ctst, res=10000)
r <- rasterize(ctst, r)
newpts <- as.points(r)


## ----loca17-------------------------------------------------------------------
g <- gwr(pan ~ ALT, data=as.data.frame(spt), coords=geom(spt)[,c("x", "y")], bandwidth=bw, fit.points=geom(newpts)[,c("x", "y")])
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)


## -----------------------------------------------------------------------------
houses <- spat_data("houses1990.csv")
dim(houses)
head(houses)


## -----------------------------------------------------------------------------
hvect <- vect(houses, c("longitude", "latitude"))


## ----gwr1---------------------------------------------------------------------
plot(hvect, cex=0.5, pch=1, axes=TRUE)


## -----------------------------------------------------------------------------
crs(hvect) <- crs(counties)


## -----------------------------------------------------------------------------
cnty <- extract(counties, hvect)
head(cnty)


## -----------------------------------------------------------------------------
hd <- cbind(data.frame(houses), cnty)


## -----------------------------------------------------------------------------
totpop <- tapply(hd$population, hd$NAME, sum)
totpop


## -----------------------------------------------------------------------------
# total income
hd$suminc <- hd$income * hd$households
# now use aggregate (similar to tapply)
csum <- aggregate(hd[, c('suminc', 'households')], list(hd$NAME), sum)
# divide total income by number of housefholds
csum$income <- 10000 * csum$suminc / csum$households
# sort
csum <- csum[order(csum$income), ]
head(csum)
tail(csum)


## -----------------------------------------------------------------------------
hd$roomhead <- hd$rooms / hd$population
hd$bedroomhead <- hd$bedrooms / hd$population
hd$hhsize <- hd$population / hd$households


## -----------------------------------------------------------------------------
# OLS
m <- glm( houseValue ~ income + houseAge + roomhead + bedroomhead + population, data=hd)
summary(m)
coefficients(m)


## -----------------------------------------------------------------------------
 hd2 <- hd[!is.na(hd$NAME), ]


## -----------------------------------------------------------------------------
regfun <- function(x)  {
  dat <- hd2[hd2$NAME == x, ]
  m <- glm(houseValue~income+houseAge+roomhead+bedroomhead+population, data=dat)
  coefficients(m)
}


## -----------------------------------------------------------------------------
countynames <- unique(hd2$NAME)
res <- sapply(countynames, regfun)


## ----gwr3, fig.height=10------------------------------------------------------
dotchart(sort(res["income", ]), cex=0.65)


## -----------------------------------------------------------------------------
resdf <- data.frame(NAME=colnames(res), t(res))
head(resdf)


## -----------------------------------------------------------------------------
dim(counties)
dcounties <- aggregate(counties[, "NAME"], "NAME")
dim(dcounties)


## ----gwr5---------------------------------------------------------------------
cnres <- merge(dcounties, resdf, by="NAME")
plot(cnres, "income")


## ----gwr6, fig.width=10-------------------------------------------------------
# a copy of the data
cnres2 <- cnres

# scale all variables, except the first one (county name)
values(cnres2) <- as.data.frame(scale(as.data.frame(cnres)[,-1]))
plot(cnres2, names(cnres2)[1:6], plg=list(x="topright"), mar=c(1,1,1,1))


## ----gwr10--------------------------------------------------------------------
lw <- adjacent(cnres2, pairs=FALSE)
autocor(cnres$income, lw)
autocor(cnres$houseAge, lw)


## -----------------------------------------------------------------------------
TA <- "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=WGS84 +units=m"
countiesTA <- project(counties, TA)


## -----------------------------------------------------------------------------
r <- rast(countiesTA)
res(r) <- 50000


## -----------------------------------------------------------------------------
xy <- xyFromCell(r, 1:ncell(r))


## -----------------------------------------------------------------------------
housesTA <- project(hvect, TA)
crds <- geom(housesTA)[, c("x", "y")]


## -----------------------------------------------------------------------------
regfun2 <- function(d)  {
 m <- glm(houseValue~income+houseAge+roomhead+bedroomhead+population, data=d)
 coefficients(m)
}


## -----------------------------------------------------------------------------
res <- list()
for (i in 1:nrow(xy)) {
	d <- sqrt((xy[i,1]-crds[,1])^2 + (xy[i,2]-crds[,2])^2)
	j <- which(d < 50000)
	if (length(j) > 49) {
		d <- hd[j,]
		res[[i]] <- regfun2(d)
	} else {
		res[[i]] <- NA
	}
}


## -----------------------------------------------------------------------------
inc <- sapply(res, function(x) x['income'])


## ----gwr20--------------------------------------------------------------------
rinc <- setValues(r, inc)
plot(rinc)
plot(countiesTA, add=T)
autocor(rinc)


## -----------------------------------------------------------------------------
r <- rast(countiesTA)


## -----------------------------------------------------------------------------
res(r) <- 25000


## -----------------------------------------------------------------------------
ca <- rasterize(countiesTA, r)


## -----------------------------------------------------------------------------
fitpoints <- crds(ca)