#Read in data

excerpt_data <- read_xlsx('D:/Google Drive/IURIF/Projects/2020_4_COVID-19 Budgets/dedoose/updated_DeDoose_coding.xlsx')
article_data <- read_xlsx('D:/Google Drive/IURIF/Projects/2020_4_COVID-19 Budgets/dedoose/dedoose_articles.xlsx')

us_states <- usa_sf() %>%
  mutate(ID= as.character(name)) %>%
  select(ID, geometry)%>%
  filter(ID != "District of Columbia")

themes_by_state <- article_data %>%
  mutate(
    impact= if_else(impact > 1, 1, impact),
    reaction= if_else(reaction > 1, 1, reaction)) %>% 
  group_by(state) %>%
  summarise(
    impact= sum(impact),
    reaction= sum(reaction),
    total= impact + reaction) %>%
  na_if(0)

theme_map_data <- us_states %>%
  left_join(themes_by_state, by= c('ID'= 'state'))

theme_map_data2 <- theme_map_data %>%
  mutate(impact= replace_na(impact, 0),
         impact_bins= case_when(
           impact== 0 ~ "0",
           impact >0 & impact <= 2 ~ "1 - 2",
           impact > 2 & impact <= 6 ~ "3 - 6",
           impact > 6 & impact <= 12 ~ "7 - 12",
           impact > 12 & impact <= 27 ~ "13 - 27"),
         impact_bins= factor(impact_bins, levels= c("0","1 - 2","3 - 6","7 - 12","13 - 27")))

This document is tested out a bunch of different ways to add shapes to a ggplot map. Overall, the problem is getting the points on the same coordinate system as the map. First, here’s an example using geom_sf and using manually converted centroids in the lat-long format. This does exactly what we want (mulitple shapes, in slightly different positions), the only problem is when we convert the map projection to the Albers projection that we normally use.

states_centers <- as.data.frame(state.center)
states_centers$name <- state.name
coordinates(states_centers) <- ~x+y
proj4string(states_centers) <- CRS(us_longlat_proj)
states_centers <- spTransform(states_centers, CRSobj = CRS(us_aeqd_proj))
states_centers <- as.data.frame(coordinates(states_centers))

centroids <- data.frame(state=(state.name), long=state.center$x, lat=state.center$y)

impact_pal <- c('#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

ggplot()+
  geom_sf(data=theme_map_data2, aes(fill= impact_bins))+
  #coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100") +
  scale_fill_manual(values= impact_pal) +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  theme_minimal() +
 geom_point(data = centroids, aes(color = "blue", x=long+0.5, y=lat), shape=15)+
 geom_point(data = centroids, aes(color = "red", x=long-0.5, y=lat), shape= 16)+
 geom_point(data = centroids, aes(color = "green", x=long, y=lat+0.5), shape= 17)+
 geom_point(data = centroids, aes(color = "yellow", x=long, y=lat-0.5), shape= 18)+
 theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_text(size=16),
        legend.position = "right",
        #legend.position = c(0.87, 0.37),
        legend.text = element_text(size=16),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))

The Albers projection changes the coordinates of the state, but not the points, as shown below.

ggplot()+
  geom_sf(data=theme_map_data2, aes(fill= impact_bins))+
  coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100") +
  scale_fill_manual(values= impact_pal) +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  theme_minimal() +
 geom_point(data = centroids, aes(color = "blue", x=long+0.5, y=lat), shape=15)+
 geom_point(data = centroids, aes(color = "red", x=long-0.5, y=lat), shape= 16)+
 geom_point(data = centroids, aes(color = "green", x=long, y=lat+0.5), shape= 17)+
 geom_point(data = centroids, aes(color = "yellow", x=long, y=lat-0.5), shape= 18)+
 theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_text(size=16),
        legend.position = "right",
        legend.text = element_text(size=16),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))

We can also try use state centers.

ggplot()+
  geom_sf(data=theme_map_data2, aes(fill= impact_bins))+
  coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100") +
  scale_fill_manual(values= impact_pal) +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  theme_minimal() +
 geom_point(data = states_centers, aes(color = "blue", x=x+0.5, y=y), shape=15)+
 geom_point(data = states_centers, aes(color = "red", x=x-0.5, y=y), shape= 16)+
 geom_point(data = states_centers, aes(color = "green", x=x, y=y+0.5), shape= 17)+
 geom_point(data = states_centers, aes(color = "yellow", x=x, y=y-0.5), shape= 18)+
  #stat_sf_coordinates() + #extracts center point, misses michigan
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_text(size=16),
        legend.position = "right",
        legend.text = element_text(size=16),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))

