Space

Load and look at data

First we’ll load the libraries we’re going to use:

library(tidyverse)  # For ggplot, dplyr, and friends
library(sf)         # For GIS magic

Next we’ll load all the different shapefiles we downloaded using read_sf():

# Download "Admin 0 – Countries" from
# https://www.naturalearthdata.com/downloads/110m-cultural-vectors/
world_map <- read_sf("data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp")

# Download cb_2018_us_state_20m.zip under "States" from
# https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
us_states <- read_sf("data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp")

# Download cb_2018_us_county_5m.zip under "County" from
# https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
us_counties <- read_sf("data/cb_2018_us_county_5m/cb_2018_us_county_5m.shp")

# Download "Admin 1 – States, Provinces" from
# https://www.naturalearthdata.com/downloads/10m-cultural-vectors/
us_states_hires <- read_sf("data/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp")

# Download "Rivers + lake centerlines" from
# https://www.naturalearthdata.com/downloads/10m-physical-vectors/
rivers_global <- read_sf("data/ne_10m_rivers_lake_centerlines/ne_10m_rivers_lake_centerlines.shp")

# Download "Rivers + lake centerlines, North America supplement" from
# https://www.naturalearthdata.com/downloads/10m-physical-vectors/
rivers_na <- read_sf("data/ne_10m_rivers_north_america/ne_10m_rivers_north_america.shp")

# Download "Lakes + Reservoirs" from
# https://www.naturalearthdata.com/downloads/10m-physical-vectors/
lakes <- read_sf("data/ne_10m_lakes/ne_10m_lakes.shp")

# Download from https://data.georgiaspatial.org/index.asp?body=preview&dataId=41516
# after creating an account and logging in
ga_schools <- read_sf(file.path("data", "schools_2009", "DOE Schools 2009.shp"))

Basic plotting

If you look at the world_map dataset in RStudio, you’ll see it’s just a standard data frame with 177 rows and 95 columns. The last column is the magical geometry column with the latitude/longitude details for the borders for every country. RStudio only shows you 50 columns at a time in the RStudio viewer, so you’ll need to move to the next page of columns with the » button in the top left corner.

Because this is just a data frame, we can do all our normal dplyr things to it. Let’s get rid of Antarctica, since it takes up a big proportion of the southern hemisphere:

world_sans_antarctica <- world_map %>% 
  filter(ISO_A3 != "ATA")

Ready to plot a map? Here’s all you need to do:

ggplot() + 
  geom_sf(data = world_sans_antarctica)

If you couldn’t tell from the lecture, I’m completely blown away by how amazingly easy this every time I plot a map :)

Because this a regular ggplot geom, all our regular aesthetics and themes and everything work:

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  theme_void()

The Natural Earth dataset happens to come with some columns with a coloring scheme with 7–13 colors (MAPCOLOR7, MAPCOLOR9, etc.) so that no countries with a shared border share a color. We can fill by that column:

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          aes(fill = as.factor(MAPCOLOR7)),
          color = "#401D16", size = 0.25) +
  scale_fill_viridis_d(option = "plasma") +
  guides(fill = FALSE) +
  theme_void()

World map with different projections

Changing projections is trivial: add a coord_sf() layer where you specify the CRS you want to use.

Here’s Robinson (yay):

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = 54030) +  # Robinson
  # Or use the name instead of the number
  # coord_sf(crs = "+proj=robin")
  theme_void()

Here’s sinusoidal:

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = 54008) +  # Sinusoidal
  theme_void()

And here’s Mercator (ewww):

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = 3785) +  # Mercator
  # Or use the name instead of the number
  # coord_sf(crs = "+proj=merc")
  theme_void()

US map with different projections

This same process works for any shapefile. The map of the US can also be projected differently—two common projections are NAD83 and Albers. We’ll take the us_states dataset, remove Alaska, Hawaii, and Puerto Rico (they’re so far from the rest of the lower 48 states that they make an unusable map—if you want to include them, it’s easiest to plot them as separate plots and use patchwork to stitch them together), and plot it.

