Country Level Data

Load some United States demographic data.

usa <- map_data("usa")

Before you look at it by running usa, Census data sets are typically huge for obvious reasons. Let’s see how big the file actually is.

dim(usa)

## [1] 7243    6

With 43,458 entries, that’s pretty big! Maybe we shouldn’t load it up (even in View() it takes awhile). Let’s instead just see the first and last six rows to get an idea of what it looks like to get an idea.

head(usa) 

##        long      lat group order region subregion
## 1 -101.4078 29.74224     1     1   main      <NA>
## 2 -101.3906 29.74224     1     2   main      <NA>
## 3 -101.3620 29.65056     1     3   main      <NA>
## 4 -101.3505 29.63911     1     4   main      <NA>
## 5 -101.3219 29.63338     1     5   main      <NA>
## 6 -101.3047 29.64484     1     6   main      <NA>

tail(usa) 

##           long      lat group order         region subregion
## 7247 -122.6187 48.37482    10  7247 whidbey island      <NA>
## 7248 -122.6359 48.35764    10  7248 whidbey island      <NA>
## 7249 -122.6703 48.31180    10  7249 whidbey island      <NA>
## 7250 -122.7218 48.23732    10  7250 whidbey island      <NA>
## 7251 -122.7104 48.21440    10  7251 whidbey island      <NA>
## 7252 -122.6703 48.17429    10  7252 whidbey island      <NA>

We have columns that can be used! The most obvious choice for map data is using the coordinates for longitude and latitude. So let’s plot using the longitudinal (long) and latitudinal (lat) columns.

ggplot() + 
  geom_polygon(data = usa, 
               aes(x=long, y = lat, group = group))

Well that looks terrible. So what’s happening here? In a nutshell, it’s plotting exactly how you would see the US if you flattened (aka projected) the coordinates from a sphere to a flat surface. Most maps account for this using a correction via an aspect ratio. In R, we use a command called coord_fixed

ggplot() + 
  geom_polygon(data = usa, 
               aes(x=long, y = lat, group = group)) + 
  coord_fixed(1.3) # 1.3 is a decent standard estimate

There we go. That looks better but of course its just a grey silhouette. Let’s first try adding some base color.

ggplot() + 
  geom_polygon(data = usa, aes(x=long, y = lat, group = group), color = "#1b85b8", 
               fill = "#559e83") + 
  coord_fixed(1.3)

Its looking better but something is missing. Hmmm….Oh!…Its the “S” part of “USA” - To get that, we need to look at state level data.

State Level Data

We have two choices: we can either bring in a particular state or bring them all in and then filter later. Unless you specifically know which state you want, its a good rule of thumb to bring in all state level data.

state <- map_data("state")

As a side note, if you know which state you want to look at, then just do the following. I’ll use California as an example.

California <- map_data("state", region = "CA")

A good thing to note is that region = is a bit misleading. At first you may think that it is a state’s abbreviation but it is not. In fact, you can’t even use abbreviations in that field. For example, if we use the abbreviation for Kansas, you’ll get nonsense

Kansas <- map_data("state", region = "KS")

## Error in map.poly(database, regions, exact, xlim, ylim, boundary, interior, : no recognized region names

(Yes the output is going off the page and no I can’t do anything about it)

The command region = is actually a filter so you need to put as much information as you can to ensure that you only get the state you want. For example, let’s see what happens when I do the following

Michigan <- map_data("state", region = "MI")

Well MI is certainly the standard abbreviation for Michigan but again, that is not what the command is doing. Other states have those letters:

unique(Michigan$region)

## [1] "michigan"    "minnesota"   "mississippi" "missouri"

That is not what we wanted so please make sure you use the region filter correctly. In a worst case scenario or if you don’t mind the extra effort, just type in the name of the state.

Michigan <- map_data("state", 
                     region = "Michigan")

unique(Michigan$region)

## [1] "michigan"

That’s much better! OK back to mapping. I’m actually going to stick with Michigan. Notice here that I am using the entire state data set and filtering Michigan out with the command subset and the == (double equals sign).

mi <- subset(state, 
             region=="michigan")

For a basic plot, we do the following:

ggplot() + 
  geom_polygon(data = mi, 
                        aes(x=long, y = lat, group = group)) + 
  coord_fixed(1.3)

Michigan’s shape is pretty unique and that sure looks like it but the plot is pretty bland. Let’s give it some color and a legend.

ggplot() + 
  geom_polygon(data = mi, aes(x=long, y = lat, group = group, fill=region)) +
  scale_fill_manual(values=c("#6A0032")) +
  coord_fixed(1.3)

