ulfelder/womens_march_attendance_2017.R

## womens_march_attendance_2017.R

# metadata
data_source <- "https://docs.google.com/spreadsheets/d/1xa0iLqYKz8x9Yc_rfhtmSOJQ2EGgeUVjvV4A8LsIaxY/htmlview?sle=true#gid=0"
data_collectors <- "Jeremy Pressman (@djpressman, U of Connecticut) and\nErica Chenoweth (@EricaChenoweth, U of Denver)"
code_gist <- "https://gist.github.com/benmarwick/a1ac9c7235ebef542824512162ff2f44"


#  ------------------------------------------------------------------------

# read in data from google sheets to get a data frame

library(rvest)
suppressMessages(library(tidyverse))

crowdestimates_raw <-
  read_html(data_source)  %>%
  html_table()
# it's a list, we just want the first item
crowdestimates <- crowdestimates_raw[[1]]
names(crowdestimates) <- make.names(names(crowdestimates), unique = TRUE)


#  ------------------------------------------------------------------------

# clean a few names so that they match with map data sources

crowdestimates$A <- gsub("Oregon", "OR", crowdestimates$A )
crowdestimates$A <- gsub("Oklahoma City", "Oklahoma City", crowdestimates$A )
crowdestimates$A <- gsub("Maine", "ME", crowdestimates$A )
crowdestimates$A <- gsub("Pequannock Township / Pompton Plains, NJ", "Pequannock, NJ", crowdestimates$A )
crowdestimates$A <- gsub("Asbury, Park, NJ", "Asbury Park, NJ", crowdestimates$A )
crowdestimates$A <- gsub("Chillicothe OH", "Chillicothe, OH", crowdestimates$A )
crowdestimates$A <- gsub("Champaign IL", "Champaign, IL", crowdestimates$A )
crowdestimates$A <- gsub("Washington DC" , "Washington, DC" , crowdestimates$A )
crowdestimates$A <- gsub("San Antonio" , "San Antonio, CA" , crowdestimates$A )
crowdestimates$A <- gsub("San Fransisco" , "San Francisco" , crowdestimates$A )
crowdestimates$A <- gsub("St. " , "Saint " , crowdestimates$A )
crowdestimates$A <- gsub("Saint Paul/Minnesota, MN", "Saint Paul, MN" , crowdestimates$A )
crowdestimates$A <- gsub("Saint Paul/Minneapolis, MN", "Saint Paul, MN" , crowdestimates$A )
crowdestimates$A <- gsub("Sault Saint Marie MI", "Sault Sainte Marie, MI" , crowdestimates$A )
crowdestimates$A <- gsub("Orange County", "Santa Ana, CA" , crowdestimates$A )
crowdestimates$A <- gsub("Ft. Worth", "Fort Worth" , crowdestimates$A )
crowdestimates$A <- gsub("Maui", "Lahaina" , crowdestimates$A )
crowdestimates$A <- gsub("Terre Haute, Indiana", "Terre Haute, IN" , crowdestimates$A )
crowdestimates$A <- gsub("Martha's Vineyard", "Edgartown" , crowdestimates$A )
crowdestimates$A <- gsub("Chesapeake Bay", "Havre de Grace" , crowdestimates$A )
crowdestimates$A <- gsub("Kennbunk", "Kennebunk" , crowdestimates$A )
crowdestimates$A <- gsub("San Antonio, CA", "San Antonio, TX" , crowdestimates$A )
crowdestimates$A <- gsub("Ridgeway", "Ridgway" , crowdestimates$A )
crowdestimates$A <- gsub("Melbourne/Brevard County", "Melbourne" , crowdestimates$A )
crowdestimates$A <- gsub("Mt. ", "Mount " , crowdestimates$A )
crowdestimates$A <- gsub("Monterey Bay", "Monterey" , crowdestimates$A )
crowdestimates$A <- gsub("Sand Point", "Sandpoint" , crowdestimates$A )
crowdestimates$A <- gsub("Glen Falls", "Glens Falls" , crowdestimates$A )
crowdestimates$A <- gsub("Pendelton", "Pendleton" , crowdestimates$A )
crowdestimates$A <- gsub("Milhelm", "Millheim" , crowdestimates$A )
crowdestimates$A <- gsub("San Juan Island", "Friday Harbor" , crowdestimates$A )
crowdestimates$A <- gsub("Jackson Hole", "Jackson" , crowdestimates$A )
crowdestimates$A <- gsub("Utqiagvik", "Barrow" , crowdestimates$A )
crowdestimates$A <- gsub("Klamath Falls, Oregon", "Klamath Falls, OR" , crowdestimates$A )