lower_48 <- us_states %>% 
  filter(!(NAME %in% c("Alaska", "Hawaii", "Puerto Rico")))

ggplot() + 
  geom_sf(data = lower_48, fill = "#192DA1", color = "white", size = 0.25) +
  coord_sf(crs = 4269) +  # NAD83
  theme_void()

ggplot() + 
  geom_sf(data = lower_48, fill = "#192DA1", color = "white", size = 0.25) +
  coord_sf(crs = 102003) +  # Albers
  theme_void()

Individual states

Again, because these shapefiles are really just fancy data frames, we can filter them with normal dplyr functions. Let’s plot just Georgia:

only_georgia <- lower_48 %>% 
  filter(NAME == "Georgia")

ggplot() +
  geom_sf(data = only_georgia, fill = "#EC8E55") +
  theme_void()

We can also use a different projection. If we look at spatialreference.org, there’s a version of NAD83 that’s focused specifically on Georgia.

ggplot() +
  geom_sf(data = only_georgia, fill = "#EC8E55") +
  theme_void() +
  coord_sf(crs = 2239)  # NAD83 focused on Georgia

There’s one small final issue though: we’re missing all the Atlantic islands in the southeast like Cumberland Island and Amelia Island. That’s because we’re using the Census’s low resolution (20m) data. That’s fine for the map of the whole country, but if we’re looking at a single state, we probably want better detail in the borders. We can use the Census’s high resolution (500k) data, but even then it doesn’t include the islands for whatever reason, but Natural Earth has high resolution US state data that does have the islands, so we can use that:

only_georgia_high <- us_states_hires %>% 
  filter(iso_3166_2 == "US-GA")

ggplot() +
  geom_sf(data = only_georgia_high, fill = "#EC8E55") +
  theme_void() +
  coord_sf(crs = 2239)  # NAD83 focused on Georgia

Perfect.

Plotting multiple shapefile layers

The state shapefiles from the Census Bureau only include state boundaries. If we want to see counties in Georgia, we need to download and load the Census’s county shapefiles (which we did above). We can then add a second geom_sf() layer for the counties.

First we need to filter the county data to only include Georgia counties. The counties data doesn’t include a column with the state name or state abbreviation, but it does include a column named STATEFP, which is the state FIPS code. Looking at lower_48 we can see that the state FIPS code for Georgia is 13, so we use that to filter.

ga_counties <- us_counties %>% 
  filter(STATEFP == 13)

Now we can plot just the counties:

ggplot() +
  geom_sf(data = ga_counties) +
  theme_void()

Technically we can just draw the county boundaries instead of layer the state boundary + the counties, since the borders of the counties make up the border of the state. But there’s an advantage to including both: we can use different aesthetics on each, like adding a thicker border on the state:

ggplot() +
  geom_sf(data = only_georgia_high, color = "#EC8E55", size = 3) +
  geom_sf(data = ga_counties, fill = "#A5D46A", color = "white") +
  theme_void()

It’s also useful if we want to only show some counties, like metropolitan Atlanta:

atl_counties <- ga_counties %>% 
  filter(NAME %in% c("Cherokee", "Clayton", "Cobb", "DeKalb", "Douglas",
                     "Fayette", "Fulton", "Gwinnett", "Henry", "Rockdale"))
ggplot() +
  geom_sf(data = only_georgia_high, fill = "#EC8E55") +
  geom_sf(data = atl_counties, fill = "#A5D46A", color = "white") +
  theme_void()

Plotting multiple shapefile layers when some are bigger than the parent shape

So far we’ve been able to filter out states and counties that we don’t want to plot using filter(), which works because the shapefiles have geometry data for each state or county. But what if you’re plotting stuff that doesn’t follow state or county boundaries, like freeways, roads, rivers, or lakes?