Not bad but we can get even more granular. Let’s look at county level data.

County Level Data

Much like before, we’re going to pull in all of the county level data,

county <- map_data("county")

then we’ll filter it for Michigan,

mi_count <- subset(county, region=="michigan")

and finally we’ll plot it.

ggplot() + 
  geom_polygon(data = mi_count, aes(x=long, y = lat, group = group, fill=subregion), 
               color = "#FFFFFF") +
  scale_fill_viridis(discrete = TRUE, 
                     alpha=0.7, 
                     option="inferno") +
  theme_bw() +
  theme(plot.title = element_text(size = 26, color ="#105456",  vjust = -1),
        legend.position ="bottom",
        legend.direction = 'vertical',
        legend.text = element_text(size=15, color = "#59595B"),
        legend.title = element_blank(),
        legend.title.align = 0.5,
        legend.spacing.x = unit(0.5, "cm"), 
        legend.spacing.y = unit(0.5, "cm"),
        legend.background = element_rect(linetype = 0, size = 0.5, colour = 1),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  coord_fixed(1.3)

Your output may look a little different depending on your system settings.

That legend is taking up way too much space AND we don’t need the counties color coded by their names. Let’s get rid of it with guide = FALSE within scale_fill_viridis to get rid of that particular one or to turn all fill legends off, you can use guides(fill=FALSE) which is switch off in the plot above. In the example below, let’s try the one with guide = FALSE.

ggplot() + 
  geom_polygon(data = mi_count, aes(x=long, y = lat, group = group, fill=subregion), 
               color = "#FFFFFF") +
  scale_fill_viridis(discrete = TRUE, 
                     alpha=0.7, 
                     option="inferno", 
                     guide = FALSE) + #HERE IT IS
  theme_bw() +
  theme(plot.title = element_text(size = 26, color ="#105456",  vjust = -1),
        legend.position ="bottom",
        legend.direction = 'vertical',
        legend.text = element_text(size=15, color = "#59595B"),
        legend.title = element_blank(),
        legend.title.align = 0.5,
        legend.spacing.x = unit(0.5, "cm"), 
        legend.spacing.y = unit(0.5, "cm"),
        legend.background = element_rect(linetype = 0, size = 0.5, colour = 1),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  coord_fixed(1.3)

Well now we definitely have county level data for Michigan.

Please take a bit to review and play around with every layer. Many of the typical things you would want to change or amend about a plot are included.

OK now we should probably do something with this map. Let get some Census data.

Gathering and Making Important Decisions

Choose an endpoint aka a set of data that you want.

At this point it is completely understandable if you are wondering what an API is. Without getting overtly technical, an API is the acronym for Application Programming Interface. This is essentially a text based software that allows two applications to talk to each other. Why text based? Well its really the only aspect that all computers share. In fact, each time you use an app like YouTube, send a text message, or check the weather on your phone, etc. from any device worldwide, you’re using an API. As an example, take a look at the Twitter API https://developer.twitter.com/en/docs/tweets/search/api-reference.html that we’ll be using later in the term to perform some network data visualizations.

Now when looking at our endpoint, what should we choose? As of 2018, there are over 200 Census API endpoints are available, including the Decennial Census, American Community Survey, Poverty Statistics, and Population Estimates. You can find a list here: https://api.census.gov/data.html. Please look them over!

The censusapi package is designed to let you get data from all of those APIs using the same main function getCensus and the same syntax for each data set. In this example case, let’s look at the 2017 Density Estimates (variables can be found here: https://api.census.gov/data/2017/pep/population/variables.html) for all states within the US. Alternatively you can get that information directly through the API we have been using:

apis <- listCensusApis()

View(apis)

Now certain repositories require certain basic variables to be run (for example if years are a factor). To get an idea, the 2017 American Community Survey which is used to update some data within the primary Census.

1. Here are the available variables

acs2017_vars <- listCensusMetadata(name = "2017/pep/population/", 
                                 type = "variables")

head(acs2017_vars)

##      name                          label                concept predicateType
## 1    MDIV     Metropolitan Division code Selectable Geographies           int
## 2  COUSUB Minor Civil Division FIPS code              Estimates           int
## 3  SUMLEV                  Summary Level              Estimates          <NA>
## 4   STATE                State FIPS code              Estimates          <NA>
## 5 GEONAME                Geographic Name              Estimates          <NA>
## 6   DATE_                           Date              Estimates           int
##   group limit          required
## 1   N/A     0              <NA>
## 2   N/A     0              <NA>
## 3   N/A     0              <NA>
## 4   N/A     0              <NA>
## 5   N/A     0              <NA>
## 6   N/A     0 default displayed

2. Here are the available geographies

acs2017_geos <- listCensusMetadata(name = "2017/pep/population/", 
                                   type = "geography")

View(acs2017_geos)

If you want to see which variables are required for mapping, then there is a problem and that is actually due to bad luck as a result of timing. The Census is going through a new update with a new API format that was supposed to go live two weeks ago. It has, however, been pushed back to October 31. Until then, you will have to utilize the website. However if you are adamant about using R for the whole process as I typically am, then you can get table, aka tabular data from any website. If this interests you, please seethe section below. Otherwise move on to the next section.

Scraping the Web using R

In order to grab data off of the web, you will need to several things. At first it may seem arduous but having lengthy tables as data within R can be very beneficial, especially when filtering or if you just don’t want to manually load or type in entries. Included are steps to get the data you need in a usable form within R.

1. Download and run the rvest package (though there are others like a package called XML that may be useful in other circumstances)

library(rvest)

## Loading required package: xml2

## 
## Attaching package: 'rvest'

## The following object is masked from 'package:purrr':
## 
##     pluck

## The following object is masked from 'package:readr':
## 
##     guess_encoding

2. Define the website you want to pull the data from

webpage <- read_html("https://api.census.gov/data/2017/pep/population/variables.html")

3. Tell rvest that you want the tables using the command html_nodes. This does a lot of things but in a nutshell, it recognizes if a page is a webpage by scanning the HTML code and then finds what you need. In our case, we’re asking it to recognize the Census webpage from above and to identify the HTML codes where the term table exists.

tbls <- html_nodes(webpage, "table") 

4. See information regarding all of the tables

tbls # or if there are multiple tables, consider using head(tbls)

## {xml_nodeset (1)}
## [1] <table>\n<caption>Census Data API: Variables in /data/2017/pep/population ...

By the output, we only have one table which is indicated by [1].

5. Now do some piping via dplyr to get a list of everything within that table

tbls_ls <- webpage %>%
        html_nodes("table") %>% # Creates a multidimensional list (i.e. a list 
                                # with layers and depth akin to folders within 
                                # folders within folders, ect.)
        .[1] %>% # We only have table so we put.[1] here. If we had multiple 
                 # tables and we wanted say the first three, we could have 
                 # used .[1:3]
        html_table(fill = TRUE) # Convert the data to a table format and populate 
                                # the entries

6. See the result

View(tbls_ls)

Notice that its a list and can go pretty deep down the rabbit hole depending on how many tables you are looking at. In this case, we have one so ts pretty simple.

7. Now narrow down the results. Say you only wanted the ones that were required. You can figure that out by filtering for the term default displayed under the column Required. To do this, we first run the following which splits up all of the tables and then attaches them together by corresponding entries.

tbls_ls <- do.call(rbind.data.frame, tbls_ls) # Concatenates a list of data frames 
                                              # into a single data frame

8. (Optional) Take a look to see that it looks like an actual table

View(tbls_ls)

9. Now we get rid of the rows that that don’t have default displayed under Required

filtered_tbls_ls <- tbls_ls[tbls_ls$Required == "default displayed", ] 
                                                              # Filter out the rows that
                                                              # don't have `default 
                                                              # displayed` under the 
                                                              # column **Required**

10. Just show the names of the required variables

filtered_tbls_ls$Name # Just shows the Names of the required variables

## [1] "DATE_"    "UNIVERSE"

So here we see that the variables DATE and UNIVERSE are required.

Alternatively and if you understand the above, you can simply use the following in one fell swoop using dplyr:

library(rvest)

read_html("https://api.census.gov/data/2017/pep/population/variables.html") %>%
html_nodes("table")  %>% # Creates a multidimensional list (i.e. a list with layers 
                         # and depth akin to folders within folders within folders, 
                         # ect.)
  .[1] %>% # We only have table so we put.[1] here. If we had multiple tables and 
           # we wanted say the first three, we could have used .[1:3]
  html_table(fill = TRUE) %>% # Convert the data to a table format and populate the 
                              # entries
  bind_rows() %>% # Concatenates a list of data frames into a single data frame
  filter(Required == "default displayed") %>% # Filter out the rows that don't have
                                              # `default displayed` under the column 
                                              # **Required**
  select(Name) # Just shows the Names of the required variables

## New names:
## * NA -> ...9

##       Name
## 1    DATE_
## 2 UNIVERSE

Life is much easier with pipes! Note that this and the other method should work with many of the Census web tables but may not for all, especially for older tables that may have not been updated and do not have the column Required. Part of understanding what is going on in each line affords you the ability to change commands and variables as necessary.

OK now back to the mapping!

Getting Census Tract Level Data

Its not so obvious to get tract data sets. While there are multiple approaches, the easiest is to just ask for it. Here we’ll use the median household income from the 2015-2019 ACS in Monongalia County, WV.

mon <- get_acs(state = "WV", 
               county = "Monongalia", 
               geography = "tract", 
               variables = "B19013_001", 
               geometry = TRUE)

datatable(mon)

Here is the discrete scale with five colors…

wes_palette("Rushmore1")

…that has been converted to a continuous scale with 105,474 colors.

pal <- wes_palette("Rushmore1", 
                   max(mon$estimate), 
                   type = "continuous")

Single Map

We can plot the median household income using

mon %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf(color = NA) + 
  scale_fill_gradientn(colors = pal) +
  theme_map() 

Faceting Map

We can also break down a map using facet_wrap(). For example, let’s pull data for race, in particular for White, Black, Asian and Hispanic Populations.

racevars <- c(White = "P005003", 
              Black = "P005004", 
              Asian = "P005006", 
              Hispanic = "P004003")

mon_race <- get_decennial(geography = "tract",
                          variables = racevars,
                          state = "WV", 
                          county = "Monongalia County", 
                          geometry = TRUE,
                          summary_var = "P001001") 

Again we can use the datatable() command from the DT package to see an interactive view of the output

datatable(mon_race)

mon_race %>%
  mutate(pct = 100 * (value / summary_value)) %>%
  ggplot(aes(fill = pct)) +
  facet_wrap(~variable) +
  geom_sf(color = NA) + 
  scale_fill_gradientn(colors = pal) +
  theme_map() 

Well we know that West Virginia is primarily White resulting in a skewed distribution. To see how the other populations look proportionally by facet rather than overall, we can use a transformation in our scaling. In particular using a logarithmic transformation via trans="log10" when within a scale_ is sufficient. Remember that maps are generally considered to be approximations of the actual data, rather than an exact representation.

mon_race %>%
  mutate(pct = 100 * (value / summary_value)) %>%
  ggplot(aes(fill = pct)) +
  facet_wrap(~variable) +
  geom_sf(color = NA) + 
  scale_fill_gradientn(colors = pal,
                       trans="log10") +
  theme_map() 

Getting Census State and County Level Data

We need to first grab the data set

poverty_data <- getCensus(name = "acs/acs5/profile", # same as above
                          vintage = 2018,
                          vars = c("NAME", "DP03_0126PE"), # column/variable name
                          key = Sys.getenv("CENSUS_KEY"), # using the census key from before
                          region = "county:*",# any county
                          regionin = "state:*") # any state

Wrangling

A majority of your time is spent actually wrangling data. This often includes cleaning and restructuring it into a form you can us. For this example, we’ll do just that

First we’ll clean the data set

poverty_data_cleaned <- poverty_data %>%
                        separate(NAME, into = paste0('thing', 1:2), sep = '[,]') %>%
                        mutate(thing1 = str_replace_all(thing1, " County", "")) %>% # Remove the term County (with a space in front) 
                        mutate(thing2 = str_squish(thing2)) %>%
                        select(-c(state, county)) %>% # Get rid of these because we'll be using state and county names to merge
                        rename(county = thing1, 
                               state = thing2,
                               poverty_percentage = DP03_0126PE) %>%
                        mutate(poverty_percentage = replace(poverty_percentage,
                                                            which(poverty_percentage <= 0),
                                                            NA)) %>%
                        arrange(state, county)

Let’s take a look

poverty_data_cleaned %>% 
  head()

##    county   state poverty_percentage
## 1 Autauga Alabama               52.3
## 2 Baldwin Alabama               35.0
## 3 Barbour Alabama               60.6
## 4    Bibb Alabama               50.6
## 5  Blount Alabama               46.1
## 6 Bullock Alabama               66.5

So now that we have the figures, we need data about the mapping.

county_data <- map_data("county") %>% 
               rename(state = region, # rename because we'll use this to merge with the income_data_cleaned set
                      county = subregion) %>% # rename because we'll use this to merge with the income_data_cleaned set
               mutate(state = str_to_title(state)) %>% # make into a title case
               mutate(county = str_to_title(county)) # make into a title case 

Again, let’s take a look

county_data %>% 
  head()

##        long      lat group order   state  county
## 1 -86.50517 32.34920     1     1 Alabama Autauga
## 2 -86.53382 32.35493     1     2 Alabama Autauga
## 3 -86.54527 32.36639     1     3 Alabama Autauga
## 4 -86.55673 32.37785     1     4 Alabama Autauga
## 5 -86.57966 32.38357     1     5 Alabama Autauga
## 6 -86.59111 32.37785     1     6 Alabama Autauga

Now we can merge them together using joins

map_data <- left_join(county_data, 
                      poverty_data_cleaned, 
                      by = c("state","county")) 

map_data %>% 
  head()

##        long      lat group order   state  county poverty_percentage
## 1 -86.50517 32.34920     1     1 Alabama Autauga               52.3
## 2 -86.53382 32.35493     1     2 Alabama Autauga               52.3
## 3 -86.54527 32.36639     1     3 Alabama Autauga               52.3
## 4 -86.55673 32.37785     1     4 Alabama Autauga               52.3
## 5 -86.57966 32.38357     1     5 Alabama Autauga               52.3
## 6 -86.59111 32.37785     1     6 Alabama Autauga               52.3

Mapping Using Estimated Curvatures

Let’s go ahead and map the data

estimate_map <- 
  ggplot() + 
  geom_polygon(data = map_data, aes(x = long, 
                                    y = lat, 
                                    fill = poverty_percentage, 
                                    group = group),
               color = "#f8f8fa", 
               size = 0.20, 
               show.legend = TRUE) +
  geom_polygon(data=subset(map_data("state")), 
               aes(x = long, 
                   y = lat, 
                   group = group), 
               size = 0.5,
               color="white", 
               fill = NA) +
  scale_fill_gradient(name="Estimates", 
                      low = "#fa9fb5", 
                      high = "#7a0177",
                      breaks = c(0, 25, 50, 75, 100),
                      limits = c(0,100), 
                      na.value = "#c9d3d9",
                      labels = c(str_wrap("0.00% - 24.99%", 10), 
                                 str_wrap("25.00 - 49.99%", 10),
                                 str_wrap("50.00% - 74.99%", 10),
                                 str_wrap("75.00% - 99.99%", 10),
                                 "100%"),
                      guide = "legend") +
  theme(plot.title = element_text(size = 13,
                                  color ="#4682b4", 
                                  vjust = -1,
                                  hjust = 0.5,
                                  face = "bold"),
        plot.caption = element_text(size = 8,
                                    color ="#4682b4",
                                    face = "italic",
                                    margin = margin(t = 1, 
                                                    unit = "cm")),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.key = element_rect(size = 5, 
                                  color = NA),
        legend.key.width = unit(2, 
                               "cm"),
        legend.key.height = unit(0.3, 
                               "cm"),
        legend.text = element_text(size = 8, 
                                   color = "#4682b4"),
        legend.title.align = 0.5,
        legend.title = element_text(size = 10, 
                                    color = "#4682b4"),
        legend.background = element_rect(linetype = 0, 
                                         size = 0.5, 
                                         color = 1),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  labs(title = str_to_title(str_wrap("2014-2018 estimate of families with female householder, no husband present with related children of the householder under 18 years", 50)),
       caption = "Source: 2018 American Community Survey: 5-Year Estimates") +
  guides(fill = guide_legend(title.position = "top",
                             label.position = "bottom", 
                             label.hjust = 0.5)) 

estimate_map

Well that’s the default map projection and its just terrible. You can fool around with the representation, or projection of the map using any number of coordinate systems via coord_map(). Below are just a few but a full list can be found on the mapproj CRAN site.

Albers projection

estimate_map + 
  coord_map("albers", 
            lat0 = 39, 
            lat1 = 45)

projection with an approximated aspect ratio

estimate_map + 
  coord_quickmap()

Lambert azimuthal equal-area projection

estimate_map + 
  coord_map("azequalarea")

Bonne projection

estimate_map + 
  coord_map("bonne", 
            lat0 = 50)

cylindrical projection

estimate_map + 
  coord_map("cylindrical")

Gilbert projection

estimate_map + 
  coord_map("gilbert")

LaGrange projection

estimate_map + 
  coord_map("lagrange")

estimate_map + 
  coord_map("mercator")

estimate_map + 
  coord_map()

orthographic projection

estimate_map + 
  coord_map("orthographic")

stereographic projection

estimate_map + 
  coord_map("stereographic")

Saving

You can even save any of the plots as a pdf on standard A4 paper in landscape mode

poverty_estimates_albers <- estimate_map + 
                            coord_map("mercator")

ggsave("poverty_estimates.pdf",
       poverty_estimates_albers,
       width = 297, 
       height = 210, 
       units = "mm")