library(ggmap)
## Loading required package: ggplot2
library(tmap)
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-6, (SVN revision 773)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1
library(sp)
library(ggplot2)
library(raster)
library(rnaturalearth)
library(ggrepel)
Ok. Let’s start with the low-hanging fruit: using ggplot to make maps. The map_data() verb is ggplot verb that lets you pull in pre-defined map sets. I want to make a basic map of places I’ve lived, just to get the juices flowing, so I’ll grap a contiguous USA map and a map of State Boundaries.
usa <- map_data("usa")
states <- map_data("state")
I’ll have a peek, since I’ve learned this is a pretty good habit to get into.
head(states)
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
## 4 -87.53076 30.33239 1 4 alabama <NA>
## 5 -87.57087 30.32665 1 5 alabama <NA>
## 6 -87.58806 30.32665 1 6 alabama <NA>
I want to select a subset of the states I’ve lived in, so I grab the relevant ones using the subset() verb.
t_states <- subset(states, region %in% c("california", "oregon", "illinois"))
Then, I make a data frame with latitude and longitude information for the towns I’ve been affiliated with. Much like the state and country data, the plotting function from ggplot uses the longitude values on the x-axis (which includes negative values, as opposed to E/W), so it’s important to make sure they’re lined up properly. I would have made this map also include places I worked or have been more loosely affiliated with, to give a clearer picture of my history, but that meant more information in the data frame. Nope.
t_homes <- data.frame(
long = c(-122.2711, -123.0868, -122.6587, -120.1833, -87.6298, -118.2437),
lat = c(37.8044, 44.0521, 45.5122, 39.3280, 41.8781, 34.0522),
names = c("Oakland", "Eugene", "Portland", "Truckee", "Chicago", "Los Angeles"),
stringsAsFactors = FALSE
)
To the plot!
ggplot() +
geom_polygon(data = states, aes(x = long, y = lat, group = group), fill = "cornsilk2", color = "grey80") +
geom_polygon(data = t_states, aes (x =long, y = lat, group = group), fill = "darkseagreen3", color = "grey80") +
geom_point(data = t_homes, aes(x = long, y = lat), color = "lightcoral", size = 2) +
geom_text(data = t_homes, aes(x = long, y = lat, label = names), nudge_x = 2, nudge_y = .75) +
theme(panel.background = element_rect(fill = "lightcyan2", color = "black", size = 1),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
ggtitle("A few of the places I've lived...")
Probably best not to use cyan as the map background, because that implies water when there isn’t neceesarily any. But this is good enough for this exercise, I think.
And now for my next trick: I’d like to pull maps from Google or Stamen…
…But it turns out I can’t. I understand there have been issues with the Google API (or at least there are hoops to jump through, and I don’t feel like jumping). I thought I would try to grab it from Stamen, as shown in the video, but that didn’t seem to work (somehow it still wants to use Google as the source). And maybe I’m an idiot, but I couldn’t figure out how to download a shapefile from Stamen either…they seem to let you make images from their maps but not much else. So I’ll move on, I guess.
Let’s try shapefiles, downloaded from Natural Earth.
us_states <- readOGR("Shapefiles/ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/travisbott/Desktop/Work/2018-19/1 Fall/GEOG 208/Week 7/Intro_to_Mapping/Shapefiles/ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces.shp", layer: "ne_50m_admin_1_states_provinces"
## with 100 features
## It has 83 fields
## Integer64 fields read as strings: ne_id
populated_cites <- readOGR("Shapefiles/ne_50m_populated_places/ne_50m_populated_places.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/travisbott/Desktop/Work/2018-19/1 Fall/GEOG 208/Week 7/Intro_to_Mapping/Shapefiles/ne_50m_populated_places/ne_50m_populated_places.shp", layer: "ne_50m_populated_places"
## with 1249 features
## It has 119 fields
## Integer64 fields read as strings: wof_id ne_id
urban_areas <- readOGR("Shapefiles/ne_50m_urban_areas/ne_50m_urban_areas.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/travisbott/Desktop/Work/2018-19/1 Fall/GEOG 208/Week 7/Intro_to_Mapping/Shapefiles/ne_50m_urban_areas/ne_50m_urban_areas.shp", layer: "ne_50m_urban_areas"
## with 2143 features
## It has 4 fields
## Integer64 fields read as strings: scalerank
The states shapefile (which I started by naming us_states) is of the administrative units of the US, Brazil, and Australia. So I need to select by the US, and then drop (sorry, guys) Alaska and Hawaii, because they’re going to screw everything up.
us_only <- us_states$admin == "United States of America"
USA <- us_states[us_only, ]
USA_AKHI <- USA$name %in% c("Hawaii","Alaska")
USA_continental <- USA[!USA_AKHI, ]
And then I need to select the world cities by the same.
us_cities_only <- populated_cites$ADM0NAME == "United States of America"
USA_Cities <- populated_cites[us_cities_only, ]
The plot:
I know it’s supposed to use a lot of the same ‘grammar’ as ggplot, but not all of it’s the same, and there are certain arguments I’m having trouble figuring out. It’s the end of the night for me on Monday, so I’m giving up, but I need to figure out: how to re-scale the population numbers, how to re-name the legend (this one seemed impossible), and how to change the title position. There are arguments for the last one but they don’t seem to do anything! So I just left it off, because the default was in a stupid place. Also, Michigan gets the shaft.
tm_shape(USA_continental) +
tm_fill(col = "darkolivegreen3", legend.show = FALSE) +
tm_borders("white") +
tm_shape(USA_Cities) +
tm_dots(size = "POP_MAX", scale = 3, col = "coral1", alpha = .75) +
tm_layout(title = "Relative Populations of US Cities", title.position = 1)
## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
lakes <- map_data("lakes")
Honestly, I keep having trouble here. I’ve tried several times to get us.cities and canada.cities from the maps package but they don’t seem to work, so I’m giving up. I wanted to make a map of “Cities of the Great Lakes”, which I thought would be easiest by downlaoding a directory, but I’ll make another data frame. For you.
great_lakes_cites_df <- data.frame(
long = c(-87.7319, -92.1109, -89.2405, -87.9904, -84.4238, -80.0362, -79.0087, -81.7058, -83.0250, -84.6269),
lat = c(41.8340, 46.7651, 48.4026, 44.5229, 46.5317, 45.3432, 43.0995, 41.4976, 42.3294, 45.8657),
names = c("Chicago", "Duluth", "Thunder Bay", "Green Bay", "Sault Ste. Marie", "Parry Sound", "Niagra", "Cleveland", "Detroit", "Mackinac Island"),
stringsAsFactors = FALSE
)
head(lakes)
## long lat group order region subregion
## 1 17.97978 59.32905 1 1 Malaren <NA>
## 2 17.87617 59.27081 1 2 Malaren <NA>
## 3 17.57051 59.26763 1 3 Malaren <NA>
## 4 17.47451 59.29151 1 4 Malaren <NA>
## 5 17.37070 59.29492 1 5 Malaren <NA>
## 6 17.30459 59.27217 1 6 Malaren <NA>
Chop ’em up! Gimme them lakes.
great_lakes <- subset(lakes, region %in% c("Great Lakes"))
Mise en plot. Oui, Chef.
ggplot() +
geom_polygon(data = great_lakes, aes(x = long, y = lat, group = group), fill = "cyan4") +
geom_point(data = great_lakes_cites_df, aes(x = long, y = lat), color = "red", size = 3) +
theme(
panel.background = element_rect(fill = "cornsilk3", color = "black", size = 1),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
geom_text(data = great_lakes_cites_df, aes(x = long, y = lat, label = names, fontface = "bold"), color = "white", nudge_y = .2, nudge_x = .75) +
ggtitle("A Handful of Cities around the Great Lakes")
Sorry, ya’ll. This is what you get this week.