At the beginning we loaded a shapefile for all large and small rivers in the US. Look at the first few rows of rivers_na:

head(rivers_na)

## Simple feature collection with 6 features and 37 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: -100 ymin: 29 xmax: -86 ymax: 36
## CRS:            4326
## # A tibble: 6 x 38
##   featurecla scalerank rivernum dissolve name  name_alt note  name_full min_zoom strokeweig uident min_label label wikidataid name_ar name_bn name_de
##   <chr>          <dbl>    <dbl> <chr>    <chr> <chr>    <chr> <chr>        <dbl>      <dbl>  <dbl>     <dbl> <chr> <chr>      <chr>   <chr>   <chr>  
## 1 River             10    22360 22360Ri… Colo… <NA>     ID i… Colorado…      6         0.3  1.99e6       7   Colo… Q847785    <NA>    <NA>    Colora…
## 2 River             10    22572 22572Ri… Cima… <NA>     ID i… Cimarron…      6         0.25 2.15e6       7   Cima… Q1092055   <NA>    <NA>    Cimarr…
## 3 River             10    22519 22519Ri… Wash… <NA>     ID i… Washita …      6         0.25 1.95e6       7   Wash… Q2993598   <NA>    <NA>    Washita
## 4 River             10    22519 22519Ri… Wash… <NA>     ID i… Washita …      6         0.15 1.95e6       7   Wash… Q2993598   <NA>    <NA>    Washita
## 5 River             11    22422 22422Ri… Cone… <NA>     ID i… Conecuh …      6.7       0.15 2.17e6       7.7 Cone… Q5159475   <NA>    <NA>    <NA>   
## 6 River             10    22421 22421Ri… Pea   <NA>     ID i… Pea River      6         0.15 1.96e6       7   Pea   Q7157190   <NA>    <NA>    <NA>   
## # … with 21 more variables: name_en <chr>, name_es <chr>, name_fr <chr>, name_el <chr>, name_hi <chr>, name_hu <chr>, name_id <chr>, name_it <chr>,
## #   name_ja <chr>, name_ko <chr>, name_nl <chr>, name_pl <chr>, name_pt <chr>, name_ru <chr>, name_sv <chr>, name_tr <chr>, name_vi <chr>,
## #   name_zh <chr>, wdid_score <int>, ne_id <dbl>, geometry <MULTILINESTRING [°]>

The first row is the whole Colorado river, which flows through seven states. We can’t just use filter() to only select some parts of it based on states.

Here’s what happens if we combine our Georgia map with rivers and lakes:

ggplot() +
  geom_sf(data = only_georgia, fill = "#EC8E55") +
  geom_sf(data = rivers_na) +
  theme_void()

It plots Georgia, and it’s filled with orange, but it also plots every single river in North America. Oops.

We need to do a little GIS work to basically use only_georgia as a cookie cutter and keep only the rivers that are contained in the only_georgia boundaries. Fortunately, there’s a function in the sf package that does this: st_intersection(). Feed it two shapefile datasets and it will select the parts of the second that fall within the boundaries of the first:

ga_rivers_na <- st_intersection(only_georgia, rivers_na)

## Error in geos_op2_geom("intersection", x, y): st_crs(x) == st_crs(y) is not TRUE

Oh no! An error! It’s complaining that the reference systems used in these two datasets don’t match. We can check the CRS with st_crs():

st_crs(only_georgia)

## Coordinate Reference System:
##   User input: 4269 
##   wkt:
## GEOGCS["NAD83",
##     DATUM["North_American_Datum_1983",
##         SPHEROID["GRS 1980",6378137,298.257222101,
##             AUTHORITY["EPSG","7019"]],
##         TOWGS84[0,0,0,0,0,0,0],
##         AUTHORITY["EPSG","6269"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4269"]]

st_crs(rivers_na)

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

