038_facets_wrap_thematic_mapping.Rmd

---
title: "Faceted Thematic Mapping"
output: html_document
---


```{r load-libraries, message=FALSE, warning=FALSE}
library(tidyverse)       # Tidyverse for Tidy Data
library(readxl)
library(tmap)            # Thematic Mapping
library(tmaptools) 
library(tigris)          # Get Census Geography Poloygons
library(sf)
library(tidycensus)
library(viridis)        # Loading required package: viridisLite
library(ggplot2)        #-- version 2.2.1.90000
              #devtools::install_github("tidyverse/ggplot2")
```


## Shapefiles as sf 

Using the `tigris` package get Census Tiger shapefiles for census geographies.  Tigris will return the shapefile in the `sf`, or simple features, format.

```{r get-data, message=FALSE, warning=FALSE}
us_geo <- tigris::states(class = "sf")
# us_geo_spdf <- states()  # spdf uses @data slots -- old school
```

## Get BLS data

I've already downloaded and stored some data from the Bureau of Labor Statistics.  Thoses data are stored in an excel file in the `data` directory of the repository:  `data/OES_Report.xlsx`.  **The goal is to attach this data to the previously downloaded shapefiles.**  

But you may be intersted in how I gathered the data. below are some summary notes documenting my steps of gathering the data from the Bureau of Labor Statistics.  

https://data.bls.gov/oes/#/occGeo/One%20occupation%20for%20multiple%20geographical%20areas

-  One occupation for multiple geographical areas

    - Mental Health and Substance Abuse Social Workers  (Also, Secondary School Teacher, Waiter, Legislator, and Paralegals)
    
        - State
        
        - All States in this list
        
        - Annual Mean wage
        
            - Excel
            
- Read the Data in with the RStudio "Import Dataset" wizard available in the *Environment* tab.  This will generate the code below and ensure the import

    - Skips the first 4 lines
    - Coerces the 2nd column to numeric

```{r get-blsdata, message=FALSE, warning=FALSE}
Salary4Helpers <- 
  read_excel("data/OES_Report.xlsx",
             col_types = c("text", "numeric"), 
             skip = 4)

# Salary4Helpers
```

### Load more BLS Wages 

read in the XLXS files for the other occupations:  Secondary School Teacher, Waiter, Legislator, and Paralegals

```{r other_occupations_data, message=FALSE, warning=FALSE}
legislators <- 
  read_excel("data/OES_legislator.xlsx",
             col_types = c("text", "numeric"), 
             skip = 4)

paralegals <- 
  read_excel("data/OES_paralegals.xlsx",
             col_types = c("text", "numeric"), 
             skip = 4)

teachers <- 
  read_excel("data/OES_secondary_school_teacher.xlsx",
             col_types = c("text", "numeric"), 
             skip = 4)

waiters <- 
  read_excel("data/OES_waiters.xlsx",
             col_types = c("text", "numeric"), 
             skip = 4)

```

        

### Wrangle the data

Before we join the BLS data to the shapefile we need to transform the structure of the downloaded BLS data

```{r wrangle_bls_data, message=FALSE, warning=FALSE}
BlsWage_ToJoin <- Salary4Helpers %>% 
  rename(Area = "Area Name") %>% 
  rename(wages = "Annual mean wage(2)") %>% 
  mutate(State = gsub("\\(\\d{7}\\)", "", Area)) %>% 
  filter(wages != "NA_character_") %>% 
  select(State, wages)

#BlsWage_ToJoin
```


```{r wrangle_other_ocupations}


legislators <- legislators %>% 
  rename(Area = "Area Name") %>% 
  rename(wages = "Annual mean wage(2)") %>% 
  mutate(State = gsub("\\(\\d{7}\\)", "", Area)) %>% 
  filter(wages != "NA_character_") %>% 
  select(State, wages)

paralegals <- paralegals %>% 
  rename(Area = "Area Name") %>% 
  rename(wages = "Annual mean wage(2)") %>% 
  mutate(State = gsub("\\(\\d{7}\\)", "", Area)) %>% 
  filter(wages != "NA_character_") %>% 
  select(State, wages)

teachers <- teachers %>% 
  rename(Area = "Area Name") %>% 
  rename(wages = "Annual mean wage(2)") %>% 
  mutate(State = gsub("\\(\\d{7}\\)", "", Area)) %>% 
  filter(wages != "NA_character_") %>% 
  select(State, wages)

waiters <- waiters %>% 
  rename(Area = "Area Name") %>% 
  rename(wages = "Annual mean wage(2)") %>% 
  mutate(State = gsub("\\(\\d{7}\\)", "", Area)) %>% 
  filter(wages != "NA_character_") %>% 
  select(State, wages)

```



### Append data

Using the `append_data()` function of the `tmaptools` package, append BLS data to the previously loaded shape object

```{r make_nice_df}
paralegals <- append_data(us_geo, paralegals, 
                                 key.shp = "NAME",
                                 key.data = "State")

waiters <- append_data(us_geo, waiters, 
                                 key.shp = "NAME",
                                 key.data = "State")

teachers <- append_data(us_geo, teachers, 
                                 key.shp = "NAME",
                                 key.data = "State")

legislators <- append_data(us_geo, legislators, 
                                 key.shp = "NAME",
                                 key.data = "State")
```

