## ----getData, echo=TRUE-------------------------------------------------------
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')


## -----------------------------------------------------------------------------
library(rspatial)
library(raster)
counties <- sp_data('counties.rds')


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


## -----------------------------------------------------------------------------
rail <- sp_data('yolo-rail.rds')
rail
# removing attributes that I do not care about
rail <- geometry(rail)
class(rail)

city <- sp_data('city.rds')
class(city)
city <- geometry(city)
class(city)

projection(yolo)
projection(rail)
projection(city)


## -----------------------------------------------------------------------------
library(rgdal)
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 +ellps=GRS80")
countiesTA <- spTransform(counties, TA)
yoloTA <- spTransform(yolo, TA)
railTA <- spTransform(rail, TA)
cityTA <- spTransform(city, TA)


## -----------------------------------------------------------------------------
dav <- coordinates(cityTA)
davis <- SpatialPoints(dav, proj4string=TA)
over(davis, countiesTA)


## ---- over5-------------------------------------------------------------------
i <- intersect(cityTA, countiesTA)

data.frame(i, area=area(i, byid=TRUE))

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


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


## ---- over10------------------------------------------------------------------
buf <- buffer(railTA, width=500)
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)
box()


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


## -----------------------------------------------------------------------------
parks <- sp_data('parks.rds')
proj4string(parks)
parksTA <- spTransform(parks, TA)


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

library(rgeos)
d <- gDistance(parksTA, railTA, byid=TRUE)
dmin <- apply(d, 2, min)
parksTA$railDist <- dmin

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

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


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

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

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

# compute distance
d <- distance(r)

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

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

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


## ---- over16, fig.width=8-----------------------------------------------------
library(dismo)
centroids <- coordinates(parksTA)
v <- voronoi(centroids)
plot(v)
points(centroids, col='blue', pch=20)


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


## -----------------------------------------------------------------------------
library(raster)
# get a RasterLayer with elevation data
alt <- sp_data("elevation")
alt


## ---- over20, fig.width=7-----------------------------------------------------
yolo <- spTransform(yolo, crs(alt))
plot(alt)
plot(yolo, add=TRUE)


## ---- over22, fig.width=7-----------------------------------------------------
slope <- terrain(alt, opt='slope')
aspect <- terrain(alt, opt='aspect')
hill <- hillShade(slope, aspect, 40, 270)
plot(hill, col=grey(0:100/100), legend=FALSE, main='Elevation')
plot(alt, col=rainbow(25, alpha=0.35), add=TRUE)


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