The Georgia map uses 4269 (or NAD83), while the rivers map uses 4326 (or the GPS system of latitude and longitude). We need to convert one of them to make them match. It doesn’t matter which one.

only_georgia_4326 <- only_georgia %>% 
  st_transform(crs = 4326)

ga_rivers_na <- st_intersection(only_georgia_4326, rivers_na)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all geometries

You’ll get an ominous warning, but you should be okay—it’s just because flattening globes into flat planes is hard, and the cutting might not be 100% accurate, but it’ll be close enough for our mapping purposes.

Now we can plot our state shape and the truncated rivers:

ggplot() +
  geom_sf(data = only_georgia, fill = "#EC8E55") +
  geom_sf(data = ga_rivers_na) +
  theme_void()

Hey! It worked! Let’s put all the rivers and lakes on at once and make it a little more artsy. We’ll use the high resolution Georgia map too, which conveniently already matches the CRS of the rivers and lakes:

ga_rivers_na <- st_intersection(only_georgia_high, rivers_na)
ga_rivers_global <- st_intersection(only_georgia_high, rivers_global)
ga_lakes <- st_intersection(only_georgia_high, lakes)

ggplot() +
  geom_sf(data = only_georgia_high, 
          color = "black", size = 0.1, fill = "black") +
  geom_sf(data = ga_rivers_global, size = 0.3, color = "grey80") +
  geom_sf(data = ga_rivers_na, size = 0.15, color = "grey80") +
  geom_sf(data = ga_lakes, size = 0.3, fill = "grey80", color = NA) +
  coord_sf(crs = 4269) +  # NAD83
  theme_void()

Heck yeah. That’s a great map. This is basically what Kieran Healy did here, but he used even more detailed shapefiles from the US Geological Survey.

Plotting schools in Georgia

Shapefiles are not limited to just lines and areas—they can also contain points. I made a free account at the Georgia GIS Clearinghouse, searched for “schools” and found a shapefile of all the K–12 schools in 2009. This is the direct link to the page, but it only works if you’re logged in to their system. This is the official metadata for the shapefile, which you can see if you’re not logged in, but you can’t download anything. It’s a dumb system and other states are a lot better at offering their GIS data (like, here’s a shapefile of all of Utah’s schools and libraries as of 2017, publicly accessible without an account).

We loaded the shapefile up at the top, but now let’s look at it:

head(ga_schools)

## Simple feature collection with 6 features and 16 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 2100000 ymin: 320000 xmax: 2200000 ymax: 5e+05
## proj4string:    +proj=tmerc +lat_0=30 +lon_0=-84.16666666666667 +k=0.9999 +x_0=700000.0000001107 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs 
## # A tibble: 6 x 17
##      ID  DATA COUNTY  DISTRICT  SCHOOLNAME  GRADES  ADDRESS CITY  STATE ZIP   TOTAL SCHOOLID DOE_CONGRE CONGRESS SENATE HOUSE                 geometry
##   <dbl> <dbl> <chr>   <chr>     <chr>       <chr>   <chr>   <chr> <chr> <chr> <dbl>    <dbl> <chr>      <chr>    <chr>  <chr> <POINT [US_survey_foot]>
## 1  4313   224 Early   Early Co… Early Coun… PK,KK,… 283 Ma… Blak… GA    3982…  1175    43549 2          002      011    149           (2052182 494322)
## 2  4321   227 Early   Early Co… ETN Eckerd… 06,07,… 313 E … Blak… GA    3982…    30    47559 2          002      011    149            (2053200 5e+05)
## 3  4329   226 Early   Early Co… Early Coun… 06,07,… 12053 … Blak… GA    3982…   539    43550 2          002      011    149            (2055712 5e+05)
## 4  4337   225 Early   Early Co… Early Coun… 09,10,… 12020 … Blak… GA    3982…   716    43552 2          002      011    149            (2055712 5e+05)
## 5  4345   189 Decatur Decatur … John Johns… PK,KK,… 1947 S… Bain… GA    3981…   555    43279 2          002      011    172           (2168090 321781)
## 6  4353   192 Decatur Decatur … Potter Str… PK,KK,… 725 Po… Bain… GA    3981…   432    43273 2          002      011    172           (2168751 327375)