# get city and state in own cols, make the numbers useful
crowdestimates_sep <-
crowdestimates %>%
  separate(A, c("city", "state_abb"), sep = ", ") %>%
  mutate(lower = as.numeric(gsub(",", "", B)),
         upper = as.numeric(gsub(",", "", C)),
         city = tolower(city),
         state_abb = toupper(state_abb)) %>%
  rowwise() %>%
  dplyr::mutate(mean_estimate = mean(c(lower, upper), na.rm = TRUE))

# how many US locations do we have data for?
us_locations <-
  crowdestimates_sep %>%
  filter(state_abb %in% state.abb) %>%
  filter(mean_estimate > 0) %>%
  nrow
 # 359


#  ------------------------------------------------------------------------

# State-level maps

# aggregate by state
crowdestimates_state <-
  crowdestimates_sep %>%
  right_join(data_frame(state_abb = state.abb,
                        state = tolower(state.name))) %>%
  group_by(state_abb, state) %>%
  dplyr::summarise(state_total = sum(mean_estimate, na.rm = TRUE))

# basic map of total counts
library(viridis)
states_map <- map_data("state")

crowdestimates_raw_count_map <-
ggplot(crowdestimates_state,
       aes(map_id = state)) +
  geom_map(aes(fill = state_total),
           map = states_map)  +
  scale_fill_viridis(option = "plasma") +
  expand_limits(x = states_map$long,
                y = states_map$lat) +
  coord_map() +
  theme_bw() +
  ggtitle("Women's March, Jan 2017: \nCounts of people participating by state",
          subtitle = paste0("Data collected by \n",  data_collectors)) +
  labs(caption = paste0("R code online at \n", code_gist))


#  basic map of proprtion of state population
library(choroplethr) # has state population data
data(df_pop_state) # 2012 estimate

crowdestimates_state_prop <-
crowdestimates_state %>%
  left_join(df_pop_state,
            by = c("state" = "region")) %>%
  mutate(mean_estimate_prop = state_total / value)

crowdestimates_prop_count_map <-
  ggplot(crowdestimates_state_prop,
         aes(map_id = state)) +
  geom_map(aes(fill = mean_estimate_prop),
           map = states_map)  +
  scale_fill_viridis(option = "plasma") +
  expand_limits(x = states_map$long,
                y = states_map$lat) +
  theme_bw()  +
  coord_map() +
  ggtitle("Women's March, Jan 2017: \nPropotion of state population participating",
          subtitle = paste0("Data collected by \n",  data_collectors)) +
  labs(caption = paste0("R code online at ", code_gist))


#  ------------------------------------------------------------------------

# County-level maps

# Basic city-level map, using county populations

# here are the county polygons for the map
counties_map <- map_data("county")
counties_map$county <- counties_map$subregion

# get county names and FIPS codes
library(maps)
data(county.fips)

# gives us FIPS, state and county names
county_fips <-
  county.fips %>%
  separate(polyname, c("state", "county"), ",") %>%
  left_join(counties_map,
            by = c("county" = "subregion")) %>%
  distinct(fips, state, county)

# gives us city names, states, & county FIPS for a city (but not county names)
library(noncensus)
data(zip_codes)
zip_codes_county_fips <-
zip_codes %>%
  mutate(state = ifelse(state == "OR" & city == "Walla Walla", "WA", state)) %>%
  mutate(city = tolower(city)) %>%
  dplyr::select(-zip, -latitude, -longitude)

# get county and city names, and county FIPS, and county populations
data(df_county_demographics) # get populations for counties, about 3,144 of them
county_city_fips_map_pops <-
df_county_demographics %>%
  left_join(county_fips,
            by = c("region" = "fips"))  %>%
  left_join(zip_codes_county_fips,
            by = c("region" = "fips")) %>%
  distinct(city, county, state.y, region, total_population) %>%
  mutate(city = ifelse(city == "parcel return service" & state.y == "DC", "washington", city)) %>%
  arrange(city, county)

# join county name and population data to crowd data
crowdestimates_sep_county <-
crowdestimates_sep %>%
  mutate(state = state_abb) %>%
  left_join(county_city_fips_map_pops,
            by = c("city", "state_abb" = "state.y" )) %>%
   mutate(prop_county = mean_estimate / total_population) %>%
  filter(!is.na(mean_estimate)) %>%
  filter(state %in% state.abb)

