```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(pacman)
p_load(tidyverse, tigris, sp, sf, janitor, scales, broom, cartogram, tidycensus)
p_load_current_gh("hrbrmstr/statebins")
p_install(usmap, force = FALSE)
p_install(here, force = FALSE)

data_dir = here::here("data")

theme_set(
    theme_classic() +
    theme(axis.line = element_blank(), axis.ticks = element_blank(),
          axis.text = element_blank(), axis.title = element_blank(),
          legend.position = c(1.0,0.5), legend.justification = c(1,1))
  )
```

# Introduction

# Wrangling the data

We take county-level vote data for the 2016 presidential election from the
[MIT Election Data and Sceince Lab](https://electionlab.mit.edu/) (MEDSL),
which is conveniently available as a `.csv` file at
GitHub: <https://github.com/MEDSL/county-returns>

Load the full election data set from the `.csv` file, and filter to keep
only the rows for 2016 (the full set includes presidential elections from
2008--2016).

Throughout this project, we will use FIPS (Federal Information Processing 
Standards) codes to identify states, counties, and census tracts. FIPS
codes use two-digit codes to identify each state and territory in the US.
Within each state each county receives a three-digit code and a county is 
uniquely identified by a five-digit FIPS code, where the first two digits 
identify the state and the next three identify the county.

We have to fix up the FIPS (Federal ) code in the presidential election data 
because the `.csv` file from MEDSL represents them as integers, 
so state FIPS codes less than 10 are represented by a single-digit and need to 
be converted to character data with a leading zero so they will have two digits.
For this task, we use the formatC command.

Also, the party affiliation of third-party and independent candidates are 
represented by `NA`, and we need to convert this to a character value, which
we will set equal to `"other"` using the `tidyverse` function `replace_na()`.

```{r load_data, cache = TRUE, cache.extra = c(file.mtime("data/MEDSL/county-returns/countypres_2000-2016.csv"))}
election_2016 <-
  read_csv("data/MEDSL/county-returns/countypres_2000-2016.csv") %>%
  filter(year == 2016) %>%
    mutate(party = replace_na(party, "other"),
           FIPS = formatC(FIPS, digits = 4, format = "d", flag = "0")) %>%
  filter(! is.na(state_po))
```

Let's take a look at the first few rows of the data frame:

```{r peek_raw_data, cache=TRUE, dependson = "load_data"}
glimpse(election_2016)
```

The original data frame is in a "long" format with one row for each party in
each county, but we want a column with the two-party vote margin, so we
to calculate the margin need to pivot the data frame to a "wide" format with
one row per state and one column for each party:

```{r wrangle_2016_data}
election_2016 <- election_2016 %>%
  select(state = state_po, party, county, FIPS, candidatevotes) %>%
  pivot_wider(names_from = "party", values_from = "candidatevotes") %>%
  filter(! is.na(democrat))

head(election_2016)
```

Next, we add the counties within a state to create a state-level data set.

For both state-level and county-level data, we create three new columns:
`total`, for the total votes cast, `margin` for the margin of victory as a
fraction of the total two-party vote (positive for Democratic victory and 
negative for Republican victory), and `winner` indicating the winner in 
that state or county.

Finally, a nagging detail we must take care of 
(data cleaning is very important!):
Shannon County SD (FIPS code 46113) changed its name to Oglala Lakota County
in 2015 and received a new FIPS code (46102), but the MEDSL data files are out 
of date and still list it as Shannon County with FIPS code 46113 so we need to 
fix this to make it consistent with the map and population data.

```{r state_data, cache=TRUE, dependson="load_data"}
state_votes <- election_2016 %>% group_by(state) %>%
  summarize_at(vars(democrat, republican, other), ~sum(., na.rm = TRUE)) %>%
  ungroup() %>% # 1 row for each group
  mutate(total = democrat + republican + other,
         margin = (democrat - republican) / (democrat + republican),
         winner = ifelse(democrat > republican, "Clinton", "Trump"))

county_votes <- election_2016 %>%
  mutate(total = democrat + republican + other,
         margin = (democrat - republican) / (democrat + republican),
         winner = ifelse(democrat > republican, "Clinton", "Trump"),
         FIPS = ifelse(FIPS == "46113", "46102", FIPS))

head(state_votes)
```

Now we have the data and we want to get polygons with the shapes of the states
and counties. We use the U.S. Census's TIGR maps, which we get with the
`tigris` package and then we convert them to `sf` ("**simple feature**") data 
frames.

**Simple feature** data frames look very much like the `tibble` data frames 
you are used to from _R for Data Science_, but they have one special column 
called `geometry`, which is treated differently from the other columns. 
This column contains the shapes for the _simple features_. Each row can have
a polygon, a set of multiple polygons, a point, a set of points, 
or other types of spatial data. Each row can have a different type 
of _simple feature_. The term _simple feature_ refers to a formal standard
(ISO 19125-1:2004), that specifies a uniform digital representations 
of spatial geometry objects. This specification serves as the basis for 
geometrical representations in popular geospatial software inclusing ArcGIS,
the PostGIS database, and the GeoJSON file format.

```{r load_maps, cache=TRUE, dependson="state_data", message = FALSE}
options(tigris_use_cache = TRUE, tigris_refresh = FALSE)

state_map_file  <- file.path(data_dir, "maps", "state_maps.Rds")
county_map_file <- file.path(data_dir, "maps", "county_maps.Rds")

if (file.exists(state_map_file)) {
  state_maps <- read_rds(state_map_file)
} else {
  state_maps <- states(cb = TRUE, year = 2016, progress_bar = FALSE, 
                       class = "sf") %>%
    clean_names() %>% rename(state = stusps)
  write_rds(state_maps, state_map_file)
}

if (file.exists(county_map_file)) {
  county_maps <- read_rds(county_map_file)
} else {
  county_maps <- counties(state = c(state.abb, "DC"), cb = TRUE, year = 2016, 
                          progress_bar = FALSE, class = "sf") %>%
    clean_names()
  write_rds(county_maps, county_map_file)
}

glimpse(state_maps)
```

Now we join the map data frames to the election return data frames.
This will produce data frames that combine the voting return data with 
geometric representation of states or counties for mapping.

```{r join_data, cache=TRUE, dependson=c("state_data", "load_maps")}
state_maps <- state_maps %>% inner_join(state_votes, by = "state")
county_maps <- county_maps %>% inner_join(county_votes, by = c(geoid = "FIPS"))
```

# Maps

Now we are ready to make some maps.

For graphical simplicity, we will leave off Alaska and Hawaii.
If you want to make compact maps that display these states offset to fit
nicely around the contiguous 48 states, look at the `usmaps` package.

We also project the maps onto a plane using the NAD83 projection. 
This uses the 1983 North American Datum, an Albers equal-area conic
projection optimized to minimize distortions in the 48 contiguous states.

```{r omit_ak_hi, cache = TRUE}
state_maps <- state_maps %>% filter(! state %in% c("AK", "HI")) %>%
  st_transform("+init=epsg:5070")   # NAD83

county_maps <- county_maps %>% 
  filter(! statefp %in% c("02", "15"), ! state %in% c("AK", "HI")) %>%
  st_transform("+init=epsg:5070")  # NAD83
```

## State-Level Maps

First, we can make a simple map that shows the winner of each state as red or 
blue:

```{r state_winner, cache = TRUE, dependson="omit_ak_hi"}
  ggplot(state_maps, aes(fill = winner)) + geom_sf() +
  scale_fill_manual(values = c(Clinton = "blue", Trump = "red"), 
                    name = "Winner")
```

This gives a strong sense of how many states went red or blue in the election,
and since each state except Nebraska and Maine award all of their electoral 
votes to the winner of a plurality of votes, it is informative.

#### Mapping the Margin of Victory

But we may want to know something about the margin of victory in each state.
Here we map the margin of victory (Democratic minus Republican votes divided
by the two-party total) to a gradient that goes from bright red to bright blue
with white for a margin of zero (a tie):

```{r state_margin, cache=TRUE, dependson="omit_ak_hi"}
  ggplot(state_maps, aes(fill = margin)) + geom_sf() +
  scale_fill_gradient2(high = "blue", low = "red", mid = "white",
                       name = "Vote margin", labels = percent_format())
```

This gives a clearer sense that some states had a much greater margin than 
others, but many close states just appear in white, so we don't know who won 
them.

What may be helpful here would be to clip the vote margin at 10% or so (which
is a landslide margin already), in order to have better color discrimination
in the close states:

```{r state_margin_clipped, cache=TRUE, dependson="omit_ak_hi"}
  state_maps %>% 
  mutate(margin = raster::clamp(margin, -0.10, 0.10)) %>%
  ggplot(aes(fill = margin)) + geom_sf() +
  scale_fill_gradient2(high = "blue", low = "red", mid = "white",
                       name = "Vote margin",
                       breaks = seq(-0.1, 0.1, 0.05),
                       labels = c("< -10%", "-5%", "0", "5%", "> 10%"))
```

This gives a clearer sense of the election with most states going to one
candidate or the other by a very large margin, but some states being much 
closer.
A few states are still too close to tell with this color scheme: Michigan,
New Hampshire, Pennsylvania, and Wisconsin all had margins of victory less 
than 1%:

```{r mi_pa_nh_margin, cache=TRUE, dependson="omit_ak_hi"}
state_votes %>% filter(state %in% c("MI", "NH", "PA", "WI")) %>%
  select(state, democrat, republican, margin) %>%
  arrange(abs(margin)) %>%
  mutate(margin = percent(margin, accuracy = 0.01))
```

An alternate visualization would use a muted purple instead of white as the
intermediate color for states with small margins of victory for either party:

```{r state_margin_purple, cache=TRUE, dependson="omit_ak_hi"}
  state_maps %>% mutate(margin = raster::clamp(margin, -0.10, 0.10)) %>%
  ggplot(aes(fill = margin)) + geom_sf(color = "gray90") +
  scale_fill_gradient2(high = "blue", low = "red", mid = muted("purple"),
                       name = "Vote margin",
                       breaks = seq(-0.1, 0.1, 0.05),
                       labels = c("< -10%", "-5%", "0", "5%", "> 10%"))
```

I prefer white as the intermediate color because I have a very hard time
telling whether a shade or purple for, say, Minnesota or Colorado, is 
redder or bluer than the neutral value, but many others prefer the 
purple palette because it has a more uniform lightness/darkness across the 
different colors.

### Color Blindness and Color Palettes

This is also a good time to mention that people with _protanopia_ 
color-blindness can have difficulty distinguishing certain shades of red from 
shades of blue, so despite the popularity of drawing maps with red and blue to
indicate Republican and Democratic victories, respectively, it may be useful to
draw maps using a color palette designed to be useful for people with common 
types of color-blindhess. The `RColorBrewer` and `viridisLiote` packages 
include such palettes.

Here, the `RColorBrewer` palette `RdBu` gives a gradient from red to white to
blue, but with the specific hues chosen to be readable by someone with 
_protanopia_:

```{r state_margin_RdBu, cache=TRUE, dependson="omit_ak_hi"}
  state_maps %>% mutate(margin = raster::clamp(margin, -0.10, 0.10)) %>%
  ggplot(aes(fill = margin)) + geom_sf(color = "gray90") +
  scale_fill_distiller(palette = "RdBu", direction = 1,
                       name = "Vote margin",
                       breaks = seq(-0.1, 0.1, 0.05),
                       labels = c("< -10%", "-5%", "0", "5%", "> 10%"))
```

### Alternate Representation of States

One problem with state maps is that states with large areas have much greater
visual salience than states with smaller area, even if the smaller states have
greater population and more electoral votes.
Rhode Island is almost invisible on this map even though it has one more 
electoral vote than Montana, North Dakota, South Dakota, and Wyoming, all of 
which are very large and visible.

One approach to addressing this problem are either to represent each state
by a square of equal size, located approximately in the same relative location 
to other states as its location on a real map. Here I use the `statebins`
package by [Bob Rudis](https://rud.is/b/) 
(`hrbrmstr` on [GitHub](https://github.com/hrbrmstr), 
[GitLab](https://gitlab.com/hrbrmstr), 
[Twitter](https://twitter.com/hrbrmstr), etc.)

```{r state_bins, cache=TRUE, dependson="omit_ak_hi"}
  state_votes %>% filter(! state %in% c("AK", "HI")) %>%
  ggplot(aes(state = state, fill = winner)) + geom_statebins() +
  # scale_fill_manual(values = c(Clinton = "blue", Trump = "red"), 
  # name = "Winner") +
  scale_fill_manual(name = "Winner",
                    values = set_names(brewer_pal(palette = "RdBu")(3), 
                                       c("Trump", "", "Clinton"))) +
  coord_equal() + theme_statebins()
```

```{r state_bins_margin, cache=TRUE, dependson="omit_ak_hi"}
  state_votes %>% filter(! state %in% c("AK", "HI")) %>%
  mutate(margin = raster::clamp(margin, -0.10, 0.10)) %>%
  ggplot(aes(state = state, fill = margin)) + geom_statebins() +
  scale_fill_distiller(palette="RdBu",
                       breaks = c(-0.10, -0.05, 0.0, 0.05, 0.10),
                       labels = c("< -10%", "-5%", "0", "5%", "> 10%")) +
  coord_equal() + theme_statebins() + theme(legend.position = "right")
```

### Cartograms

Another approach is to use what are called **cartograms**, which distort the 
map to make areas of the states proportional to some quantity, such as
population or electoral votes. Here, we create a cartogram of the states
weighted so their areas are proportional to the population.

We take the population from the `usmaps` package, which contains the 2015
population of each state, based on the U.S. Census Bureau's 2015 American 
Community Survey. 

If we wanted to be more precise, we could use the `tidycensus` package to 
download the 2016 state populations from the Census Bureau, but 2015 will
be close enough for this demonstration project.

```{r state_cartogram, message = FALSE, cache=TRUE, dependson="omit_ak_hi"}
state_pop <- usmap::statepop %>% 
  select(statefp = fips, pop = pop_2015)

state_maps %>%
  left_join(state_pop, by = "statefp") %>%
  select(winner, margin, pop) -> new_state_data

# Cartogram weighted by number of votes in each state
state_cartogram_file <- file.path(data_dir, "maps", "state_cartogram.Rds")
if (file.exists(state_cartogram_file)) {
  cgram_st_c <- read_rds(state_cartogram_file)
} else {
  cgram_st_c <- cartogram_cont(new_state_data, "pop")
  write_rds(cgram_st_c, state_cartogram_file)
}

ggplot(cgram_st_c, aes(fill = margin)) + geom_sf() +
  scale_fill_gradient2(low = "red", high = "blue", mid = "white",
                       breaks = seq(-0.10, 0.10, 0.05),
                       labels = c("< -10%", "-5%", "0", "5%", "> 10%"))
```

## Making Effective Maps

There is no right answer to the question of what quantities to use in coloring
maps or what color palette to use. A lot depends on the purpose for which you 
are making the map, who your target audience is, and your personal aesthetic 
judgment. 

It is useful to read lots of maps so you can develop an informed sense of the 
range of choices others have used. The more you know about the possibilities,
the more you can make deliberate and informed choices in your own maps rather 
than just relying on the defaults of some piece of software.

Some of my recommendations are:

* Cynthia A. Brewer, _Designing Better Maps: A Guide for GIS Users_, 
  2nd edition (Esri Press, 2016)
* Gretchen N. Peterson, _GIS Cartography: A Guide to Effective Map Design_,
  2nd edition (CRC Press, 2015)
* Mark Monmonier: many books, including 
  * _Mapping It Out: Expository Cartography for the Humanities and Social 
    Sciences_ (U. Chicago Press, 1993),
  * _How to Lie with Maps_, 3rd edition (U. Chicago Press, 2018),
  * _Cartographies of Danger: Mapping Hazards in America_ (U. Chicago Press, 1997),
  * _Air Apparent: How Meteorologists Learned to Map, Predict, and Dramatize 
    Weather_ (U. Chicago Press, U. Chicago Press 1999),
  * _Coastlines: How Mapmakers Frame the World and Chart Environmental Change_
    (U. Chicago Press, 2008),
  * and _Bushmanders and Bullwinkles: How Politicians Manipulate Electronic 
    Maps and Data to Win Elections_ (U. Chicago Press, 2001).


## County-Level Mapping

All the states except Nebraska and Maine award all of their electoral votes to
the candidate who receives the most total votes. Thus, state-level maps are
the most directly relevant to understanding the winner, but states are not
homogeneous and it can be useful to examine maps at the county level to see
how votes are distributed within the states.

First, we can make a simple map that shows the winner of each state as red or 
blue:

```{r county_winner, cache = TRUE, dependson="omit_ak_hi"}
county_maps %>% 
  mutate(winner = factor(ifelse(margin > 0, "Clinton", "Trump"))) %>%
  ggplot(aes(fill = winner)) + geom_sf(color = "gray50") +
  scale_fill_manual(values = c(Clinton = "blue", Trump = "red"), 
                    name = "Winner")
```

This is a map strongly favored by President Trump because it shows that he
won the overwhelming majority of counties in the US.

!["Try to impeach this!" Trump tweet 1 Oct. 2019](assets/images/trump_tweet_impeach_this.png)

Critics quickly pointed out first that President Trump had tweeted a county 
election map from 1968 instead of 2016 and also that this way of presenting 
a map did not account for the numbers of votes, either with regard to the 
margin of victory or the different populations in different counties.

Thus, it may be useful to compare this map with one that uses a gradient color
scale to represent the margin of victory in each county, much as we did for the 
state-level maps:

```{r county_margin, cache=TRUE, dependson="omit_ak_hi"}
county_maps %>%
  ggplot(aes(fill = margin)) + geom_sf() +
  scale_fill_gradient2(high = "blue", low = "red", mid = "white",
                       name = "Vote margin", labels = percent_format())
```

And here is the same map, but using a muted purple instead of white for a margin
of 0%:

```{r county_margin_purple, cache=TRUE, dependson="omit_ak_hi"}
  county_maps %>%
  ggplot(aes(fill = margin)) + geom_sf() +
  scale_fill_gradient2(high = "blue", low = "red", mid = muted("purple"),
                       name = "Vote margin", labels = percent_format())
```

### County-Level Cartograms

And here is a cartogram, with the area of each county weighted by its 
population. This cartogram has some complications. 

First, because of the large number of counties, the algorithm becomes very 
computationally demanding and calculating the cartogram takes several hours.
Second, the algorithm does not work very effectively with an `sf` simple-feature
data frame when there is this much data. The resulting cartogram was obviously
wrong. I eventually discovered that I had to convert the data frame to a 
different R class, called `Spatial`, using the `sf` function `as_Spatial`. 
Once I did this, the `cartogram` function could process it properly, but it 
still took several hours to calculate the cartogram.

```{r county_cartogram, message = FALSE, cache=TRUE, dependson="omit_ak_hi"}
county_pop <- usmap::countypop %>% select(geoid = fips, pop = pop_2015)

county_maps %>% 
  mutate(winner = ifelse(margin > 0, "Clinton", "Trump"),
         vpop = democrat + republican + other) %>%
  inner_join(county_pop, by = "geoid") %>%
  select(winner, margin, pop) -> new_data

# Cartogram weighted by number of votes in each state
county_cartogram_file <- file.path(data_dir, "maps", "county_cartogram.Rds")
if (file.exists(county_cartogram_file)) {
  cgram_cty_csp <- read_rds(county_cartogram_file)
} else {
  cgram_cty_csp <- as_Spatial(new_data) %>% 
    cartogram_cont("pop", prepare = "remove")
  write_rds(cgram_cty_csp, county_cartogram_file)
}
```
```{r county_cartogram_winner, message=FALSE, dependson="county_cartogram"}
cgram_cty_csp %>% as("sf") %>% ggplot(aes(fill = winner)) + 
  geom_sf(color = NA) +
  scale_fill_manual(values = c(Clinton = "blue", Trump = "red"),
                    name = "Winner")
```

And the same cartogram, but colored by the margin of victory:

```{r county_cartogram_margin, message=FALSE, dependson="county_cartogram"}
cgram_cty_csp %>% as("sf") %>% ggplot(aes(fill = margin)) + 
  geom_sf(color = NA) +
  scale_fill_gradient2(low = "red", high = "blue", mid =  "white")
```

And with muted purple as the intermediate color:

```{r county_cartogram_purple, dependson="county_cartogram"}
  cgram_cty_csp %>% as("sf") %>% ggplot(aes(fill = margin)) + 
  geom_sf(color = NA) +
  scale_fill_gradient2(low = "red", high = "blue", mid =  muted("purple"))
```

The cartograms accurately represent the number of votes in each county,
but critics often point out that this is not very relevant to presidential 
elections because the Constitution does not award electoral votes in proportion
to the populaion: California has 67 times the population of Wyoming but it 
only has 18 times as many electoral votes (55 versus 3). 

Thus, maps that emphasize the number of votes in each state or county does not
accurately reflect the winner or loser in the electoral college.

## Point Density Maps

Cartograms represent the effect of population, but distort maps to the point
that features can become unrecognizable. Looking back to the state-level 
cartogram, how easy is it to recognize New Mexico, New Hampshire, 
Wyoming, or Idaho?

Another way to represent the uneven distribution of population is to use a 
point-density map in which the states and counties retain their actual
sizes and shapes but the voters are represented as points.

This is closely related to a concept called _dasymetry_ in which population data
that is measured as aggregates over certain areas is distributed as points, 
using probabilities taken from ancillary data.

Professor Katherine Nelson, now in the Department of Geography and Geospatial
Sciences at Kansas State University, developed sophisticated techniques during 
her Ph.D. research at Vanderbilt to map population data in Nashville from the 
U.S. Census to specific dwellings by cross-referencing it against maps of 
residential housing units obtained from the Nashville-Davidson County Property 
Assessor. 
This allowed Dr. Nelson to then assess vulnerability of different demographic 
groups to flooding based on how specific houses and apartment buildings were
situated relative to the flood plain and on statistical models of how likely 
different demographic groups (by age, race, socio-economic status, etc.)
were to occupy those dwellings.

Another simmilar project, led by Professor Marshall Shepherd at the 
University of Georgia and his graduate student Neil Debbage (now an 
assistant professor in the Department of Political Science and Geography
at the University of Texas at San Antonio), 
which applied dasymetric mapping to study racial disparaties in vulnerability 
to flooding throughout the _Charlanta Megaregion_ that stretches from 
Charlotte NC to Atlanta GA. 

To keep things simple, here I will compare non-dasymetric to dasymetric 
techniques to represent votes as point-density maps.

More than 136 million people voted in the 2016 presidential election, and it 
will be very difficult, time-consuming, and computationally demanding to 
draw a figure with 136 million points. Thus, I will sample a random fraction of 
the votes to make these figures.

First, I have actual vote counts for each county. The median number of 
votes in the different counties ranges from 65 to more than 3 million, with 
a median of about 11,000. For my first point-density map, I will draw one 
red point for each 1000 Republican votes in a county, 
one blue point for every 1000 Democratic votes, 
and one green point for every 1000 votes for third-party or independent 
candidates. I will place each point randomly within the boundaries of the 
county. The `st_sample()` function makes this easy.

As I was writing the code to draw the map, I discovered an interesting 
problem: Even counties that I knew had voted overwhelmingly for Hillary Clinton
appeared red. Then I realized that `ggplot()` plots points in the order they 
appear in the data frame and I had built my data frame by combining three
data frames with points corresponding to Democratic, Republican, and 
Other voters, respectively, which meant that the Democratic (blue) points were
plotted first and the Republican (red) points were plotted on top of them.

Where there was a high density of points, the red points on top would almost 
completely obscure the blue points underneath, even if there were many more 
blue points.

To solve this, I added a new step to the data wrangling: After I combined the 
three data frames, I shuffled the rows using the `tidyverse` function 
`sample_frac()` to randomly sample all of the rows without replacement
(i.e., shuffle the rows into a random order). Now the map looked sensible!

It takes a long time (an hour or more) to sample points for all the counties,
so to save time, I save the results to a file and 
in future runs, I read the file back in if it exists instead of re-sampling.

```{r county_points, cache = TRUE, dependson="omit_ak_hi"}
points_file <- file.path(data_dir, "county_points.Rds")

if (file.exists(points_file)) {
  points <- read_rds(points_file)
} else {
  d_points <- county_maps %>% 
    st_sample(., size = round(.$democrat / 1000)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "Democratic")
  r_points <- county_maps %>% 
    st_sample(., size = round(.$republican / 1000)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "Republican")
  i_points <- county_maps %>% 
    st_sample(., size = round(.$other / 1000)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "3rd Party")
  points <- bind_rows(d_points, r_points, i_points) %>% 
    mutate(party = factor(party)) %>%
    sample_frac() # scramble the rows so we don't get all the red points 
                  # plotted on top of all the blue points.
  
  write_rds(points, points_file)
}
```

Now that we have the points, we can plot them. I found that plotting 
simple feature points using `geom_sf()` didn't work very well and I could
use the `ggplot2` function `geom_point()` more effectively. However, to do
this I had to turn the `geometry` column of the data frame into regular
_x_ and _y_ coordinates. I did this by calling `as_tibble()` on the
simple feature data frame (see the code above where I generate the points).
This produces a data frame with one row per point, 
where the coordinates of the point are in columns `X` and `Y`.
Now I can use the `ggplot2` `geom_point()` function to plot the points using
`size`, `color`, and `alpha` (transparency) aesthetics.

```{r plot_county_point_density, cache=TRUE, dependson="county_points"}
ggplot(points) + 
  geom_sf(data = state_maps, color = "gray30", fill = NA) +
  geom_point(aes(x = X, y = Y, color = party), size = 0.01, alpha = 0.1) + 
  scale_color_manual(values = c(Democratic = "blue", Republican = "red", 
                                "3rd Party" = "dark green"),
                     name = "Party") +
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
  theme(panel.background = element_rect(fill = "black"),
        legend.background = element_rect(fill = "black", color = "white"),
        legend.text = element_text(color = "white"),
        legend.title = element_text(color = "white"),
        legend.position = c(0,0), legend.justification = c(0,0))
```

This gives a sense of the differing population density for different counties
across the US, but it doesn't accurately reflect the fact that population is not
uniformly distriuted within each county, so we see some weird-looking artifacts
where the population density abruptly changes at the boundary from one county
to another.

We can do a better job if we apply a dasymetric approach: The U.S. Census 
measures the population at a very high spatial resolution and we can use
the census figures to allocate votes within the counties. The census divides
each county into a number of _tracts_, each of which generally has 2000--8000
people. Maps of the tracts and data on their populations and demographic 
characteristics are freely available for download and a number of R packages,
such as `tigris` and `tidycensus`, make it easy to download this data.

`tigris` has functions to download maps of states, counties, census tracts, 
and finer-grained census features such as _block groups_ and _blocks_.
For this project, census tracts will be sufficient.

`tidycensus` has functions to download data from the census on population,
demographic characteristics, etc. for each census tract. 
However, to use these functions you need a key to access the automated census
application programming interface (API) that allows you to automate downloading
data instead of manually clicking on the "[American Fact Finder](https://factfinder.census.gov/faces/nav/jsf/pages/index.xhtml)" 
web site operated by the Census Bureau. 

To get a free key for accessing the Census API, go to the signup page at
<https://api.census.gov/data/key_signup.html> and enter your email and your 
institution (Vanderbilt University). You have to agree to the terms of service, 
which essentially consists of promising to be honest (give credit saying that
your data comes from the Census and don't represent false data as coming from 
the Census) and that you will respect people's privacy by not using this data
to attempt to "hack" the census in order to identify specific individuals.

You will receive your Census API key by email. 
The API key is a character string consisting of 40 hexidecimal digits.
Once you have your key, you can enter it into `tidycensus` with the command 
`census_api_key("<your api key>")`.
You have to do this in each new R session that you start.
You don't want to hard-code your API key into R scripts because then anyone 
you share your scripts with would be able to use your key and you could 
potentially be help responsible for any misconduct that they engage in using
your key.
On a personal computer that you control, you can permanently install your API 
key in your `.RProfile` file so it will be loaded automatically any time you 
start R when you're loggied into your account on that computer.
To do this, call `census_api_key("<your api key>", install=TRUE)`
You can check whether your key was correctly installed by running
`Sys.getenv("CENSUS_API_KEY")` at the R console. This should return a character
string with your API key. If it returns an empty string, you may need to restart
your R session (in RStudio, use the menu command "Session/Restart R" or press
`Ctrl+Shift+F10` on Windows or `Cmd+Shift+F10` on a Mac)

Here is some preliminary housekeeping to make sure the API key has been properly
installed and to set up some options to save tract maps in a local cache so 
you don't have to download them again every time you run this script.

```{r setup_census, cache = TRUE}
options(tigris_use_cache = TRUE, tigris_refresh = FALSE)

if (Sys.getenv("CENSUS_API_KEY") == "") {
  stop('ERROR: You need to get a census API key from ',
       'https://api.census.gov/data/key_signup.html and install it using ',
       'tidycensus::census_api_key("<your key token>").\n\n',
       'On a personal computer that you control, you can permanently install ',
       'a census API key token permanently by calling ',
       'tidycensus::census_api_key("<your key token>", install=TRUE).')
}

tract_map_file  <- file.path(data_dir, "maps", "tract_maps.Rds")
tract_pop_file <- file.path(data_dir, "tract_population.Rds")
```

Now that we have the housekeeping out of the way, download the TIGER tract maps

```{r get_tract_maps, cache = TRUE, message=FALSE, dependson="get_tract_data"}
get_tracts_for_fips <- function(fips) {
  state_fp = str_sub(fips, 1, 2)
  county_fp = str_sub(fips, 3)
  tracts(state_fp, county_fp, cb = TRUE, year = 2016, class = "sf", 
         refresh = FALSE, progress_bar = FALSE)
}

if (file.exists(tract_map_file)) {
  tract_maps <- read_rds(tract_map_file)
} else {
  tract_maps <- map(unique(county_maps$geoid), get_tracts_for_fips) %>%
        { do.call(rbind, .) } %>%
        clean_names()
  write_rds(tract_maps, tract_map_file)
}

# Filter out Alaska and Hawaii
tract_maps <- tract_maps %>% filter(! statefp %in% c("02", "15"))
```

Now we download the population data for each tract and join this with the 
maps to create a simple features data frame with the maps and the population.

```{r get_tract_pop, cache = TRUE, messages = FALSE, dependson = "get_tract_maps"}
fetch_tract_pop <- function(fips, variables = "B01003_001", year = 2016) {
  state_fp <- str_sub(fips, 1, 2)
  county_fp <- str_sub(fips, 3)
  get_acs("tract", variables = variables, year = year, cache_table = TRUE, 
          state = state_fp, county = county_fp, geometry = FALSE)
}

if (file.exists(tract_pop_file)) {
  tract_pop <- read_rds(tract_pop_file)
} else {
  tract_pop <- map_df(county_maps$geoid, fetch_tract_pop) %>%
    clean_names()
  write_rds(tract_pop, tract_pop_file)
}

tract_maps <- tract_maps %>% inner_join(tract_pop, by = "geoid")
```

Now we are ready to generate the dasymetrically sampled points.
We join the tract population maps with the county-level votes.
For each tract, we calcuate the total population of the county and
estimate that the number of votes cast in that tract is the number of 
votes cast in the county times the ratio of the tract population to the 
county population.

If we wanted to get fancier we could also look at the characteristics of the 
tract (land use (urban, suburban, rural) and demographic characteristics such 
as race and median income) and use known correlations of partisan voting 
preferences with land-use and demographic characteristics in order to allocate
the county votes to specific tracts.

Some researchers also combine tract maps with street maps and use the fact that
homes are generally located along streets in order to generate even more precise
dasymetric mapping.

But we're going to keep things simple and just allocate points to census tracts
in proportion to the tracts' populations.

When we are sampling at the tract level, we have to sample more points than one
for every 1000 votes because there are often fewer than 1000 votes per tract.
After some trial and error I figured out that I could sample one point for 
every 10 votes in each county and this would give good coverage to the tracts.
This process takes a **long** time (several hours) because spatial sampling is 
slow and there are almost 73,000 tracts and I am sampling roughly 14 million 
points. 

```{r tract_points, cache = TRUE, dependson="get_tract_pop"}
t_points_file <- file.path(data_dir, "tract_points.Rds")

if (file.exists(t_points_file)) {
  t_points <- read_rds(t_points_file)
} else {
  tm <- tract_maps %>% 
    mutate(FIPS = str_c(statefp, countyfp)) %>% 
    left_join(county_votes, by = "FIPS") %>% 
    group_by(FIPS) %>% 
    mutate(total_pop = sum(estimate)) %>%
    ungroup() %>%
    mutate(pop_weight = estimate / total_pop,
           dem_pop = democrat * pop_weight,
           gop_pop = republican * pop_weight,
           ind_pop = other * pop_weight) %>%
    st_transform("+init=epsg:5070")
  
  dt_points <- tm %>% st_sample(., size = round(.$dem_pop / 10)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "Democratic")
  rt_points <- tm %>% st_sample(., size = round(.$gop_pop / 10)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "Republican")
  it_points <- tm %>% st_sample(., size = round(.$ind_pop / 10)) %>% 
    st_coordinates() %>% as_tibble() %>% mutate(party = "3rd Party")
  
  t_points <- bind_rows(dt_points, rt_points, it_points) %>% 
    mutate(party = factor(party)) %>%
    sample_frac()
  
  write_rds(t_points, t_points_file)
}
```

Now we finally have a data frame with almost 14 million points. 
We don't want to plot 14 million points, so we sub-sample one percent of the 
points to plot. That returns us to the number of points we used in the 
county-level point-density mapping (one percent of one tenth is one in 1000).

```{r plot_tract_point_density, cache=TRUE, dependson="tract_points"}
t_points %>% sample_frac(size = 0.01) %>% 
ggplot() + 
  geom_sf(data = state_maps, color = "gray30", fill = NA) +
  geom_point(aes(x = X, y = Y, color = party), size = 0.01, alpha = 0.1) + 
  scale_color_manual(values = c(Democratic = "blue", Republican = "red", 
                                "3rd Party" = "dark green"),
                     name = "Party") +
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
  theme(panel.background = element_rect(fill = "black"),
        legend.background = element_rect(fill = "black", color = "white"),
        legend.text = element_text(color = "white"),
        legend.title = element_text(color = "white"),
        legend.position = c(0,0), legend.justification = c(0,0))
```

This map now looks like a more realistic depiction of the distribution of 
votes across the United States. However the same caveats that I raised with 
regard to county-level cartograms apply here as well. The presidency is not 
decided by the total number of votes, so maps such as this one can be justly
criticzed as misleading in the context of whether a candidate legitimately 
won the presidency.

There is no objective standard for a "right" way to draw a map. 
Each map has advantages and disadvantages that will make it better-suited for
some purposes and worse for others.

Think critically about maps. Make many maps to illustrate many aspects of a 
complex phenomenon, such as a presidential election. And when you report them,
make sure you clearly tell your audience what each map's strengths and 
shortcomings are in order to help them critically evaluate the arguments 
you are using the maps to support.