We have a bunch of columns like GRADES that has a list of what grades are included in the school, and TOTAL, which I’m guessing is the number of students. Let’s map it!

If we add a geom_sf() layer just for ga_schools, it’ll plot a bunch of points:

ggplot() +
  geom_sf(data = ga_schools)

One of these rows is wildly miscoded and ended up Indonesia! If you sort by the geometry column in RStudio, you’ll see that it’s most likely Allatoona High School in Cobb County (id = 22097). The coordinates are different from all the others, and it has no congressional district information. Let’s remove it.

ga_schools_fixed <- ga_schools %>% 
  filter(ID != 22097)

ggplot() +
  geom_sf(data = ga_schools_fixed)

That’s better. However, all we’re plotting now are the points—we’ve lost the state and/or county boundaries. Let’s include those:

ggplot() +
  geom_sf(data = only_georgia_high) +
  geom_sf(data = ga_schools_fixed) +
  theme_void()

We’re getting closer. We have some issues with overplotting, so let’s shrink the points down and make them a little transparent:

ggplot() +
  geom_sf(data = only_georgia_high) +
  geom_sf(data = ga_schools_fixed, size = 0.5, alpha = 0.5) +
  theme_void()

Neat. One last thing we can do is map the TOTAL column to the color aesthetic and color the points by how many students attend each school:

ggplot() +
  geom_sf(data = only_georgia_high) +
  geom_sf(data = ga_schools_fixed, aes(color = TOTAL), size = 0.75, alpha = 0.5) +
  scale_color_viridis_c() +
  theme_void()

Most schools appear to be under 1,000 students, except for a cluster in Gwinnett County north of Atlanta. Its high schools have nearly 4,000 students each!

ga_schools_fixed %>% 
  select(COUNTY, SCHOOLNAME, TOTAL) %>% 
  arrange(desc(TOTAL)) %>% 
  head()

## Simple feature collection with 6 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 2300000 ymin: 1400000 xmax: 2400000 ymax: 1500000
## proj4string:    +proj=tmerc +lat_0=30 +lon_0=-84.16666666666667 +k=0.9999 +x_0=700000.0000001107 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs 
## # A tibble: 6 x 4
##   COUNTY   SCHOOLNAME                  TOTAL                 geometry
##   <chr>    <chr>                       <dbl> <POINT [US_survey_foot]>
## 1 Gwinnett Mill Creek High School       3997        (2384674 1482772)
## 2 Gwinnett Collins Hill High School     3720        (2341010 1461730)
## 3 Gwinnett Brookwood High School        3455        (2334543 1413396)
## 4 Gwinnett Grayson High School          3230        (2370186 1408579)
## 5 Gwinnett Peachtree Ridge High School  3118        (2319344 1459458)
## 6 Gwinnett Berkmar High School          3095        (2312983 1421933)

Making your own geoencoded data

So, plotting shapefiles with geom_sf() is magical because sf deals with all of the projection issues for us automatically and it figures out how to plot all the latitude and longitude data for us automatically. But lots of data doesn’t some as shapefiles. The rats data from mini project 1, for instance, has two columns indicating the latitude and longitude of each rat sighting, but those are stored as just numbers. If we try to use geom_sf() with the rat data, it won’t work. We need that magical geometry column.

Fortunately, if we have latitude and longitude information, we can make our own geometry column.

Let’s say we want to mark some cities on our map of Georgia. We can make a mini dataset using tribble(). I found these points from Google Maps: right click anywhere in Google Maps, select “What’s here?”, and you’ll see the exact coordinates for that spot.