Overall, the problem is getting the centroids to use the same coordinate system. This one works, though!

ggplot() +
  geom_map(
    data = us_data, map = us_data,
    aes(x = long, y = lat, map_id = id, fill = impact_bins),
    color = "#2b2b2b", size = 0.1) +
  scale_fill_manual(values= impact_pal) +
  geom_point(data = states_centers, aes(color = "blue", x=x+(mean(abs(states_centers$x))*0.06), 
                                        y=y), shape=15, size=1.5)+
  geom_point(data = states_centers, aes(color = "red", x=x-(mean(abs(states_centers$x))*0.075),
                                        y=y), shape= 16, size=1.5)+
  geom_point(data = states_centers, aes(color = "green", x=x, y=y+(mean(abs(states_centers$y))*0.075)),
             shape= 17, size=1.5)+
  geom_point(data = states_centers, aes(color = "yellow", x=x, y=y-(mean(abs(states_centers$y))*0.075)), 
             shape= 18, size=1.5)+
  coord_equal() +
  ggthemes::theme_map() +
  labs(
    color= "Impact",
    fill= "Number of Articles")+
  theme(legend.title = element_text(size=16),
        legend.position = "right",
        legend.text = element_text(size=16),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))

Now we’ll use this format to add shape based on whether certain information was collected. (impacts first)

state_revenue <- excerpt_data %>%
  filter(theme=="impact", category=="revenue", !is.na(rev_type)) %>%
  select(state, rev_type) %>%
  distinct() %>%
  mutate(rev_type= ifelse(rev_type=="total"| rev_type=="income"|rev_type=="sales",
                          rev_type, "other"),
         values= 1) %>% #group many sources under "other"
  distinct() %>% #remove multiple others
  pivot_wider(names_from = rev_type, values_from= values, values_fill=0) %>%
  rename(id= state)

states <- us_map %>%
  select(id) %>%
  filter(id != "District of Columbia") %>%
  distinct()

extra_states <- as.data.frame(setdiff(us_map$id, state_revenue$id)) %>%
  rename("id"= "setdiff(us_map$id, state_revenue$id)")

state_revenue2 <- state_revenue %>%
  full_join(extra_states, by="id") %>%
  filter(id != "District of Columbia")

### need to order these states first
ahhhh <- us_map %>%
  select(id) %>% 
  distinct() %>% 
  filter(id != "District of Columbia") %>%
  inner_join(state_revenue2, by="id")

full_data <- states_centers %>%
  cbind(ahhhh) %>%
  replace_na(list(total= 0, income= 0, sales= 0, other= 0))

ggplot() +
  geom_map(
    data = us_data, map = us_data,
    aes(x = long, y = lat, map_id = id, fill = impact_bins),
    color = "#2b2b2b", size = 0.1) +
  scale_fill_manual(values= impact_pal) +
  geom_point(data = full_data, aes(color= "Total Tax", x=x+(mean(abs(x))*0.06), y=y), shape=15, size=2, alpha= 
               ifelse(full_data$total > 0, 1, 0))+
  geom_point(data = full_data, aes(color= "Income Tax",x=x-(mean(abs(x))*0.075), y=y), shape= 16, size=2, alpha= 
               ifelse(full_data$income > 0, 1, 0))+
  geom_point(data = full_data, aes(color= "Sales Tax",x=x, y=y+(mean(abs(y))*0.075)), shape= 17, size=2, alpha= 
               ifelse(full_data$sales > 0, 1, 0))+
  geom_point(data = full_data, aes(color= "Other Tax",x=x, y=y-(mean(abs(y))*0.075)), shape= 18, size=2, alpha= 
               ifelse(full_data$other > 0, 1, 0))+
  coord_equal() +
  ggthemes::theme_map() +
  labs(
    color= "Revenue Impact",
    fill= "Number of Articles")+
  theme(legend.title = element_text(size=15),
        legend.position = c(0.94, 0.15),
        legend.text = element_text(size=15),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))