11. Map overlay

11.1 Introduction

This document shows some example R code to do “overlays” and associated spatial data manipulation to accompany Chapter 11 in O’Sullivan and Unwin (2010). You have already seen many of this type of data manipulation in previsous labs. And we have done perhaps more advanced things using regression type models (including LDA and RandomForest). This lab is very much a review of what you have already seen: basic spatial data operations in R. Below are some of the key packages for spatial data analysis that we have been using.

sp - Defines classes (data structures) for points, lines, polygons, rasters, and their attributes, and related funcitons for e.g, plotting. The main classes are SpatialPoints, SpatialLines, SpatialPolygons, SpatialGrid (all also with ‘DataFrame’ appended if the geometries have attributes)

raster - Defines alternative classes for raster data (RasterLayer, RasterStack, RasterBrick) that can be used for very large data sets. The package also provides many functions to manipulate raster data. It also extends to sp and rgeos packages for manipulating vector type data

rgdal - Read or write spatial data files (raster uses it behind the scenes).

rgeos - Geometry manipulation for vector data (e.g. intersection of polygons) and related matters.

spatstat - The main package for point pattern analysis (here used for a density function).

Get the data

You can download the data for this tutorial here, or with the script below.

dir.create('data', showWarnings = FALSE)
files <- c('city.rds', 'counties.rds', 'parks.rds', 'yolo-rail.rds', 'elevation.tif')
for (filename in files) {
    localfile <- file.path('data', filename)
    if (!file.exists(localfile)) {
        download.file(paste0('https://biogeo.ucdavis.edu/data/rspatial/', filename), dest=localfile)
    }
}

11.2 Selection by attribute

By now, you are well aware that in R, polygons and their attributes can be represented by a ‘SpatialPolygonsDataFrame’ (a class defined in the sp package). Let’s read the shapefile of California counties.

library(raster)
counties <- readRDS('data/counties.rds')
## Warning in readRDS("data/counties.rds"): invalid or incomplete compressed
## data
## Error in readRDS("data/counties.rds"): error reading from connection

Selection by attribute of elements of a SpatialPolygonsDataFrame is similar to selecting rows from a data.frame. For example, to select Yolo county by its name:

yolo <- counties[counties$NAME == 'Yolo', ]
## Error in eval(expr, envir, enclos): object 'counties' not found
plot(counties, col='light gray', border='gray')
## Error in plot(counties, col = "light gray", border = "gray"): object 'counties' not found
plot(yolo, add=TRUE, density=20, lwd=2, col='red')
## Error in plot(yolo, add = TRUE, density = 20, lwd = 2, col = "red"): object 'yolo' not found

You can interactively select counties this way:

plot(counties)

s <- select(counties)

11.3 Intersection and buffer

I want to select the railroads in the city of Davis from the railroads in Yolo county. First read the data, and do an important sanity check: are the coordinate reference systems (“projections”) the same?

rail <- readRDS('data/yolo-rail.rds')
## Error in readRDS("data/yolo-rail.rds"): error reading from connection
rail
## Error in eval(expr, envir, enclos): object 'rail' not found
# removing attributes that I do not care about
rail <- geometry(rail)
## Error in geometry(rail): object 'rail' not found
class(rail)
## Error in eval(expr, envir, enclos): object 'rail' not found

city <- readRDS(file.path(datapath,'city.rds'))
## Error in file.path(datapath, "city.rds"): object 'datapath' not found
class(city)
## Error in eval(expr, envir, enclos): object 'city' not found
city <- geometry(city)
## Error in geometry(city): object 'city' not found
class(city)
## Error in eval(expr, envir, enclos): object 'city' not found

projection(yolo)
## Error in methods::extends(class(x), "BasicRaster"): object 'yolo' not found
projection(rail)
## Error in methods::extends(class(x), "BasicRaster"): object 'rail' not found
projection(city)
## Error in methods::extends(class(x), "BasicRaster"): object 'city' not found

Ay, we are dealing with two different coordinate reference systems (projections)! Let’s settle for yet another one: Teale Albers (this is really the “Albers Equal Area projection with parameters suitable for California”. This particular set of parameters was used by an California State organization called the Teale Data Center, hence the name.

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)
## Error in spTransform(counties, TA): object 'counties' not found
yoloTA <- spTransform(yolo, TA)
## Error in spTransform(yolo, TA): object 'yolo' not found
railTA <- spTransform(rail, TA)
## Error in spTransform(rail, TA): object 'rail' not found
cityTA <- spTransform(city, TA)
## Error in spTransform(city, TA): object 'city' not found

Another check, let’s see what county Davis is in, using two approaches. In the first one we get the centroid of Davis and do a point-in-polygon query.

dav <- coordinates(cityTA)
## Error in coordinates(cityTA): object 'cityTA' not found
davis <- SpatialPoints(dav, proj4string=TA)
## Error in coordinates(coords): object 'dav' not found
over(davis, countiesTA)
## Error in over(davis, countiesTA): object 'davis' not found

An alternative approach is to intersect the two polygon datasets.