ga_cities <- tribble(
  ~city, ~lat, ~long,
  "Atlanta", 33.748955, -84.388099,
  "Athens", 33.950794, -83.358884,
  "Savannah", 32.113192, -81.089350
)
ga_cities

## # A tibble: 3 x 3
##   city       lat  long
##   <chr>    <dbl> <dbl>
## 1 Atlanta   33.7 -84.4
## 2 Athens    34.0 -83.4
## 3 Savannah  32.1 -81.1

This is just a normal dataset, and the lat and long columns are just numbers. R doesn’t know that those are actually geographic coordinates. This is similar to the rats data, or any other data that has columns for latitude and longitude.

We can convert those two columns to the magic geometry column with the st_as_sf() function. We have to define two things in the function: which coordinates are the longitude and latitude, and what CRS the coordinates are using. Google Maps uses 4326, or the GPS system, so we specify that:

ga_cities_geometry <- ga_cities %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)
ga_cities_geometry

## Simple feature collection with 3 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -84 ymin: 32 xmax: -81 ymax: 34
## CRS:            EPSG:4326
## # A tibble: 3 x 2
##   city        geometry
##   <chr>    <POINT [°]>
## 1 Atlanta     (-84 34)
## 2 Athens      (-83 34)
## 3 Savannah    (-81 32)

The longitude and latitude columns are gone now, and we have a single magical geometry column. That means we can plot it with geom_sf():

ggplot() +
  geom_sf(data = only_georgia_high, fill = "#EC8E55") +
  geom_sf(data = ga_cities_geometry, size = 3) +
  theme_void()

We can use geom_sf_label() (or geom_sf_text()) to add labels in the correct locations too. It will give you a warning, but you can ignore it—again, it’s complaining that the positioning might not be 100% accurate because of issues related to taking a globe and flattening it. It’s fine.

ggplot() +
  geom_sf(data = only_georgia_high, fill = "#EC8E55") +
  geom_sf(data = ga_cities_geometry, size = 3) +
  geom_sf_label(data = ga_cities_geometry, aes(label = city),
                nudge_y = 0.2) +
  theme_void()

Plotting other data on maps

So far we’ve just plotted whatever data the shapefile creators decided to include and publish in their data. But what if you want to visualize some other variable on a map? We can do this by combining our shapefile data with any other kind of data, as long as the two have a shared column. For instance, we can make a choropleth map of life expectancy with data from the World Bank.

First, let’s grab some data from the World Bank for just 2015:

library(WDI)  # For getting data from the World Bank

indicators <- c("SP.DYN.LE00.IN")  # Life expectancy

wdi_raw <- WDI(country = "all", indicators, extra = TRUE, 
               start = 2015, end = 2015)

Let’s see what we got:

head(wdi_raw)

## # A tibble: 6 x 11
##   iso2c country                                SP.DYN.LE00.IN  year iso3c region             capital         longitude latitude income    lending     
##   <chr> <chr>                                           <dbl> <dbl> <chr> <chr>              <chr>               <dbl>    <dbl> <chr>     <chr>       
## 1 1A    Arab World                                       71.2  2015 ARB   Aggregates         <NA>                NA        NA   Aggregat… Aggregates  
## 2 1W    World                                            71.9  2015 WLD   Aggregates         <NA>                NA        NA   Aggregat… Aggregates  
## 3 4E    East Asia & Pacific (excluding high i…           74.5  2015 EAP   Aggregates         <NA>                NA        NA   Aggregat… Aggregates  
## 4 7E    Europe & Central Asia (excluding high…           72.7  2015 ECA   Aggregates         <NA>                NA        NA   Aggregat… Aggregates  
## 5 8S    South Asia                                       68.6  2015 SAS   Aggregates         <NA>                NA        NA   Aggregat… Aggregates  
## 6 AD    Andorra                                          NA    2015 AND   Europe & Central … Andorra la Vel…      1.52     42.5 High inc… Not classif…