# some cities join with more than one county
# just get the one with the highest prop
# for simplicity. This could be an error-introducing
# step, probably checking each duplicate by
# hand would be better
crowdestimates_sep_per_county <-
crowdestimates_sep_county %>%
  group_by(city, county, state) %>%
  dplyr::summarise(prop_county = max(prop_county)) %>%
  left_join(crowdestimates_sep_county) %>%
  ungroup() %>%
  mutate(county = gsub(":.*", "", county)) %>%
  filter(!is.na(county))

# join with other counties we don't have data for
# that can show on the map instead of a hole
crowdestimates_sep_per_county_all_counties <-
  crowdestimates_sep_per_county %>%
  right_join(counties_map,
             by = c("county" = "subregion"))


# plot raw counts for counties, we'll use the log values because of the huge range in county-level march attendence.

crowdestimates_county_raw_count_map <-
  ggplot() +
  geom_polygon(data = crowdestimates_sep_per_county_all_counties,
               colour = "grey40",
           aes(long,
               lat,
               group = group,
               fill = log(mean_estimate)))  +
  scale_fill_viridis(option = "plasma") +
  coord_map() +
  theme_bw()  +
  ggtitle("Women's March, Jan 2017: \nCount of participants by county",
          subtitle = paste0("Data collected by \n",  data_collectors)) +
  labs(caption = paste0("R code online at ", code_gist))


# plot proportion of county pop for counties

crowdestimates_county_prop_map <-
  ggplot() +
  geom_polygon(data = crowdestimates_sep_per_county_all_counties,
               colour = "grey40",
               aes(long,
                   lat,
                   group = group,
                   fill = log(prop_county)))  +
  scale_fill_viridis(option = "plasma") +
  coord_map() +
  theme_bw()  +
  ggtitle("Women's March, Jan 2017: \nPropotion of county population participating",
          subtitle = paste0("Data collected by \n",  data_collectors)) +
  labs(caption = paste0("R code online at ", code_gist))


# Copyright 2017 Ben Marwick

# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

	# metadata
	data_source <- "https://docs.google.com/spreadsheets/d/1xa0iLqYKz8x9Yc_rfhtmSOJQ2EGgeUVjvV4A8LsIaxY/htmlview?sle=true#gid=0"
	data_collectors <- "Jeremy Pressman (@djpressman, U of Connecticut) and\nErica Chenoweth (@EricaChenoweth, U of Denver)"
	code_gist <- "https://gist.github.com/benmarwick/a1ac9c7235ebef542824512162ff2f44"


	# ------------------------------------------------------------------------

	# read in data from google sheets to get a data frame

	library(rvest)
	suppressMessages(library(tidyverse))

	crowdestimates_raw <-
	read_html(data_source) %>%
	html_table()
	# it's a list, we just want the first item
	crowdestimates <- crowdestimates_raw[[1]]
	names(crowdestimates) <- make.names(names(crowdestimates), unique = TRUE)



	# ------------------------------------------------------------------------

	# clean a few names so that they match with map data sources

	crowdestimates$A <- gsub("Oregon", "OR", crowdestimates$A )
	crowdestimates$A <- gsub("Oklahoma City", "Oklahoma City", crowdestimates$A )
	crowdestimates$A <- gsub("Maine", "ME", crowdestimates$A )
	crowdestimates$A <- gsub("Pequannock Township / Pompton Plains, NJ", "Pequannock, NJ", crowdestimates$A )
	crowdestimates$A <- gsub("Asbury, Park, NJ", "Asbury Park, NJ", crowdestimates$A )
	crowdestimates$A <- gsub("Chillicothe OH", "Chillicothe, OH", crowdestimates$A )
	crowdestimates$A <- gsub("Champaign IL", "Champaign, IL", crowdestimates$A )
	crowdestimates$A <- gsub("Washington DC" , "Washington, DC" , crowdestimates$A )
	crowdestimates$A <- gsub("San Antonio" , "San Antonio, CA" , crowdestimates$A )
	crowdestimates$A <- gsub("San Fransisco" , "San Francisco" , crowdestimates$A )
	crowdestimates$A <- gsub("St. " , "Saint " , crowdestimates$A )
	crowdestimates$A <- gsub("Saint Paul/Minnesota, MN", "Saint Paul, MN" , crowdestimates$A )
	crowdestimates$A <- gsub("Saint Paul/Minneapolis, MN", "Saint Paul, MN" , crowdestimates$A )
	crowdestimates$A <- gsub("Sault Saint Marie MI", "Sault Sainte Marie, MI" , crowdestimates$A )
	crowdestimates$A <- gsub("Orange County", "Santa Ana, CA" , crowdestimates$A )
	crowdestimates$A <- gsub("Ft. Worth", "Fort Worth" , crowdestimates$A )
	crowdestimates$A <- gsub("Maui", "Lahaina" , crowdestimates$A )
	crowdestimates$A <- gsub("Terre Haute, Indiana", "Terre Haute, IN" , crowdestimates$A )
	crowdestimates$A <- gsub("Martha's Vineyard", "Edgartown" , crowdestimates$A )
	crowdestimates$A <- gsub("Chesapeake Bay", "Havre de Grace" , crowdestimates$A )
	crowdestimates$A <- gsub("Kennbunk", "Kennebunk" , crowdestimates$A )
	crowdestimates$A <- gsub("San Antonio, CA", "San Antonio, TX" , crowdestimates$A )
	crowdestimates$A <- gsub("Ridgeway", "Ridgway" , crowdestimates$A )
	crowdestimates$A <- gsub("Melbourne/Brevard County", "Melbourne" , crowdestimates$A )
	crowdestimates$A <- gsub("Mt. ", "Mount " , crowdestimates$A )
	crowdestimates$A <- gsub("Monterey Bay", "Monterey" , crowdestimates$A )
	crowdestimates$A <- gsub("Sand Point", "Sandpoint" , crowdestimates$A )
	crowdestimates$A <- gsub("Glen Falls", "Glens Falls" , crowdestimates$A )
	crowdestimates$A <- gsub("Pendelton", "Pendleton" , crowdestimates$A )
	crowdestimates$A <- gsub("Milhelm", "Millheim" , crowdestimates$A )
	crowdestimates$A <- gsub("San Juan Island", "Friday Harbor" , crowdestimates$A )
	crowdestimates$A <- gsub("Jackson Hole", "Jackson" , crowdestimates$A )
	crowdestimates$A <- gsub("Utqiagvik", "Barrow" , crowdestimates$A )
	crowdestimates$A <- gsub("Klamath Falls, Oregon", "Klamath Falls, OR" , crowdestimates$A )


	# get city and state in own cols, make the numbers useful
	crowdestimates_sep <-
	crowdestimates %>%
	separate(A, c("city", "state_abb"), sep = ", ") %>%
	mutate(lower = as.numeric(gsub(",", "", B)),
	upper = as.numeric(gsub(",", "", C)),
	city = tolower(city),
	state_abb = toupper(state_abb)) %>%
	rowwise() %>%
	dplyr::mutate(mean_estimate = mean(c(lower, upper), na.rm = TRUE))

	# how many US locations do we have data for?
	us_locations <-
	crowdestimates_sep %>%
	filter(state_abb %in% state.abb) %>%
	filter(mean_estimate > 0) %>%
	nrow
	# 359


	# ------------------------------------------------------------------------

	# State-level maps

	# aggregate by state
	crowdestimates_state <-
	crowdestimates_sep %>%
	right_join(data_frame(state_abb = state.abb,
	state = tolower(state.name))) %>%
	group_by(state_abb, state) %>%
	dplyr::summarise(state_total = sum(mean_estimate, na.rm = TRUE))

	# basic map of total counts
	library(viridis)
	states_map <- map_data("state")

	crowdestimates_raw_count_map <-
	ggplot(crowdestimates_state,
	aes(map_id = state)) +
	geom_map(aes(fill = state_total),
	map = states_map) +
	scale_fill_viridis(option = "plasma") +
	expand_limits(x = states_map$long,
	y = states_map$lat) +
	coord_map() +
	theme_bw() +
	ggtitle("Women's March, Jan 2017: \nCounts of people participating by state",
	subtitle = paste0("Data collected by \n", data_collectors)) +
	labs(caption = paste0("R code online at \n", code_gist))


	# basic map of proprtion of state population
	library(choroplethr) # has state population data
	data(df_pop_state) # 2012 estimate

	crowdestimates_state_prop <-
	crowdestimates_state %>%
	left_join(df_pop_state,
	by = c("state" = "region")) %>%
	mutate(mean_estimate_prop = state_total / value)

	crowdestimates_prop_count_map <-
	ggplot(crowdestimates_state_prop,
	aes(map_id = state)) +
	geom_map(aes(fill = mean_estimate_prop),
	map = states_map) +
	scale_fill_viridis(option = "plasma") +
	expand_limits(x = states_map$long,
	y = states_map$lat) +
	theme_bw() +
	coord_map() +
	ggtitle("Women's March, Jan 2017: \nPropotion of state population participating",
	subtitle = paste0("Data collected by \n", data_collectors)) +
	labs(caption = paste0("R code online at ", code_gist))


	# ------------------------------------------------------------------------

	# County-level maps

	# Basic city-level map, using county populations

	# here are the county polygons for the map
	counties_map <- map_data("county")
	counties_map$county <- counties_map$subregion

	# get county names and FIPS codes
	library(maps)
	data(county.fips)

	# gives us FIPS, state and county names
	county_fips <-
	county.fips %>%
	separate(polyname, c("state", "county"), ",") %>%
	left_join(counties_map,
	by = c("county" = "subregion")) %>%
	distinct(fips, state, county)

	# gives us city names, states, & county FIPS for a city (but not county names)
	library(noncensus)
	data(zip_codes)
	zip_codes_county_fips <-
	zip_codes %>%
	mutate(state = ifelse(state == "OR" & city == "Walla Walla", "WA", state)) %>%
	mutate(city = tolower(city)) %>%
	dplyr::select(-zip, -latitude, -longitude)

	# get county and city names, and county FIPS, and county populations
	data(df_county_demographics) # get populations for counties, about 3,144 of them
	county_city_fips_map_pops <-
	df_county_demographics %>%
	left_join(county_fips,
	by = c("region" = "fips")) %>%
	left_join(zip_codes_county_fips,
	by = c("region" = "fips")) %>%
	distinct(city, county, state.y, region, total_population) %>%
	mutate(city = ifelse(city == "parcel return service" & state.y == "DC", "washington", city)) %>%
	arrange(city, county)

	# join county name and population data to crowd data
	crowdestimates_sep_county <-
	crowdestimates_sep %>%
	mutate(state = state_abb) %>%
	left_join(county_city_fips_map_pops,
	by = c("city", "state_abb" = "state.y" )) %>%
	mutate(prop_county = mean_estimate / total_population) %>%
	filter(!is.na(mean_estimate)) %>%
	filter(state %in% state.abb)

	# some cities join with more than one county
	# just get the one with the highest prop
	# for simplicity. This could be an error-introducing
	# step, probably checking each duplicate by
	# hand would be better
	crowdestimates_sep_per_county <-
	crowdestimates_sep_county %>%
	group_by(city, county, state) %>%
	dplyr::summarise(prop_county = max(prop_county)) %>%
	left_join(crowdestimates_sep_county) %>%
	ungroup() %>%
	mutate(county = gsub(":.*", "", county)) %>%
	filter(!is.na(county))

	# join with other counties we don't have data for
	# that can show on the map instead of a hole
	crowdestimates_sep_per_county_all_counties <-
	crowdestimates_sep_per_county %>%
	right_join(counties_map,
	by = c("county" = "subregion"))


	# plot raw counts for counties, we'll use the log values because of the huge range in county-level march attendence.

	crowdestimates_county_raw_count_map <-
	ggplot() +
	geom_polygon(data = crowdestimates_sep_per_county_all_counties,
	colour = "grey40",
	aes(long,
	lat,
	group = group,
	fill = log(mean_estimate))) +
	scale_fill_viridis(option = "plasma") +
	coord_map() +
	theme_bw() +
	ggtitle("Women's March, Jan 2017: \nCount of participants by county",
	subtitle = paste0("Data collected by \n", data_collectors)) +
	labs(caption = paste0("R code online at ", code_gist))



	# plot proportion of county pop for counties

	crowdestimates_county_prop_map <-
	ggplot() +
	geom_polygon(data = crowdestimates_sep_per_county_all_counties,
	colour = "grey40",
	aes(long,
	lat,
	group = group,
	fill = log(prop_county))) +
	scale_fill_viridis(option = "plasma") +
	coord_map() +
	theme_bw() +
	ggtitle("Women's March, Jan 2017: \nPropotion of county population participating",
	subtitle = paste0("Data collected by \n", data_collectors)) +
	labs(caption = paste0("R code online at ", code_gist))


	# Copyright 2017 Ben Marwick

	# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

	# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

	# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.