This and other labs in this course are based on those developed by C.J. Gabbe and Iris Stewart Frey.
Create a Module 6 Lab folder
Load required packages. You may need to install any packages that you
have not used so far before you load them, such as terra,
sp, tigris, and/or FedData.
#install.packages("FedData")
#install.packages("terra")library(sf)
library(tidyverse)
library(tmap)
library(lwgeom)
library(terra)
library(sp)
library(tigris)
library(FedData)Set working directory:
setwd("/Users/istewartfrey/Google Drive/Teaching/Module6_Lab")Read the statewide Mobile Home Parks data that we worked with in the
Module 5 lab, st_transform the data to California State
Plane Zone 3, and filter so only Tuolumne County MHPs are
included. (You should already have the code for this.)
We’re going to download elevation data for Tuolumne County, CA using
the extremely useful FedData package. We first need to give
the FedData package a boundary (in this case Tuolumne
County).
We can first use the tigris package (which you’ll learn
about in the next DataCamp homework) to download the Tuolumne County
boundary and the Sonora City boundary from the Census Bureau,
CA_counties <- counties(state = "CA") # Download county boundaries for CA
Tuolumne <- filter(CA_counties, NAME == "Tuolumne")
Tuolumne <-st_transform(Tuolumne, crs = 2227)# We will also extract the boundary for the city of Sonora, within Tuolumne County. We will use this boundary to visualize where MHPs are located with respect to city boundaries (and city services).
CA_places <- places(state = "CA") # Download city boundaries for CA
Sonora <- filter(CA_places, NAME == "Sonora")
Sonora <-st_transform(Sonora, crs = 2227)
#Sonora.sp <- as_Spatial(Sonora) # This converts to an sp object
#Sonora.sp <- sp::spTransform(Sonora.sp, CRS("+init=epsg:2227")) # project to California State Plane Zone 3Then, use the FedData package to download the elevation
data.
ned_tuo <- get_ned(template = Tuolumne, label = "Tuolumne", res = "1")Check the raster’s CRS, extent, and resolution. You can do this by
typing ned_tuo into the console.
When we’re working with vector or raster data, we need to make sure
all datasets are in the same CRS (coordinate reference system). Try
using project to project your data to California State
Plane Zone 3. Type ?project if you need help with the
syntax. The only arguments that you need are the raster istelf and then
in quotes you can specify the string for the target projection.
ned_tuo <- project(ned_tuo, "+init=epsg:2227")Tip: sometimes if the raster file is large and/or you face another issue with projecting, you might consider projecting the vector (sf) datasets into the same coordinate system of the raster instead.
There are many different ways to start exploring the values in your raster dataset.
First try typing the name of the data object into the console. Then
try using hist and plot to get a better
understanding of your data. (What units of measure do you think these
values are in?)
plot(ned_tuo)Masking essentially means clipping our raster. Let’s use our Tuolumne
boundary (that you downloaded and converted to sp above) to
mask the raster. After you’ve applied the mask, create a
plot to verify that it worked. You should see the raster
now in the shape of Tuolumne County. If masking doesn’t work, the most
common explanation is that there is a mismatch in CRS
between the datasets.
ned_tuo_mask <- mask(ned_tuo, mask = Tuolumne) # mask
plot(ned_tuo_mask) # make sure it worked
writeRaster(ned_tuo_mask, "NED_Tuolumne.tif", overwrite = TRUE) # write the raster to file - we will use it for plottingYou’re going to use extract from the terra
package to connect each MHP with its underlying elevation value. First
we need to read in the MHPS that we have used in previous labs and
select those for Toulumne County
MHPs <-st_read("MHPs_HCD_022818_CLEAN.shp")
MHPs_TUO <-subset(MHPs ,Conty_Nm == "TUOLUMNE")Now we can use extract
MHPs_TUO.elev <- MHPs_TUO %>% mutate(Elev_m = raster::extract(ned_tuo_mask, MHPs_TUO, fun = mean, na.rm = TRUE))
st_write(MHPs_TUO.elev, dsn = "MHPs_Elev_Tuolumne.shp", driver = "ESRI Shapefile", append=FALSE)There are a few things going on in the code above, so try explaining what each step does in plain language.
Create a map to visually assess the data. Your map should have a
basemap, the elevation data for Tuolumne County, the Mobile Home parks
classified by lot number, and the outline of the city of Sonora. You can
use the tmap package to map our raster. Try adding
tm_raster to the syntax we’ve developed previously. Here’s
a piece of code you can add to your interactive mapping syntax:
tmap_mode("view")
tm_basemap(leaflet::providers$Esri.WorldImagery, group = "Esri.WorldImagery") +
tm_shape(ned_tuo_mask) +
tm_raster(col = "USGS_1_n38w120", alpha = 0.8, palette = "Reds", title = "Elevation (m)") +
tm_tiles(leaflet::providers$CartoDB.PositronOnlyLabels, group = "CartoDB labels") +
tm_shape(Sonora)+
tm_borders() +
tm_shape(MHPs_TUO.elev) +
tm_symbols(col = "MH_Lots", size = 0.1, n = 4, style = "quantile")And here is a piece that adds the polygon outline only:
tmap_mode("view")
tm_basemap(leaflet::providers$Esri.WorldImagery, group = "Esri.WorldImagery") +
tm_shape(ned_tuo_mask) +
tm_raster(col = "USGS_1_n38w120", alpha = 0.8, palette = "Reds", title = "Elevation (m)") +
tm_tiles(leaflet::providers$CartoDB.PositronOnlyLabels, group = "CartoDB labels") +
tm_shape(Sonora)+
tm_borders() +
tm_shape(MHPs_TUO.elev) +
tm_symbols(col = "MH_Lots", size = 0.1, n = 4, style = "quantile")+
tm_shape(Sonora)+
tm_borders()