Mapping in R tutorial

Set up

Load packages

#clear environment
rm(list = ls()) 

# #To install packages:
# install.packages("package_name")

#load packages
library(sf) #spatial analysis
library(dplyr) #tidy data manipulation
library(ggplot2) #tidy figures

Get some data (source)
We will be looking at locations of hospitals in the city of Philadelphia.

url <-'https://opendata.arcgis.com/datasets/df8dc18412494e5abbb021e2f33057b2_0.zip'

#download zipfile and unzip
temp = tempfile()
temp2 = tempfile()
download.file(url, temp, mode='wb')
unzip(zipfile=temp, exdir=temp2)

# Read the shapefiles
sf.hospitals <- read_sf(temp2)

What are we working with?

head(sf.hospitals)

## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -75.2154 ymin: 39.94977 xmax: -74.98149 ymax: 40.07901
## Geodetic CRS:  WGS 84
## # A tibble: 6 x 9
##   OBJECTID HOSPITAL_N      STREET_ADD CITY  STATE ZIP_CODE PHONE_NUMB HOSPITAL_T
##      <int> <chr>           <chr>      <chr> <chr>    <int> <chr>      <chr>     
## 1        1 Aria Health- F~ 4900 Fran~ Phil~ PA       19124 215-831-2~ General m~
## 2        2 Aria Health- T~ 10800 Kni~ Phil~ PA       19114 215-612-4~ General m~
## 3        3 Belmont Center~ 4200 Monu~ Phil~ PA       19131 215-877-2~ Behaviora~
## 4        4 Chestnut Hill ~ 8835 Germ~ Phil~ PA       19118 215-248-8~ General m~
## 5        5 The Children's~ 34th Stre~ Phil~ PA       19104 215-590-1~ General m~
## 6        6 Einstein Medic~ 5501 Old ~ Phil~ PA       19141 215-456-7~ General m~
## # ... with 1 more variable: geometry <POINT [°]>

ggplot()+
  theme_classic()+
  geom_sf(data=sf.hospitals)

We have locations of hospitals in Philadelphia with various attributes about the hospital (street address, hospital type, zip code etc.) The visualization above is simply the coordinates plotted without any background or references at all. Even to those very familiar with Philadelphia, this visualization is lacking.

Method 1: Create maps by layering shapefiles

The first method of adding a background to our hospital points, is to add another layer of a relevant shapefile. One of the advantages of this method is that you can do further analysis such as creating choropleth maps (see my previous tutorial on choropleth maps). A potential disadvantage is that you don’t always have access to the shapefile of interest or may want a more detailed backaground.

To apply this method in the Philadelphia hospital example, we could obtain shapefiles of Philly zip codes and append the hospital information to plot:

#####################################################################
#Retrieve zip code tract shapefiles
######################################################################
url <- "https://opendata.arcgis.com/api/v3/datasets/b54ec5210cee41c3a884c9086f7af1be_0/downloads/data?format=shp&spatialRefId=4326"

#download zipfile and unzip
temp = tempfile()
temp2 = tempfile()
download.file(url, temp, mode='wb')
unzip(zipfile=temp, exdir=temp2)

# Read the shapefiles
sf.zip <- read_sf(temp2)

######################################################################
#Append hospital counts to zip code
######################################################################
HospCount <- plyr::count(sf.hospitals$ZIP_CODE) %>%
  rename("CODE" = "x", "HospCount"="freq") %>%
  mutate(CODE = as.character(CODE))

sf.zip <- left_join(sf.zip, HospCount)%>%
  mutate(HospCount = ifelse(is.na(HospCount), 0, HospCount))

######################################################################
#Visualize
######################################################################
map1 <- ggplot()+
  theme_void()+
  geom_sf(data=sf.zip, fill=NA, color="gray", size=0.1)+
  geom_sf(data=sf.hospitals)+
  ggtitle("Hospitals w/ zip background")


Pal <- RColorBrewer::brewer.pal(length(c(0:max(sf.zip$HospCount))), "Blues")
map2 <- ggplot()+
  theme_void()+
  geom_sf(data=sf.zip, aes(fill=as.factor(HospCount)), color="gray", size=0.1)+
  ggtitle("Zip codes colored by hospital")+
  scale_fill_manual(values = Pal, name="Hospital Count") 