i <- intersect(cityTA, countiesTA)
## Error in intersect(cityTA, countiesTA): object 'cityTA' not found

data.frame(i, area=area(i, byid=TRUE))
## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'area' for signature '"integer"'

plot(cityTA, col='blue')
## Error in plot(cityTA, col = "blue"): object 'cityTA' not found
plot(yoloTA, add=TRUE, border='red', lwd=3)
## Error in plot(yoloTA, add = TRUE, border = "red", lwd = 3): object 'yoloTA' not found

So we have a little sliver of Davis inside of Solano.

Everything looks OK. Now we can intersect rail and city, and make a buffer.

davis_rail <- intersect(railTA, cityTA)
## Error in intersect(railTA, cityTA): object 'railTA' not found

Compute a 500 meter buffer around railroad inside Davis:

buf <- buffer(railTA, width=500)
## Error in buffer(railTA, width = 500): object 'railTA' not found
rail_buf <- intersect(buf, cityTA)
## Error in intersect(buf, cityTA): object 'buf' not found

plot(cityTA, col='light gray')
## Error in plot(cityTA, col = "light gray"): object 'cityTA' not found
plot(rail_buf, add=TRUE, col='light blue', border='light blue')
## Error in plot(rail_buf, add = TRUE, col = "light blue", border = "light blue"): object 'rail_buf' not found
plot(railTA, add=TRUE, lty=2, lwd=6)
## Error in plot(railTA, add = TRUE, lty = 2, lwd = 6): object 'railTA' not found
plot(cityTA, add=TRUE)
## Error in plot(cityTA, add = TRUE): object 'cityTA' not found
plot(davis_rail, add=TRUE, col='red', lwd=6)
## Error in plot(davis_rail, add = TRUE, col = "red", lwd = 6): object 'davis_rail' not found
box()
## Error in box(): plot.new has not been called yet

What is the percentage of the area of the city of Davis that is within 500 m of a railroad?

round(100 * area(rail_buf) / area(cityTA))
## Error in area(rail_buf): object 'rail_buf' not found

11.4 Proximity

Which park in Davis is furthest, and which is closest to the railroad? First get the parks data.

parks <- readRDS('data/parks.rds')
## Warning in readRDS("data/parks.rds"): invalid or incomplete compressed data
## Error in readRDS("data/parks.rds"): error reading from connection
proj4string(parks)
## Error in proj4string(parks): object 'parks' not found
parksTA <- spTransform(parks, TA)
## Error in spTransform(parks, TA): object 'parks' not found

Now plot the parks that are the furthest and the nearest from a railroad.

plot(cityTA, col='light gray', border='light gray')
## Error in plot(cityTA, col = "light gray", border = "light gray"): object 'cityTA' not found
plot(railTA, add=T, col='blue', lwd=4)
## Error in plot(railTA, add = T, col = "blue", lwd = 4): object 'railTA' not found
plot(parksTA, col='dark green', add=TRUE)
## Error in plot(parksTA, col = "dark green", add = TRUE): object 'parksTA' not found

d <- gDistance(parksTA, railTA, byid=TRUE)
## Error in is.projected(spgeom1): object 'parksTA' not found
dmin <- apply(d, 2, min)
## Error in apply(d, 2, min): object 'd' not found
parksTA$railDist <- dmin
## Error in eval(expr, envir, enclos): object 'dmin' not found

i <- which.max(dmin)
## Error in which.max(dmin): object 'dmin' not found
data.frame(parksTA)[i,]
## Error in data.frame(parksTA): object 'parksTA' not found
plot(parksTA[i, ], add=TRUE, col='red', lwd=3, border='red')
## Error in plot(parksTA[i, ], add = TRUE, col = "red", lwd = 3, border = "red"): object 'parksTA' not found

j <- which.min(dmin)
## Error in which.min(dmin): object 'dmin' not found
data.frame(parksTA)[j,]
## Error in data.frame(parksTA): object 'parksTA' not found
plot(parksTA[j, ], add=TRUE, col='red', lwd=3, border='orange')
## Error in plot(parksTA[j, ], add = TRUE, col = "red", lwd = 3, border = "orange"): object 'parksTA' not found

Another way to approach this is to first create a raster with distance to the railroad values. Here we compute the average distance to any place inside the park, not to its border. You could also compute the distance to the centroid of a park.

library(raster)
# use cityTA to set the geogaphic extent
r <- raster(cityTA)
## Error in raster(cityTA): object 'cityTA' not found

# arbitrary resolution
dim(r) <- c(50, 100)
## Error in dim(r) <- c(50, 100): object 'r' not found

# rasterize the railroad lines
r <- rasterize(railTA, r, field=1)
## Error in rasterize(railTA, r, field = 1): object 'railTA' not found

# compute distance
d <- distance(r)
## Error in distance(r): object 'r' not found

# extract distance values for polygons
dp <- extract(d, parksTA, fun=mean, small=TRUE)
## Error in extract(d, parksTA, fun = mean, small = TRUE): object 'd' not found

