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


## ----warning=FALSE------------------------------------------------------------
library(terra)
library(rspat)
counties <- spat_data('counties')


## ----over1--------------------------------------------------------------------
yolo <- counties[counties$NAME == 'Yolo', ]
plot(counties, col='light gray', border='gray')
plot(yolo, add=TRUE, density=20, lwd=2, col='red')


## -----------------------------------------------------------------------------
rail <- spat_data("yolo-rail")
rail
class(rail)
city <- spat_data("city")

crs(yolo, TRUE)
crs(rail, TRUE)
crs(city, TRUE)


## -----------------------------------------------------------------------------
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)
yoloTA <- project(yolo, TA)
railTA <- project(rail, TA)
cityTA <- project(city, TA)


## -----------------------------------------------------------------------------
davis <- centroids(cityTA)
i <- which(relate(davis, countiesTA, "intersects"))
countiesTA$NAME[i]


## ----over5--------------------------------------------------------------------
i <- intersect(cityTA, countiesTA)
i$area <- expanse(i)
as.data.frame(i)

plot(cityTA, col='blue')
plot(yoloTA, add=TRUE, border='red', lwd=3)


## -----------------------------------------------------------------------------
davis_rail <- intersect(railTA, cityTA)


## ----over10-------------------------------------------------------------------
buf <- buffer(railTA, width=500)
buf <- aggregate(buf)
rail_buf <- intersect(buf, cityTA)

plot(cityTA, col='light gray')
plot(rail_buf, add=TRUE, col='light blue', border='light blue')
plot(railTA, add=TRUE, lty=2, lwd=6)
plot(cityTA, add=TRUE)
plot(davis_rail, add=TRUE, col='red', lwd=6)


## -----------------------------------------------------------------------------
round(100 * expanse(rail_buf) / expanse(cityTA))


## -----------------------------------------------------------------------------
parks <- spat_data("parks")
crs(parks, TRUE)
parksTA <- project(parks, TA)


## ----over12-------------------------------------------------------------------
plot(cityTA, col="light gray", border="light gray")
plot(railTA, add=TRUE, col="blue", lwd=4)
plot(parksTA, col="dark green", add=TRUE)

d <- distance(parksTA, railTA)
dmin <- apply(d, 1, min)
parksTA$railDist <- dmin

i <- which.max(dmin)
as.data.frame(parksTA)[i,]
plot(parksTA[i, ], add=TRUE, col="red", lwd=3, border="red")

j <- which.min(dmin)
as.data.frame(parksTA)[j,]
plot(parksTA[j, ], add=TRUE, col="red", lwd=3, border="orange")


## ----over14, fig.width=8, fig.height=4.5--------------------------------------
# use cityTA to set the geographic extent
r <- rast(cityTA)

# arbitrary resolution
dim(r) <- c(50, 100)

# rasterize the railroad lines
r <- rasterize(railTA, r)

# compute distance
d <- distance(r)
names(d) <- "dist"

# extract distance values for polygons
dp <- extract(d, parksTA, fun=mean, exact=TRUE)

dp <- data.frame(park=parksTA$PARK, distance=dp[,2])
dp <- dp[order(dp$distance), ]

plot(d)
plot(parksTA, add=TRUE)
plot(railTA, add=T, col="blue", lty=2)


## ----over16, fig.width=8------------------------------------------------------
centr <- centroids(parksTA)
v <- voronoi(centr)
plot(v)
points(centr, col="blue", pch=20)


## ----over18, fig.width=8------------------------------------------------------
vc <- intersect(v, cityTA)
plot(vc, border="red")
plot(parksTA, add=T, col="green")


## -----------------------------------------------------------------------------
alt <- spat_data("elevation")
alt


## ----over24-------------------------------------------------------------------
v <- extract(alt, yolo)
hist(unlist(v), main="Elevation in Yolo county")


## ----over26-------------------------------------------------------------------
# cut out a rectangle (extent) of Yolo
yalt <- crop(alt, yolo)
# "mask" out the values outside Yolo
ymask <- mask(yalt, yolo)

# summary of the raster cell values
summary(ymask)

plot(ymask)
plot(yolo, add=T)