The new maps that have been created using the polygon data start on Line 416

#This code chunk is going to be where I dissolve the counties of each state so that I am left with only the outline of the state.

ca_counties$area <- st_area(ca_counties)
# ca_counties

ca <- ca_counties %>%
  summarise(area = sum(area))

or_counties$area <- st_area(or_counties)
or <- or_counties %>%
  summarise(area = sum(area))

wa_counties$area <- st_area(wa_counties)
wa <- wa_counties %>%
  summarise(area = sum(area))

nv_counties$area <- st_area(nv_counties)
nv <- nv_counties %>%
  summarise(area = sum(area))

Cleaning and wrangling the score data

#Cleaning and wrangling the data

data_scores <- full_join(data_goals, scores_clean) %>%
  clean_names()


estuary_sf <- data_scores %>%
  drop_na("long") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  clean_names()
estuary_reactive_format <- estuary_sf %>%
  gather(score_type, score, -estuary_or_subbasin, -geometry)

# I have put all of the scores into a single column. This should help with the interactive map or if I make a shiny app of the map.

estuary_reactive_lat_lon <- data_scores %>%
  drop_na("long") %>%
  gather(score_type, score, -estuary_or_subbasin, - long, -lat)
#Making an interactive map with tmap

snapp_estuary_map <- tm_shape(estuary_sf) +
  tm_dots(labels = "estuary_or_subbasin", col = "green", size = 0.1)

basemap <- tm_basemap("Esri.WorldImagery")

tmap_mode("view")
snapp_estuary_map +
  basemap

Interactive Map of Estuaries: Hovering the mouse over an estuary will show the estuary’s name. Clicking on an estuary will show the estuaries scores.

The code below will be using the estuary shapefiles provided by the project.

# Making individual map of the ecology estuary scores

SNAPP_estuary_ecology_polygons <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_polygons, aes(fill = Ecol1), lwd = 0) +
  scale_fill_gradientn(colors = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_void()

# SNAPP_estuary_ecology_polygons

SNAPP_estuary_ecology_points <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Ecol1), size = 3) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_void()

# SNAPP_estuary_ecology_points
#Making individual map of the harvest estuary scores

SNAPP_estuary_harvest_polygons <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_polygons, aes(fill = Harvest1), lwd = 0) +
  scale_fill_gradientn(
    colors = c(
    "#005a32",
    "#74c476",
    "#c7e9c0"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_classic()
# SNAPP_estuary_harvest_polygons

SNAPP_estuary_harvest_points <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Harvest1), size = 3) +
  scale_color_gradientn(
    colors = c(
    "#005a32",
    "#74c476",
    "#c7e9c0"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_classic()


# SNAPP_estuary_harvest_points
# Making individual map of the restoration estuary scores

SNAPP_estuary_restoration_polygons <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_polygons, aes(fill = Restor1), lwd = 0) +
  scale_fill_gradientn(colors = c(
    "violet",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme(panel.background = element_blank())

# SNAPP_estuary_restoration_polygons

# Points
SNAPP_estuary_restoration_points <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Resto1), size = 3) +
  scale_color_gradientn(colors = c(
    "violet",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme(panel.background = element_blank())


# SNAPP_estuary_restoration_points
#Making individual map of the community esutary scores


SNAPP_estuary_community_polygons <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_polygons, aes(fill = Comm1), lwd = 0) +
  scale_fill_gradientn(colors = c(
    "green",
    "blue",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()

SNAPP_estuary_community_points <- ggplot() +
  geom_sf(data = us_boundaries())  +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Comm1), size = 3) +
  scale_color_gradientn(colours = c(
    "green",
    "blue",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()
  

# SNAPP_estuary_community_polygons
# SNAPP_estuary_community_points

Estuaries

Polygons based

Dots based

#Here I am going to be working on displaying multiple attributes on one map I will be using points for these maps

ecology_harvest_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Harvest1, size = Ecol1)) +  
  scale_color_gradientn(
    colors = c(
    "#005a32",
    "#74c476",
    "#c7e9c0"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()

ecology_restoration_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Resto1, size = Ecol1)) +  
  scale_color_gradientn(colors = c(
    "violet",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) 

ecology_comm_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = SNAPP_estuary_points, aes(color = Comm1, size = Ecol1)) +  
  scale_color_gradientn(colors = c(
    "green",
    "blue",
    "purple"
  )) +
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) 

ecology_harvest_map

ecology_restoration_map

ecology_comm_map

#Here I will focus on maps with ecology scores over 0.5

high_ecology_polygons <- SNAPP_estuary_polygons %>%
  filter(Ecol1 > 0.5)

high_ecology_points <- SNAPP_estuary_points %>%
  filter(Ecol1 >= 0.5)

high_ecology_harvest_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_points, aes(color = Harvest1), size = 4) +  
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()

high_ecology_restoration_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_points, aes(color = Resto1), size = 4) +  
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()

high_ecology_comm_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_points, aes(color = Comm1), size = 4) +  
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()


