Mapping data with tmap

Reading the .shp file

In this exercise, we will use data from Chicago-IL communities. The source is the GeoDa Data and Lab. The dataset Airbnb (same information used in HW #1) refers to Airbnb rentals, socioeconomic indicators, and crime by community area in Chicago. You can find more information about it here. Here, you can download the related shapefile.

Again, the first step is to set up your working directory using setwd(), placing all the necessary data at the same place.

setwd("C:/Users/User/Desktop/Rlabs/Lab3")

After that, you need to require() all the packages you need to perform the task. If it is the first time you are using them, you need to install.packages() first.

libs <- c("tigris","rgdal", "maptools", "reshape", "sp", "spdep", "GISTools", "ggplot2", "tidyverse", "tmap", "viridis")
## If you do not have those packages, run 
## lapply(libs, install.packages, character.only = TRUE)
## before this
lapply(libs, require, character.only = TRUE)

To read the .shp file airbnb_Chicago 2015, use the function readOGR().

chitown<-readOGR("airbnb_Chicago 2015.shp", layer="airbnb_Chicago 2015")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\User\Desktop\Rlabs\Lab3\airbnb_Chicago 2015.shp", layer: "airbnb_Chicago 2015"
## with 77 features
## It has 20 fields
## Integer64 fields read as strings:  AREAID num_spots income_pc harship_in num_crimes num_theft population

Let’s check the file first:

class(chitown)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

That means your spatial object has some data attached to it. To see how the shape looks like:

plot(chitown)

To check the data attached to the shapefile (try View(chitown@data)),

head(chitown@data)

##         community    shape_area     shape_len AREAID response_r accept_r rev_rating
## 0         DOUGLAS 46004621.1581 31027.0545098     35   98.77143 94.51429   87.77778
## 1         OAKLAND 16913961.0408 19565.5061533     36   99.20000 90.10526   88.81250
## 2     FULLER PARK 19916704.8692          <NA>     37   68.00000       NA   91.75000
## 3 GRAND BOULEVARD 48492503.1554 28196.8371573     38   94.03704 83.61539   92.75000
## 4         KENWOOD 29071741.9283 23325.1679062     39   92.54286 88.14286   90.65625
## 5  LINCOLN SQUARE 71352328.2399 36624.6030848      4   90.00990 87.89691   95.24176
##    price_pp room_type num_spots poverty crowded dependency without_hs unemployed
## 0  78.15789  1.789474        38    29.6     1.8       30.7       14.3       18.2
## 1  53.77500  1.850000        20    39.7     1.3       40.4       18.4       28.7
## 2  84.00000  1.833333         6    51.2     3.2       44.9       26.6       33.9
## 3 119.53333  1.533333        30    29.3     3.3       39.5       15.9       24.3
## 4  77.99145  1.615385        39    21.7     2.4       35.4       11.3       15.7
## 5  80.18826  1.436364       110    10.9     3.4       25.5       13.4        8.2
##   income_pc harship_in num_crimes num_theft population
## 0     23791         47       5013      1241      18238
## 1     19252         78       1306       311       5918
## 2     10432         97       1764       383       2876
## 3     23472         57       6416      1428      21929
## 4     35911         26       2713       654      17841
## 5     37524         17       3670      1014      39493

Mapping data with tmap

Using the tmap package, one can easily make static quantile maps. In this one, we see the distribution of households below poverty (in percent) in Chicago communities:

tm_shape(chitown)+tm_borders(lwd = 2, col = "white", alpha = .4)+
  tm_polygons("poverty", style="quantile", alpha=0.7,palette = "-viridis", title="Poverty in Chicago Communities")+
  tm_legend(legend.title.size = 1.7,legend.text.size = 1, legend.position = c("left", "bottom"))

Besides, maps can be viewed interactively by switching to view mode, using the command tmap_mode("view")

tmap_mode("view")

tm_shape(chitown)+tm_borders(lwd = 2, col = "white", alpha = .4)+
  tm_polygons("poverty", style="quantile", alpha=0.7,palette = "-viridis", title="Poverty in Chicago Communities")+
  tm_legend(legend.title.size = 1.7,legend.text.size = 1, legend.position = c("left", "bottom"))

It is also possible to show more information in your map by choosing variables using popu.vars. Finally, if you want to change the colors, choose a different palette (here you have some options)

tm_shape(chitown)+tm_borders(lwd = 2, col = "white", alpha = .4)+
  tm_polygons("poverty", style="quantile", alpha=0.7,palette = "YlGnBu", 
              title="Poverty in Chicago Communities",id="community",
              popup.vars=c("Airbnb Prices"="price_pp","Number of Airbnb spots"="num_spots",
                           "Total Population"="population"))+
  tm_legend(legend.title.size = 1.7,legend.text.size = 1, legend.position = c("left", "bottom"))+
 tm_basemap(server="OpenStreetMap",alpha=0.5)

Tmap + Tidycensus

Using tidycensus with the option geometry = T, one can get coordinates of the desired areas together with Census information. Let’s try the commuting times per counties in the U.S. for 2009-2013, from the American Community Survey - 5 year. First, check the dataset acs5/profile.

require(tidycensus)
var<-load_variables(2013, dataset = "acs5/profile")
head(var)

## # A tibble: 6 x 3
##   name     label                                            concept                       
##   <chr>    <chr>                                            <chr>                         
## 1 DP02_00~ Estimate!!HOUSEHOLDS BY TYPE!!Total households   SELECTED SOCIAL CHARACTERISTI~
## 2 DP02_00~ Percent!!HOUSEHOLDS BY TYPE!!Total households    SELECTED SOCIAL CHARACTERISTI~
## 3 DP02_00~ Estimate!!HOUSEHOLDS BY TYPE!!Total households!~ SELECTED SOCIAL CHARACTERISTI~
## 4 DP02_00~ Percent!!HOUSEHOLDS BY TYPE!!Total households!!~ SELECTED SOCIAL CHARACTERISTI~
## 5 DP02_00~ Estimate!!HOUSEHOLDS BY TYPE!!Total households!~ SELECTED SOCIAL CHARACTERISTI~
## 6 DP02_00~ Percent!!HOUSEHOLDS BY TYPE!!Total households!!~ SELECTED SOCIAL CHARACTERISTI~

You are looking for COMMUTING TO WORK and mean travel time. The variable that corresponds to the average travel time to work is DP03_0025. Load it with geometry=T to get the coordinates for counties. Also, use shift_geo = T to place Alaska at the bottom of the map.

workers_time<-get_acs(geography = "county", 
                      variables = c(average_time="DP03_0025"), 
                      sumfile="acs5/profile",
                      year=2013,
                      geometry = T,
                      shift_geo = T,
                      output="wide")

Now, use tmap to plot average_timeE in a quantile map:

tmap_mode("plot")
tm_shape(workers_time)+tm_borders(lwd = 2, col = "white", alpha = .4)+
tm_polygons("average_timeE", style="quantile",palette = "BuPu",
              title="Commuting Time")+
tm_legend(legend.title.size = 1,legend.text.size = .8, 
              legend.position =c("left", "top"))+
tm_layout(title="Average Commuting Time in U.S. Counties", 
            title.position = c("center", "top"),
            title.size = 1.1)

Finally, you can adjust the margins inside the frame and have more space between the legend and the map. Just add one more line to calibrate the inner.margins. The order is bottom, left, top and right. The default value is .02, and I am increasing the space on the left and at the top.

tm_shape(workers_time)+tm_borders(lwd = 2, col = "white", alpha = .4)+
tm_polygons("average_timeE", style="quantile",palette = "YlGnBu",
              title="Commuting Time")+
tm_legend(legend.title.size = 1.2,legend.text.size = 1, 
              legend.position =c("left", "top"))+
tm_layout(title="Average Commuting Time in U.S. Counties", 
            title.position = c("center", "top"),
            title.size = 1.1)+
tm_layout(inner.margins = c(0.02, 0.15, 0.10, 0.02))