This and other labs in this course are based on those developed by C.J. Gabbe and Iris Stewart Frey
Create a ‘Module4_Lab’ folder in the drive of your choice (preferably somewhere that’s backed up into the cloud)
Create a new RScript by going to File > New File > RScript. Save the RScript to your ‘Module4_Lab’ folder.
(As a reminder it’s good coding form to leave the install.packages code out of your RScript) We’re going to be using some new packages today. Google them to see what each one does:
install.packages("tmap")
install.packages("tmaptools")
install.packages("shinyjs")
install.packages("leaflet")
install.packages("lwgeom")Load the packages we’ll need (again, you have to load packages this every time R restarts)
library(tidyverse)
library(sf)
library(tmap)
library(tmaptools)
library(shinyjs)
library(leaflet)
library(lwgeom)Your working directory will be the Module 4 lab folder, similar to:
setwd("~/Google Drive/My Drive/ENVS 117/04_Labs/Module4_Lab")or you could set it up from your desktopsetwd(“~/Desktop/ENVS117/lab4”)
Again, you can always check your working directory by typing
getwd() into the console.
We’re going to be exploring the state of California’s Fire Hazard Severity Zone Maps. Scroll down to see where these maps are given by county. More fire-risk related data can be found here.
Please scroll down on this page to sonoma county [https://osfm.fire.ca.gov/what-we-do/community-wildfire-preparedness-and-mitigation/fire-hazard-severity-zones/fire-hazard-severity-zones-maps] and download the shapefile for the State Responsibility Area for Sonoma County. Unzip and move it into your Module 4 lab folder. Read the shapefile into R:
SRA_Sonoma <- st_read("fhszs06_3_49.shp")summary(SRA_Sonoma)What’s the name of the most relevant variable if we’re looking to understand wildfire risk? What are the three classes (categories) of this variable? What type of variable is it?
Let’s make a very basic plot of this variable:
plot(SRA_Sonoma["HAZ_CLASS"])Just as in with ArcGIS, we must be mindful of the geographic and projected coordinate systems for our data. In “R lingo” this is called the coordinate reference system (CRS). First, let’s check the CRS for our wildfire hazard data, you can do so by running the file in R chunk:
SRA_SonomaThe EPSG and proj4string codes give us an indication of the CRS of our data. From the proj4string value, we know that our data include some projection info (which is good!) but we want our data to contain an EPSG code too (another posotive of this is proj4strings are being slowly decommisoned in favor of strictly EPSG. Remember in ENVS 116 when we projected data into different projected coordinate systems? (For example, California State Plane Zone 3 for the Bay Area.) Each one of those named coordinate systems has a corresponding EPSG number and proj4string code. Just like in lab 3, you can search online for the right codes. What is EPSG:3310? You can also read more about proj4 here.
We’re going to project our data to EPSG:2227 (what is this CRS?):
SRA_Sonoma <- st_transform(SRA_Sonoma, crs = 2227) # project to California State Plane Zone 3Try explaining in simple English what the above code did. If you want
to learn more, try typing ?st_transform into the
console.
You can check the CRS again to make sure the projection worked:
SRA_SonomaWe’re going to use the tmap package to make some maps.
First, we’re going to make a simple QTM:. We may have some topological
errors to fix this, try this from the lwgeom package and then re-run the
plot: SRA_Sonoma <- st_make_valid(SRA_Sonoma) If you get an error
message about invalid polygons try the code above again
tmap_mode("view")
SRA_Sonoma <- st_make_valid(SRA_Sonoma)
quickmap1 <- qtm(SRA_Sonoma, fill = "HAZ_CLASS", check.and.fix = TRUE) # This creates the map object
quickmap1 # This will display the map in the plots section of RStudioWe’re going to use the tmap package to create some nicer
maps with lots of different colors. There’s a good primer to
tmap here.
Try this code to see some of your color palette options:
tmaptools::palette_explorer() # Use this to find your color palette optionsUsing tmap is a little like making a sandwich, where we
add different components and change the ingredients through the code.
We’ll create a first version of our map by taking our shapefile and then
adding polygons that display the ‘HAZ_CLASS’ variable with the red color
palette.
tmap_mode("view")
hazmap_v1 <- tm_shape(SRA_Sonoma) + tm_fill(col = "HAZ_CLASS", palette = "Reds")
hazmap_v1So what do you think? Looks good, but there’s something wrong - what is it?
The levels aren’t in the right order! Let’s reorder them so they go Moderate, High, Very High
SRA_Sonoma$HAZ_CLASS <- ordered(SRA_Sonoma$HAZ_CLASS, levels = c("Moderate", "High", "Very High"))
summary(SRA_Sonoma$HAZ_CLASS) # This will let you check the levels to make sure it workedAgain, try describing in plain English what you just did. Make sure to mention what ‘SRA’, ‘HAZ_CLASS’, ‘ordered’, and ‘levels’ mean.
There are a lot of basemaps you can choose from; first, let’s figure out what ones are available (from the Leaflet package):
providersIf you want to check them out online, there’s a website you can use to preview the options.
Let’s try
#tmap_mode("view")
#hazmap_v2 <- tm_shape(SRA_Sonoma) + tm_fill(col = "HAZ_CLASS", palette = "Reds") + tm_basemap(server = "CartoDB.Positron")
#hazmap_v2What does tmap_mode do? What is the difference between
“view” and “plot” options? Try using a different basemap too.
TIP: If you ever get an error that a geometry is invalid, you can
test it using st_is_valid and fix it using
lwgeom::st_make_valid
#10 Adding map elements Remember in ENVS 116 we had the acronym TALNS. To add map elements we have to switch to plot which yields a non interactive map.
Lets start easy and add a scale bar (keep in mind if you are copying and pasting these codes you need to remove the hashtags (#))
#tmap_mode("plot")
#tm_shape(SRA_Sonoma) + tm_fill(col = "HAZ_CLASS", palette = "Reds")+
# tm_scale_bar(breaks = c(0, 5, 10), text.size = 2, position = c("right", "bottom"))Now lets try and add a title (keep in mind if you are copying and pasting these codes you need to remove the hashtags (#))
#tmap_mode("plot")
#tm_shape(SRA_Sonoma) + tm_fill(col = "HAZ_CLASS", palette = "Reds")+
# tm_layout(main.title = "Fire Severity in Sonoma County",
# main.title.size = 0.95, main.title.position="center",
# legend.outside = TRUE) +
# tm_scale_bar(breaks = c(0, 5, 10), text.size = 2, position = c("right", "bottom"))Now lets try to add a north arrow and increase the title size(keep in mind if you are copying and pasting these codes you need to remove the hashtags (#))
#tmap_mode("plot")
#tm_shape(SRA_Sonoma) + tm_fill(col = "HAZ_CLASS", palette = "Reds")+
# tm_layout(main.title = "Fire Severity in Sonoma County",
# main.title.size = 1.95, main.title.position="center",
# legend.outside = TRUE) +
#tm_scale_bar(breaks = c(0, 5, 10), text.size = 2, position = c("right", "bottom"))+
#tm_compass(type = "4star", position = c("left", "bottom")) We will learn to add map more elements later on this quarter while making an interactive map but for now lets subset some of our data
Sometimes we want to create a subset of our dataset by filtering out the data that we’re interested in. Let’s say we want to create a shapefile that only include polygons with very high wildfire risk:
SRA_veryhigh <- filter(SRA_Sonoma, HAZ_CLASS == "Very High")
summary(SRA_veryhigh) # Quick check to make sure that it workedMap an interactive version on of our subset. (keep in mind if you are copying and pasting these codes you need to remove the hashtags (#))
#tmap_mode("view")
#hazmap_v3 <- tm_shape(SRA_veryhigh) + tm_fill(col = "HAZ_CLASS", palette = "Reds") + tm_basemap(server = "CartoDB.Positron")
#hazmap_v3It’s always a good idea to clean things up and save at the end of every session. For example, you might want to create some categories like these:
# Set up
# Data import
# Data management
# MappingThen, you could move the code where we reordered the HAZ_CLASS variable to the #Data import section rather than having it in the middle of your mapping.
```