# high_ecology_harvest_map
# high_ecology_restoration_map
# high_ecology_comm_map
#This will be for a close up of estuarys with high ecology score

Zoom_high_ecology1 <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = high_ecology_polygons, fill = "red", color = "red") +  
  coord_sf(xlim = c(-125, -121.5), ylim = c(45, 49.5)) +
  theme_bw()

Zoom_high_ecology2 <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = high_ecology_polygons, fill = "red", color = "red") +  
  coord_sf(xlim = c(-123.5, -118.5), ylim = c(34, 38.5)) +
  theme_bw()


Zoom_high_ecology3 <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_polygons, fill = "red", color = "red") +  
  coord_sf(xlim = c(-120, -116), ylim = c(30, 34.5)) +
  theme_bw()


Zoom_high_ecology1

Zoom_high_ecology2

Zoom_high_ecology3

#this will be trying to make the close up maps with the point data

Zoom_high_ecology1_points <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = high_ecology_points, color = "red", size = 3) +  
  coord_sf(xlim = c(-125, -121.5), ylim = c(45, 49.5)) +
  theme_bw()

Zoom_high_ecology2_points <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = high_ecology_points, fill = "red", color = "red", size = 3) +  
  coord_sf(xlim = c(-123.5, -118.5), ylim = c(34, 38.5)) +
  theme_bw()


Zoom_high_ecology3_points <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_points, fill = "red", color = "red", size = 3) +  
  coord_sf(xlim = c(-120, -116), ylim = c(30, 34.5)) +
  theme_bw()


Zoom_high_ecology1_points

Zoom_high_ecology2_points

Zoom_high_ecology3_points

#Using scatterpie to put piecharts over map

#need wide format data
pie_data <- high_ecology_points %>%
  select(-NCEASmap) %>%
  rename(ecology = Ecol1, restoration = Resto1, harvest = Harvest1, commercial = Comm1) %>%
  as.data.frame()
#using the scores from the 10 estuary/subbasins with radius proportional to Ecological score
pie_data$radius <- pie_data$ecology/2.5

Zoom_pie_1 <-  ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                  cols = c("restoration", "harvest", "commercial"), color=NA, alpha=.8, show.legend = FALSE) +
  # geom_scatterpie_legend(pie_data$radius, x=-121.5, y=49) +
  coord_sf(xlim = c(-125.5, -121), ylim = c(45, 49.5)) +
  theme_bw()


Zoom_pie_2 <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                    cols = c("restoration", "harvest", "commercial"), color=NA, alpha=.8, show.legend = FALSE) +
  # geom_scatterpie_legend(pie_data$radius, x=-119.5, y=38) +
  coord_sf(xlim = c(-123.5, -119), ylim = c(34.85, 38.5)) +
  theme_bw()

Zoom_pie_3 <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = mexico) +
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                    cols = c("restoration", "harvest", "commercial"), color=NA, alpha=.8, legend_name = "Score") +
  geom_scatterpie_legend(pie_data$radius, x=-119.5, y=30.5, ) +
  coord_sf(xlim = c(-120, -115), ylim = c(30, 34.5)) +
  theme_bw()


# put it together
Zoom_pie_ecol_rad <- Zoom_pie_1 + Zoom_pie_2 + Zoom_pie_3
ggsave("figures/final_map_pie_ecol_rad.png", Zoom_pie_ecol_rad, width = 12, height = 6, dpi = 300)

#using a fixed radius

pie_data$radius <- .225

Zoom_pie_1_fix <-  ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                  cols = c("ecology", "restoration", "harvest", "commercial"), color=NA, alpha=.8, show.legend = FALSE) +
  coord_sf(xlim = c(-125.5, -120), ylim = c(45, 49.5)) +
  theme_bw()


Zoom_pie_2_fix <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                    cols = c("ecology", "restoration", "harvest", "commercial"), color=NA, alpha=.8, show.legend = FALSE) +
  coord_sf(xlim = c(-123.5, -119), ylim = c(34.8, 38.5)) +
  theme_bw()

Zoom_pie_3_fix <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = mexico) +
  geom_scatterpie(aes(x = Longitude, y = Latitude, group = Name, r = radius), data = pie_data, 
                    cols = c("ecology", "restoration", "harvest", "commercial"), color=NA, alpha=.8, legend_name = "Score") +
  coord_sf(xlim = c(-120, -115), ylim = c(30, 34.5)) +
  theme_bw()

# put it together
Zoom_pie_ecol_rad_fixed <- Zoom_pie_1_fix + Zoom_pie_2_fix + Zoom_pie_3_fix
ggsave("figures/final_map_pie_ecol_rad_fixed.png", Zoom_pie_ecol_rad_fixed, width = 12, height = 6, dpi = 300)

#testing changing alpha to better view overlapping sites

high_ecology_rest_alpha <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = high_ecology_points, aes(color = Resto1), size = 4, alpha = 0.75) +  
  coord_sf(xlim = c(-130, -115), ylim = c(30, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  theme_bw()

high_ecology_rest_alpha