dp <- data.frame(parksTA$PARK, dist=dp)
## Error in data.frame(parksTA$PARK, dist = dp): object 'parksTA' not found
dp <- dp[order(dp$dist), ]
## Error in eval(expr, envir, enclos): object 'dp' not found

plot(d)
## Error in plot(d): object 'd' not found
plot(parksTA, add=TRUE)
## Error in plot(parksTA, add = TRUE): object 'parksTA' not found
plot(railTA, add=T, col='blue', lty=2)
## Error in plot(railTA, add = T, col = "blue", lty = 2): object 'railTA' not found

Thiessen polygons

Here I compute Thiessen (or Voronoi) polygons for the Davis parks. Each polygon shows the area that is closest to (the centroid of) a particular park.

library(dismo)
centroids <- coordinates(parksTA)
## Error in coordinates(parksTA): object 'parksTA' not found
v <- voronoi(centroids)
## Loading required namespace: deldir
## Error in voronoi(centroids): object 'centroids' not found
plot(v)
## Error in plot(v): object 'v' not found
points(centroids, col='blue', pch=20)
## Error in points(centroids, col = "blue", pch = 20): object 'centroids' not found

To keep the polygons within Davis.

proj4string(v) <- TA
## Error in proj4string(v) <- TA: object 'v' not found
vc <- intersect(v, cityTA)
## Error in intersect(v, cityTA): object 'v' not found
plot(vc, border='red')
## Error in plot(vc, border = "red"): object 'vc' not found
plot(parksTA, add=T, col='green')
## Error in plot(parksTA, add = T, col = "green"): object 'parksTA' not found

11.5 Fields

Raster data

raster data can be read with readGDAL to get a SpatialGridDataFrame. But I prefer to use the raster package to create Raster* objects. Raster* stands for RasterLayer, RasterStack, or RasterBrick. See the vignette for the raster package for more details: http://cran.r-project.org/web/packages/raster/vignettes/Raster.pdf

library(raster)
# create a RasterLayer object from a file
alt <- raster("data/elevation.tif")
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file.
alt
## Error in eval(expr, envir, enclos): object 'alt' not found

RasterLayer ‘alt’ is in lon/lat, and so is ‘yolo’. However, ‘yolo’ has a different datum (NAD83) than ‘alt’ (WGS84). While there is no real difference between these, this will lead to errors, so we first transform yolo. It is generally better to match vector data to raster data than vice versa.

yolo <- spTransform(yolo, crs(alt))
## Error in spTransform(yolo, crs(alt)): object 'yolo' not found
plot(alt)
## Error in plot(alt): object 'alt' not found
plot(yolo, add=TRUE)
## Error in plot(yolo, add = TRUE): object 'yolo' not found

Shaded relief is always nice to look at.

slope <- terrain(alt, opt='slope')
## Error in nlayers(x): object 'alt' not found
aspect <- terrain(alt, opt='aspect')
## Error in nlayers(x): object 'alt' not found
hill <- hillShade(slope, aspect, 40, 270)
## Error in compareRaster(slope, aspect): object 'slope' not found
plot(hill, col=grey(0:100/100), legend=FALSE, main='Elevation')
## Error in plot(hill, col = grey(0:100/100), legend = FALSE, main = "Elevation"): object 'hill' not found
plot(alt, col=rainbow(25, alpha=0.35), add=TRUE)
## Error in plot(alt, col = rainbow(25, alpha = 0.35), add = TRUE): object 'alt' not found

You can also try the plot3D function in the rasterVis package.

Query

Now extract elevation data for Yolo county.

v <- extract(alt, yolo)
## Error in extract(alt, yolo): object 'alt' not found
hist(unlist(v), main='Elevation in Yolo county')
## Error in unlist(v): object 'v' not found

Another approach:

# cut out a rectangle (extent) of Yolo
yalt <- crop(alt, yolo)
## Error in crop(alt, yolo): object 'alt' not found
# 'mask' out the values outside Yolo
ymask <- mask(yalt, yolo)
## Error in mask(yalt, yolo): object 'yalt' not found

# summary of the raster cell values
summary(ymask)
## Error in summary(ymask): object 'ymask' not found

plot(ymask)
## Error in plot(ymask): object 'ymask' not found
plot(yolo, add=T)
## Error in plot(yolo, add = T): object 'yolo' not found

You can also get values (query) by clicking on the map (use click(alt))

11.6 Exercise

We want to travel by train in Yolo county and we want to get as close as possible to a hilly area; whenever we get there we’ll jump from the train. It turns out that the railroad tracks are not all connected, we will ignore that inconvenience.

Define “hilly” as some vague notion of a combination of high elevation and slope (in degrees). Use some of the functions you have seen above, as well as function ‘distance’ to create a plot of distance from the railroad against hillyness. Make a map to illustrate the result, showing where you get off the train, where you go to, and what the elevation and slope profile would be if you follow the shortest (as the crow flies) path.

Bonus: use package gdistance to find the least-cost path between these two points (assign a cost to slope, perhaps using Tobler’s hiking function).