We have a bunch of columns here, but we care about two in particular: life expectancy, and the ISO3 code. This three-letter code is a standard system for identifying countries (see the full list here), and that column will let us combine this World Bank data with the global shapefile, which also has a column for the ISO3 code.

(We also have columns for the latitude and longitude for each capital, so we could theoretically convert those to a geometry column with st_as_sf() and plot world capitals, which would be neat, but we won’t do that now.)

Let’s clean up the WDI data and shrink it down substantially:

wdi_clean_small <- wdi_raw %>% 
  select(life_expectancy = SP.DYN.LE00.IN, iso3c)
wdi_clean_small

## # A tibble: 264 x 2
##    life_expectancy iso3c
##              <dbl> <chr>
##  1            71.2 ARB  
##  2            71.9 WLD  
##  3            74.5 EAP  
##  4            72.7 ECA  
##  5            68.6 SAS  
##  6            NA   AND  
##  7            77.3 ARE  
##  8            63.4 AFG  
##  9            76.5 ATG  
## 10            78.0 ALB  
## # … with 254 more rows

Next we need to merge this tiny dataset into the world_map_sans_antarctica shapefile data we were using earlier. To do this we’ll use a function named left_join(). We feed two data frames into left_join(), and R will keep all the rows from the first and include all the columns from both the first and the second wherever the two datasets match with one specific column. That’s wordy and weird—stare at this animation here for a few seconds to see what’s really going to happen. We’re essentially going to append the World Bank data to the end of the world shapefiles and line up rows that have matching ISO3 codes. The ISO3 column is named ISO_A3 in the shapefile data, and it’s named iso3c in the WDI data, so we tell left_join() that those are the same column:

world_map_with_life_expectancy <- world_sans_antarctica %>% 
  left_join(wdi_clean_small, by = c("ISO_A3" = "iso3c"))

If you look at this dataset in RStudio now and look at the last column, you’ll see the WDI life expectancy right next to the magic geometry column.

We technically didn’t need to shrink the WDI data down to just two columns—had we left everything else, all the WDI columns would have come over to the world_sans_antarctica, including columns for region and income level, etc. But I generally find it easier and cleaner to only merge in the columns I care about instead of making massive datasets with a billion extra columns.

Now that we have a column for life expectancy, we can map it to the fill aesthetic and fill each country by 2015 life expectancy:

ggplot() + 
  geom_sf(data = world_map_with_life_expectancy, 
          aes(fill = life_expectancy),
          size = 0.25) +
  coord_sf(crs = 54030) +  # Robinson
  scale_fill_viridis_c(option = "viridis") +
  labs(fill = "Life expectancy") +
  theme_void() +
  theme(legend.position = "bottom")

Voila! Global life expectancy in 2015!

(Sharp-eyed readers will notice that France and Norway are grayed out because they’re missing data. That’s because the ISO_A3 code in the Natural Earth data is missing for both France and Norway for whatever reason, so the WDI data didn’t merge with those rows. To fix that, we can do some manual recoding before joining in the WDI data)

world_sans_antarctica_fixed <- world_sans_antarctica %>% 
  mutate(ISO_A3 = case_when(
    # If the country name is Norway or France, redo the ISO3 code
    ADMIN == "Norway" ~ "NOR",
    ADMIN == "France" ~ "FRA",
    # Otherwise use the existing ISO3 code
    TRUE ~ ISO_A3
  )) %>% 
  left_join(wdi_clean_small, by = c("ISO_A3" = "iso3c"))

ggplot() + 
  geom_sf(data = world_sans_antarctica_fixed, 
          aes(fill = life_expectancy),
          size = 0.25) +
  coord_sf(crs = 54030) +  # Robinson
  scale_fill_viridis_c(option = "viridis") +
  labs(fill = "Life expectancy") +
  theme_void() +
  theme(legend.position = "bottom")