Load required packages:

#install.packages("rgdal")
#install.packages("maps")
#install.packages("maptools")
#install.packages("sp")
#install.packages("RColorBrewer")
library(rgdal)
library(maptools)
library(sp)
library(tidyverse)
library(maps)
library(RColorBrewer)
library(cowplot)

Exercise 1: Mapping in base R

This map of the U.S. was found here.

# First, we need to read in a shapefile using the readOGR command
# dsn = data source name; folder/directory your shapefile is in
# layer is shapefile name without .shp
usa <- readOGR(dsn = "shapefiles", layer = "cb_2015_us_state_500k")
## OGR data source with driver: ESRI Shapefile 
## Source: "shapefiles", layer: "cb_2015_us_state_500k"
## with 49 features
## It has 9 fields
plot(usa)

Exercise 2: Mapping in ggplot2

# Code from presentation to build simple world map in ggplot:
ggworld <- map_data("world") # built-in maps within ggplot; can specify anything within the quotation marks, e.g. "usa", "italy"

# Take a look at what the map looks like using ggplot (notice that we retain general plotting formatting, regardless of the fact that we're plotting a map):
ggplot(ggworld, aes(x=long, y=lat, group=group)) +
  geom_polygon()

# We can add on theme_nothing to get rid of axes (since they are somewhat extraneous in maps)
ggplot(ggworld, aes(x=long, y=lat, group=group)) +
  geom_polygon() +
  theme_nothing()

Exercise 3: Federally protected lands

Let’s begin by constructing a map of the U.S. where we’ll be able to overlay information about the federally protected lands.

# Using the template above, build a map with state lines ("state") using the map_data function
# YOUR CODE HERE
ggstate <- map_data("state")
ggplot(ggstate, aes(x=long, y=lat, group=group)) + 
  geom_polygon() +
  theme_nothing()

We now have our background state map ready to go. We eventually want to colorize each state by how many visitors it brings in for the national parks, and also by how large an area of federally protected land it has. First, we need to import and mold the data to map it to the U.S. map.

# Read in the data and take a look at it:
NPSdata <- read.table("data/NPS_data.txt", header=T)
head(NPSdata)
##   visitation_2015                  name  region   area_km2 category
## 1           19430       Russell_Cave_NM alabama     1.2563   NatMon
## 2           17818  Lake_Clark_NP_&_PRES  alaska 10601.6006  NatPark
## 3          551353 Glacier_Bay_NP_&_PRES  alaska 13044.5699  NatPark
## 4           37818      Katmai_NP_&_PRES  alaska 14869.6398  NatPark
## 5          560757      Denali_NP_&_PRES  alaska 19185.7868  NatPark
## 6               0   Cape_Krusenstern_NM  alaska  2626.9138   NatMon
# We need a single value for each state for both visitation and area:
NPSdata %>% 
  group_by(region) %>% 
  summarise(visitation = sum(visitation_2015), area = sum(area_km2)) -> NPSsumdata

# Now, we need to append the data to our map using a column shared between our map (ggstate) and our data (NPSsumdata):
allNPS <- left_join(ggstate, NPSsumdata, by="region")
head(allNPS)
##        long      lat group order  region subregion visitation   area
## 1 -87.46201 30.38968     1     1 alabama      <NA>      19430 1.2563
## 2 -87.48493 30.37249     1     2 alabama      <NA>      19430 1.2563
## 3 -87.52503 30.37249     1     3 alabama      <NA>      19430 1.2563
## 4 -87.53076 30.33239     1     4 alabama      <NA>      19430 1.2563
## 5 -87.57087 30.32665     1     5 alabama      <NA>      19430 1.2563
## 6 -87.58806 30.32665     1     6 alabama      <NA>      19430 1.2563

Now it’s time to build the map!

# Colorized states based on size of federal lands:
ggplot(allNPS, aes(x=long, y=lat, group=group, fill=area)) + 
  geom_polygon()

# Let's add custom coloration using RColorBrewer and add black outlines around states:
ggplot(allNPS, aes(x=long, y=lat, group=group, fill=area)) + 
  geom_polygon(color="black") +
  scale_fill_gradient2(low="dodgerblue1", high="tomato2", midpoint = 12500)

# We can also remove the axes but retain the legend:
ggplot(allNPS, aes(x=long, y=lat, group=group, fill=area)) + 
  geom_polygon(color="black") +
  scale_fill_gradient2(low="dodgerblue1", high="tomato2", midpoint = 12500) +
  theme(axis.title=element_blank(),
        axis.text=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank())

# Colorize states based on visitation to that state's parks:
# YOUR CODE HERE
ggplot(allNPS, aes(x=long, y=lat, group=group, fill=visitation)) + 
  geom_polygon(color="black") +
  scale_fill_gradient2(low="dodgerblue1", high="tomato2", midpoint = 6000000)

Exercise 4: Using Instagram in Yosemite

We can query Instagram based on a hashtag (#yosemite) this year and determine where photos were taken within a given vicinity (Yosemite National Park). I have obtained these data from John Farrell (see this post for his code and additional maps).

For this, we’ll need three things: 1. Shapefile of Yosemite NP (obtained here) 1. Point data (coordinates) of where photos were taken 1. Roads of U.S. to overlay and compare

Let’s prepare shapefiles for entry into ggplot using the fortify() command.

# Read in the shapefile using readOGR, as before:
yosemite <- readOGR(dsn = "shapefiles", layer = "nps_boundary")
## OGR data source with driver: ESRI Shapefile 
## Source: "shapefiles", layer: "nps_boundary"
## with 1 features
## It has 11 fields
# Get shapefile into format readable by ggplot:
yosemite <- fortify(yosemite, region="REGION")

# Plot the park as we did before with the U.S.; make outline black and center whatever color you want:
# YOUR CODE HERE
ggplot(yosemite, aes(x=long, y=lat, group=group)) + 
         geom_polygon(fill="white", color="black")

# Import our Instagram data as coordinates:
coordinates <- read_csv("data/coordinates.csv")

# Add Instagram point data onto map of Yosemite (notice that we need to remove grouping):
ggplot(yosemite, aes(x=long, y=lat)) + 
  geom_polygon(fill="white", color="black") +
  geom_point(data=coordinates, aes(x=long, y=lat)) +
  coord_fixed(ratio=1.3)

# Import road data using readOGR (file is named "roadtrl020"), and fortify it (using "STATE" as region); call the object "roads"
# YOUR CODE HERE
roads <- readOGR(dsn="shapefiles", layer="roadtrl020")
## OGR data source with driver: ESRI Shapefile 
## Source: "shapefiles", layer: "roadtrl020"
## with 1429 features
## It has 10 fields
# Add road data onto our plot:
ggplot(yosemite, aes(x=long, y=lat)) + 
  geom_polygon(fill="skyblue", color="black") +
  geom_point(data=coordinates, aes(x=long, y=lat)) +
  coord_fixed(ratio=1.3) +
  geom_path(data=roads, aes(x=long, y=lat, group=group), color="grey") # need group to retain order

# Zoom in (set x and y limits) and remove axes:
ggplot(yosemite, aes(x=long, y=lat)) + 
  geom_polygon(fill="darkslategray4", color="black") +
  geom_point(data=coordinates, aes(x=long, y=lat), alpha=0.3, color="coral4") +
  geom_path(data=roads, aes(x=long, y=lat, group=group), color="black") +
  xlim(-119.9,-119.1) + ylim(37.45, 38.25) +
  coord_fixed(ratio=1.3) +
  theme_nothing()