```{r append_bls_data}
HelperShapeObject <- append_data(us_geo, BlsWage_ToJoin, 
                                 key.shp = "NAME",
                                 key.data = "State")
```



### Contiguous 48 states

Filter to only the contiguous 48 states + D.C.

```{r subset-data}
contiguous_states <- HelperShapeObject %>% 
  filter(REGION != 9) %>% 
  filter(STUSPS != "AK") %>% 
  filter(STUSPS != "HI")


paralegals <- paralegals %>% 
  filter(REGION != 9) %>% 
  filter(STUSPS != "AK") %>% 
  filter(STUSPS != "HI")

waiters <- waiters %>% 
  filter(REGION != 9) %>% 
  filter(STUSPS != "AK") %>% 
  filter(STUSPS != "HI")

teachers <- teachers %>% 
  filter(REGION != 9) %>% 
  filter(STUSPS != "AK") %>% 
  filter(STUSPS != "HI")

legislators <- legislators %>% 
  filter(REGION != 9) %>% 
  filter(STUSPS != "AK") %>% 
  filter(STUSPS != "HI")
```


## Get Population Data -- Tidycensus

### identify and pick census variables

B01003_001E = Total Population
B06011_001E = Median income in the past 12 months --!!Total:

```{r identify_ACS_variables}
variables_census <- load_variables(2015, "acs5", cache = TRUE)
```

### Get Median Income

```{r get_median_income}
us_pop <- 
  get_acs(geography = "state",
          variables = "B06011_001E",
          # state = "NC",
          geometry = FALSE)

us_pop <- us_pop %>% 
  filter(NAME != "Alaska") %>% 
  filter(NAME != "Hawaii") %>% 
  filter(NAME != "Puerto Rico")
```

### Append Census Data

```{r make_nice_df2}
appended_population_and_states_shapefile <- append_data(contiguous_states, us_pop, 
                                 key.shp = "GEOID",
                                 key.data = "GEOID")
```


## Make Tall

Wrangle the data into a tall format and single data frame

```{r make_tall}
contiguous_states_slim <- contiguous_states %>% 
  select(GEOID, NAME, wages, geometry) %>% 
  mutate(category = "Substance Abuse Soc.Wrk") %>% 
  rename(variable = wages)
```

```{r wrangle_pop}
population_states_slim <- appended_population_and_states_shapefile %>% 
  select(GEOID, NAME, estimate, geometry) %>% 
  mutate(category = "Median Income") %>% 
  rename(variable = estimate)
```


```{r wrangle_other_bls}
paralegals <- paralegals %>% 
  select(GEOID, NAME, wages, geometry) %>% 
  mutate(category = "Paralegals") %>% 
  rename(variable = wages)

waiters <- waiters %>% 
  select(GEOID, NAME, wages, geometry) %>% 
  mutate(category = "Waiters") %>% 
  rename(variable = wages)

teachers <- teachers %>% 
  select(GEOID, NAME, wages, geometry) %>% 
  mutate(category = "Teachers") %>% 
  rename(variable = wages)

legislators <- legislators %>% 
  select(GEOID, NAME, wages, geometry) %>% 
  mutate(category = "Legislators") %>% 
  rename(variable = wages)
```

### Bind Rows 

To combine all the smaller data frames into one tall format, bind the rows.

base::`rbind` works.  `dplyr::bind_rows` did not work.

```{r bind_tall}
tall <- rbind(contiguous_states_slim,
              population_states_slim,
              paralegals,
              waiters,
              teachers,
              legislators)
```

### More Munging

Make the category variable a categorical factor with levels.  This will improve the order of the facets when displayed.

```{r}
display_levels <- c("Median Income", "Legislators", 
                    "Paralegals", "Substance Abuse Soc.Wrk", 
                    "Teachers", "Waiters")

tall <- tall %>% 
  mutate(category = factor(category, display_levels)) %>% 
  rename(Wages = variable)
```



## Facet_grid

`facet_wrap` works better than `facet_grid` because there is more control over the display of the rows and columns.  The first example, below, does not run.  I keep it only to document something that could work.

```{}
tall %>% 
  ggplot(aes(fill = variable, color = variable)) +
  geom_sf() +
  coord_sf(crs = 5070) +
  scale_fill_viridis(option = "viridis") + 
  scale_color_viridis(option = "viridis") +
  facet_grid(. ~ category)
```

### Display Map

```{r facet_wrap_it_up}
tall %>% 
  ggplot(aes(fill = Wages, color = Wages)) +
  geom_sf() +
  coord_sf(crs = 5070, datum = NA) +
  scale_fill_viridis(option = "viridis") + 
  scale_color_viridis(option = "viridis") +
  facet_wrap(~ category, nrow = 3, ncol = 2)
```

### Save Map

```{r save_it}
ggsave("facet_map.png", width = 8, height = 8, units = "in")
```