map1;map2

Method 2: Add a background of imported tiles

There are several sources of map data that can be pulled in to form a background for your maps. My favorite source is OpenStreetMap data via stamen maps.

Pull in some background tiles from open street maps

library(ggmap)

## Zoom- the higher the zoom the more detail. Zoom=10 is pretty fast less detail. Zoom =18 is max detail but takes a couple minutes to load

#Get a 'bounding box' for philly
coords <- as.data.frame(st_coordinates(sf.zip)) #pull out lat long coords
bbox_PHL = c(left=min(coords$X),right=max(coords$X), #longitude
             top=max(coords$Y),bottom=min(coords$Y)) #latitude

#Open street maps
OSM1 <- ggmap(get_stamenmap(maptype='terrain', 
                            bbox=bbox_PHL,
                            zoom=10)) 

OSM2 <- ggmap(get_stamenmap(maptype='toner-lines', 
                            bbox=bbox_PHL,
                            zoom=10))

OSM3 <-  ggmap(get_stamenmap(maptype='watercolor', 
                            bbox=bbox_PHL,
                            zoom=10))


map4 <- OSM3 + 
  theme_void()+
  geom_sf(data=sf.hospitals, inherit.aes = FALSE)+
  ggtitle("stamen map 'water color'")

map5 <- OSM2 + 
  theme_void()+
  geom_sf(data=sf.zip,  
          fill= "dodgerblue3", color="dodgerblue3", size=.1, alpha=0.1,
          inherit.aes = FALSE)+
  geom_sf(data=sf.hospitals, color="red3", pch=17, size=2, inherit.aes = FALSE)+
  ggtitle("stamen map 'toner'")

map6 <- OSM1 + 
  theme_void()+
  geom_sf(data=sf.zip, aes(fill=as.factor(HospCount)), size=0.1, alpha=0.7,
          inherit.aes = FALSE)+
  scale_fill_manual(values = Pal) +
  theme(legend.position="none")+ 
  ggtitle("stamen map 'terrain'")


ggpubr::ggarrange(map4, map5, map6, ncol=3)

A couple of ‘good practices’ when mapping:

Include a scale bar
Include a North arrow

map5 +
  ggsn::north(sf.zip, location = "topleft")+
  ggsn::scalebar(sf.zip, 
                 dist=2, dist_unit="mi", transform=TRUE, 
                 anchor=c(x=mean(coords$X) + .2, y=mean(coords$Y))-.12, 
                 location = "topleft",
                 height=0.03, st.dist=0.05,
                 st.color="goldenrod3", box.fill=c("goldenrod3", "white"),
                 st.bottom=FALSE)

Method 3: Interactive maps with leaflet

One powerful package for mapping in R is leaflet. Leaflet works similar to the get_map function above where you call a base map and plot over it. However, these maps can be made to be interactive where you can zoom in and out and click on spatial objects to pop up information about different features.

library(leaflet)

lat_center <- mean(coords$Y, na.rm=TRUE)
long_center <- mean(coords$X, na.rm=TRUE)

#Set a color palette
pal1 <- colorFactor(c("red3", "goldenrod3", "dodgerblue3", "slateblue4"),
                        sf.hospitals$HOSPITAL_T) 

leaflet() %>%
    setView(lng = long_center, lat = lat_center, zoom = 11)%>% #center in AQP
    addProviderTiles(providers$CartoDB.Positron)%>% # add third party provider tile
    addCircleMarkers(data=sf.hospitals, 
                     popup = ~paste(HOSPITAL_N),
                     color=~pal1(HOSPITAL_T))%>%
    addLegend("bottomright", data= sf.hospitals,
              pal=pal1,
              values=~HOSPITAL_T,    
              title = 'Hospital type',
              opacity = 0.8) %>%
  addScaleBar()

Click on the hospital of interest to pull up the name of the hospital.

Mapping in R tutorial

B. Raynor

6/7/2022

Set up

Method 1: Create maps by layering shapefiles

Method 2: Add a background of imported tiles

Method 3: Interactive maps with leaflet