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


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

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

library(rspatial)
counties <- sp_data('counties')
p <- sp_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 <- 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(counties, 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)


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


## -----------------------------------------------------------------------------
library(sp)
coordinates(houses) <- ~longitude+latitude


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


## -----------------------------------------------------------------------------
library(raster)
crs(houses) <- crs(counties)


## -----------------------------------------------------------------------------
cnty <- over(houses, counties)
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, by='NAME')
dim(dcounties)


## ---- gwr5--------------------------------------------------------------------
cnres <- merge(dcounties, resdf, by='NAME')
spplot(cnres, 'income')


## ---- gwr6--------------------------------------------------------------------
# a copy of the data
cnres2 <- cnres

# scale all variables, except the first one (county name)
# assigning values to a "@data" slot is risky, but (I think) OK here
cnres2@data = data.frame(scale(data.frame(cnres)[, -1]))
spplot(cnres2)


## ---- gwr10-------------------------------------------------------------------
library(spdep)
nb <- poly2nb(cnres)
plot(cnres)
plot(nb, coordinates(cnres), add=T, col='red')

lw <- nb2listw(nb)
moran.test(cnres$income, lw)
moran.test(cnres$roomhead, lw, na.action=na.omit)


## -----------------------------------------------------------------------------
TA <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000
              +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
countiesTA <- spTransform(counties, TA)


## -----------------------------------------------------------------------------
library(raster)
r <- raster(countiesTA)
res(r) <- 50000


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


## -----------------------------------------------------------------------------
housesTA <- spTransform(houses, TA)
crds <- coordinates(housesTA)


## -----------------------------------------------------------------------------
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)
Moran(rinc)


## -----------------------------------------------------------------------------
r <- raster(countiesTA)


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


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


## -----------------------------------------------------------------------------
fitpoints <- rasterToPoints(ca)


## -----------------------------------------------------------------------------
fitpoints <- fitpoints[,-3]