final_maps.Rmd

---
title: "Aquacutlure Mapping: pie & barcharts"
author: "Maggie Klope"
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'final_maps.html'))})
output: html_document
---

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


```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(ggplot2)
library(patchwork)
library(here)
library(janitor)
library(cowplot)
library(ggpubr)
library(ggrepel)
library(USAboundaries)
library(ggmap)
library(maps)
library(sf)

```



## Using Final Data
```{r message=FALSE, warning=FALSE}
# New Data from Ted 
data_scores <- read_csv(here("data", "Aquaculture_Score_for_Figs.csv")) %>% 
  clean_names() %>% 
  rename(short_names = estuary) %>% 
  rename(estuary = polygon_names)

#spatial data from working group
SNAPP_estuary_points <- read_sf(dsn = here("locations"), layer = "FINAL_SNAPP_ESTUARIES_POINTS-44") %>% 
  dplyr::select(Name, Latitude, Longitude) %>% 
  rename(estuary = Name)

# need to use polygon file to get all of the lat/long for new data set
SNAPP_estuary_polygons <- read_sf(dsn = here("locations"), layer = "FINAL_SNAPP_ESTUARIES_POLYGONS-66") %>%
  rename(estuary = Estuary_Na)

#st_centroid needs projection (using NAD83/California Albers) for st_centroid() to work
SNAPP_estuary_polygons_NAD83 <- SNAPP_estuary_polygons %>%
  st_transform(crs = 3310)

#creating centroid data
SNAPP_estuary_polygons_centroid <- SNAPP_estuary_polygons

SNAPP_estuary_polygons_centroid$centroid <- st_centroid(SNAPP_estuary_polygons_NAD83$geometry) %>%
  st_transform(crs = 4326) #converting to WSG84 to get back in lat/long

#joining new data from Ted with the estuary centroids
estuary_sf <- left_join(data_scores, SNAPP_estuary_polygons_centroid, by = "estuary") %>% 
  dplyr::select(map_number, estuary, ecological_priority, community_restoration, community_harvest, commercial_growers, centroid, geometry)

#filtering new data for sites of ecological priority
high_ecology_points <- estuary_sf %>%
  filter(ecological_priority >= 0.5)


# ancillary map layers 
ca_counties <- read_sf(dsn = here("locations"), layer = "CA_counties")
or_counties <- read_sf(dsn = here("locations"), layer = "OR_counties")
wa_counties <- read_sf(dsn = here("locations"), layer = "WA_counties")
nv_counties <- read_sf(dsn = here("locations"), layer = "NV_counties")

#Reading in shapefiles for Canada and Mexico
canada <- read_sf(dsn = here("locations"), layer = "lpr_000b16a_e")
#st_crs(canada) #EPSG: 9001
# We need to reproject the Canada layer
canada <- st_transform(canada, crs = st_crs(or_counties))

mexico <- read_sf(dsn = here("locations"), layer = "mexstates")
# st_crs(mexico) #EPSG: 4326

```


#### Final Data Bar Chart Map
```{r message=FALSE, warning=FALSE}

coords <- as.data.frame(st_coordinates(SNAPP_estuary_polygons_centroid$centroid)) %>%
  dplyr::mutate(estuary = SNAPP_estuary_polygons_centroid$estuary) %>%
  rename(lat = Y, long = X) %>%
  arrange(desc(lat))

bar_data <- data_scores %>%
  dplyr::filter(ecological_priority >= 0.5) %>%
  left_join(coords, by = "estuary")

coords_high_ecol <- bar_data %>% 
  select(estuary, lat, long)

bar_data_long <- bar_data %>% 
   pivot_longer(-c("map_number", "short_names","estuary", "lat", "long", "ecological_priority"), names_to = "category", values_to = "score")


#' Given a Set of Points and Box sizes, find locations
#' Written by @zachp, updated by @slowkow
findboxes <- function(
  df, xcol, ycol,
  box_padding_x, box_padding_y,
  point_padding_x, point_padding_y,
  xlim, ylim,
  force = 1e-7, maxiter = 20000
) {

  # x and y posiitons as a dataframe
  posdf <- df[c(xcol, ycol)]

  # returnd a df where columns are points
  boxdf <- apply(posdf, 1, function(row) {
    xval <- row[xcol]
    yval <- row[ycol]
    return(c(
      xval - box_padding_x / 2,
      yval - box_padding_y / 2,
      xval + box_padding_x / 2,
      yval + box_padding_y / 2
    ))
  })
  # columns are x1,y1,x2,y2
  boxmatrix <- as.matrix(t(boxdf))

  moved <- ggrepel:::repel_boxes(
    data_points = as.matrix(posdf),
    point_padding_x = point_padding_x,
    point_padding_y = point_padding_y,
    boxes = boxmatrix,
    xlim = xlim,
    ylim = ylim,
    hjust = 0.5,
    vjust = 0.5,
    force = force,
    maxiter = maxiter
  )

  finaldf <- cbind(posdf, moved)
  names(finaldf) <- c("x1", "y1", "x2", "y2")
  return(finaldf)
}

# data frame with two x and two y values
df1 <- findboxes(coords_high_ecol,
  xcol = "long", ycol = "lat",
  box_padding_x = Reduce("-", rev(range(coords$lat))) * .30,
  box_padding_y = Reduce("-", rev(range(coords$long))) * .125,
  point_padding_x = Reduce("-", rev(range(coords$long))) * .30,
  point_padding_y = Reduce("-", rev(range(coords$lat))) * .125,
  force = 1e-4, #force = repulsion
  xlim = c(-120, -112),
  ylim = c(24, 55)
)

#renaming first column to match original bar_data df
colnames(df1)[1] <- "long"

#adding offset values to bar_data
bar_data_offset <- coords_high_ecol %>%
  left_join(df1, by = "long") %>%
  select(-c(y1)) %>% #removing repeated lat values
  rename(lat_offset = y2, long_offset = x2) #renaming to be more descriptive

# taking offset lat long and creating a geometry
offset_points <- st_as_sf(bar_data_offset, crs = 4326, coords = c("long_offset", "lat_offset")) %>%
  select(estuary, geometry)

# make a separate ggplot for each estuary:
test_plots <- bar_data_long %>% 
  group_by(estuary) %>% 
  nest() %>%
  mutate(plot = map2(data, estuary,
                     ~ggplot(data = .x) +
                       aes(x = category, y = score) +
                       geom_bar(stat = "identity", show.legend = FALSE, aes(fill = category)) +
                       ylim(0, 1)+
                       theme(axis.title.x=element_blank(), #removing all labels
                             axis.text.x=element_blank(),
                             axis.ticks.x=element_blank(),
                             axis.title.y=element_blank(),
                             axis.text.y=element_blank(),
                             axis.ticks.y=element_blank(),
                             panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #removing grids
                             panel.background = element_rect(fill = "transparent",colour = NA), #making background transparent
                             plot.background = element_rect(fill = "transparent",colour = NA),
                             plot.margin = unit(c(0, 0, 0, 0), "null"),
                             axis.line = element_line())
           )
         )

# turn the plots into annotations: 
estuary_annotation <- offset_points %>% 
  bind_cols(as_tibble(st_coordinates(.))) %>% #retrieves coordinates in matrix form as a tibble, and binds into one dataframe
  st_drop_geometry() %>% #removes geometry from sf object
  left_join(test_plots, by = "estuary") %>% 
  mutate(annotation = pmap(list( X, Y, plot), #creating an annotation column, then using the coordinates and plots in the annotation_custom() function
                           ~annotation_custom(ggplotGrob(..3), # ..3 refers to the third argument, so plot 
                                              xmin = ..1 - 0, xmax = ..1 + 3, # ..1 refers to X
                                              ymin = ..2 - 0, ymax = ..2 + 2))) %>% # ..2 refers to Y
  pull(annotation)

#saving legend as separate object
simple_barchart <-ggplot(bar_data_long, aes(score, fill = category))+
         geom_bar()+
  theme(legend.margin = margin(0, 0, 0, 0))+
  theme(
    legend.title = element_text(size = 7),
    legend.text = element_text(size = 7),
    legend.background = element_rect(fill = NA))+
  scale_fill_discrete(name = "Score", labels = c("Commercial Production", "Community Engagement Harvest", "Community Engagement Restoration"))

legend <- get_legend(simple_barchart)
gg_legend <-as_ggplot(legend)

#create final map
bar_charts_offset <- ggplot(offset_points) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_sf(data = mexico) + 
  xlim(c(-130, -113))+
  ylim(c(24, 55))+
  geom_segment(data = bar_data_offset, aes(x = long, y = lat, xend = long_offset, yend = lat_offset)) +
  geom_point(data = bar_data_offset, aes(long, lat), color = "black") +
  geom_point(data = bar_data_offset, aes(long_offset, lat_offset), color = "black") +
  # geom_bar(data = bar_data_long, aes(score, fill = category))+
  estuary_annotation +
  theme_void()

# bar_charts_offset + (plot_spacer() / (gg_legend + plot_spacer() + plot_spacer())) +plot_layout(nrow=1)

layout <- c(
  area(t = 1, l = 1, b = 9, r = 7),
  area(t = 8, l = 1, b = 8, r = 6)
)

final_bar_chart_map <- bar_charts_offset + gg_legend + 
  plot_layout(design = layout)

```

```{r}
ggsave("figures/final_bar_chart_map.png", final_bar_chart_map, width = 12, height = 6, dpi = 300)
final_bar_chart_map
```

#### Final Data Point Map

```{r, final point map, warning=FALSE, message = FALSE}

centroid_points <-  SNAPP_estuary_polygons_centroid %>%
  select(estuary, centroid) %>%
  right_join(data_scores, by = "estuary") %>% 
  mutate(high_ecol = ifelse(ecological_priority >= 0.5, "yes", "no"))

Zoom_point_1 <- ggplot() +
 geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = centroid_points$centroid, aes(color = centroid_points$ecological_priority, size = 2, shape = centroid_points$high_ecol)) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  scale_shape_manual(values=c(20, 18))+
  scale_shape_discrete(name = "Ecological Score", labels = c("> 0.5", "< 0.5"))+
  labs(color = "Ecological Score", shape = "Ecology Score") + #changing names on legend
  coord_sf(xlim = c(-125.5, -121), ylim = c(46, 49.5)) +
  guides(size = FALSE)+ #removing size from legend
  theme_void()+
  theme(legend.position = "none")

Zoom_point_2 <- ggplot() +
 geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = centroid_points$centroid, aes(color = centroid_points$ecological_priority, size = 2, shape = centroid_points$high_ecol)) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  scale_shape_manual(values=c(20, 18))+
  scale_shape_discrete(name = "Ecological Score", labels = c("> 0.5", "< 0.5"))+
  labs(color = "Ecological Score", shape = "Ecology Score") + #changing names on legend
  coord_sf(xlim = c(-125.5, -121), ylim = c(39, 46)) +
  guides(size = FALSE)+ #removing size from legend
  theme_void()+
  theme(legend.position = "none")

Zoom_point_3 <- ggplot() +
 geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = centroid_points$centroid, aes(color = centroid_points$ecological_priority, size = 2, shape = centroid_points$high_ecol)) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  scale_shape_manual(values=c(20, 18))+
  scale_shape_discrete(name = "Ecological Score", labels = c("> 0.5", "< 0.5"))+
  labs(color = "Ecological Score", shape = "Ecology Score") + #changing names on legend
  coord_sf(xlim = c(-123.5, -119), ylim = c(34.85, 39)) +
  guides(size = FALSE)+ #removing size from legend
  theme_void()+
  theme(legend.position = "none")

Zoom_point_4 <- ggplot() +
 geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = centroid_points$centroid, aes(color = centroid_points$ecological_priority, size = 2, shape = centroid_points$high_ecol)) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  scale_shape_manual(values=c(20, 18))+
  scale_shape_discrete(name = "Ecological Score", labels = c("> 0.5", "< 0.5"))+
  labs(color = "Ecological Score", shape = "Ecology Score") + #changing names on legend
  coord_sf(xlim = c(-120, -112), ylim = c(24, 34.85)) +
  guides(size = FALSE)+ #removing size from legend
  theme_void()+
  theme(legend.position = "none")

point_map <- ggplot() +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) +
  geom_sf(data = mexico) +
  geom_sf(data = centroid_points$centroid, aes(color = centroid_points$ecological_priority, size = 2, shape = centroid_points$high_ecol)) +
  scale_color_gradientn(colours = c(
    "#eff3ff",
    "#6baed6",
    "#084594"
  )) +
  scale_shape_manual(values=c(20, 18))+
  scale_shape_discrete(name = "Ecological Score", labels = c("> 0.5", "< 0.5"))+
  labs(color = "Ecological Score", shape = "Ecology Score") + #changing names on legend
  coord_sf(xlim = c(-130, -112), ylim = c(24, 53)) +
  scale_x_continuous(breaks = c(-130, -116)) +
  scale_y_continuous(breaks = c(30, 40, 50)) +
  guides(size = FALSE)+ #removing size from legend
  theme_void()

legend_3 <- get_legend(point_map)
legend_3

point_composite_map <- (Zoom_point_1) + (Zoom_point_2) + (Zoom_point_3) + (Zoom_point_4) + (legend_3) +
  plot_layout(nrow = 1, ncol = 5)


```


```{r}
ggsave("figures/point_map.png", point_map, width = 12, height = 6, dpi = 300)
point_map
ggsave("figures/point_map_composite.png", point_composite_map, width = 12, height = 6, dpi = 300)
point_composite_map 
```

#### Final Data Bar Chart Supplement

```{r}
bar_data_40_long <-  centroid_points %>% 
  st_drop_geometry() %>% 
  select(map_number, short_names, ecological_priority, community_restoration, community_harvest, commercial_growers) %>% 
  pivot_longer(-c(map_number, short_names), names_to = "category", values_to = "score")

all_sites_wrap <- ggplot(bar_data_40_long, aes(x = category, y = score)) +
  geom_bar(stat = "identity", show.legend = FALSE, aes(fill = category)) +
  # ggtitle(graph_title) + #trying to add different title to each
  ylim(0,1)+
  # theme(plot.title = element_text(size = 6, face = "bold"))+
  facet_wrap(~map_number, ncol = 10, nrow = 4)+
  theme(axis.title.x = element_blank(), #removing all labels
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        # if (i != 1) {axis.title.y = element_blank()},
        # axis.title.y =  ifelse(i == 1, element_blank(), element_text(size = rel(1.5), angle = 90)),
        # axis.title.y = element_blank(),
        # axis.text.y = element_blank(),
        # axis.ticks.y = element_blank(),
        #strip.text.x = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #removing grids
        panel.background = element_rect(fill = "transparent",colour = NA), #making background transparent
        plot.background = element_rect(fill = "transparent",colour = NA),
        plot.margin = unit(c(0, 0, 0, 0), "null"),
        axis.line = element_line(),
        strip.text = element_text(size = 7),#changing title font size for face_wrap()
        axis.text.y = element_text(size = 6)) + #changing y-axis font size
  coord_fixed(ratio = 5)

zoom_data <- data_scores %>%
  left_join(coords, by = "estuary")

base_map_zoom_1_data <- zoom_data %>% 
  dplyr::filter(long > -125.5 & long < -121) %>% 
  dplyr::filter(lat > 40 & lat < 49.5 )

base_map_zoom_2_data <- zoom_data %>% 
dplyr::filter(long > -123.5 & long < -119.5) %>% 
  dplyr::filter(lat > 35 & lat < 39 )

base_map_zoom_3_data <- zoom_data %>% 
  dplyr::filter(long < -112 & long > -120) %>% 
  dplyr::filter(lat < 35 & lat > 24)

base_map_zoom_4_data <- zoom_data %>% 
  dplyr::filter(long < -112 & long > -120) %>% 
  dplyr::filter(lat < 35 & lat > 31)

base_map_zoom_5_data <- zoom_data %>% 
  dplyr::filter(long < -112 & long > -120) %>% 
  dplyr::filter(lat < 31 & lat > 24)

base_map_zoom_1 <- ggplot(base_map_zoom_1_data) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  # geom_point(aes(long, lat), color = "black")+
  coord_sf(xlim = c(-125.5, -121), ylim = c(40, 49.5))+
  # geom_text(aes(long, lat, label = base_map_zoom_1_data$ID, color = "blue"))+
  geom_text_repel(data = base_map_zoom_1_data, aes(long, lat), 
                  label = base_map_zoom_1_data$map_number, 
                  color = "red",
                  #nudge_x = -126, 
                  size = 3,
                  min.segment.length = 0,
                  point.padding = NA,
                  box.padding = 0,
                  segment.color = "black"
                 )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) #changing angle of x-axis

base_map_zoom_2 <- ggplot(base_map_zoom_2_data) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  # geom_point(aes(long, lat), color = "black")+
  coord_sf(xlim = c(-123.5, -119.5), ylim = c(35, 39))+
  # geom_text(aes(long, lat, label = base_map_zoom_2_data$ID, color = "blue"))+
  geom_text_repel(data = base_map_zoom_2_data, aes(long, lat), 
                  label = base_map_zoom_2_data$map_number, 
                  color = "red",
                  #nudge_x = -126, 
                  size = 3,
                  min.segment.length = 0,
                  point.padding = NA,
                  box.padding = 0,
                  segment.color = "black"
                 )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) #changing angle of x-axis)

base_map_zoom_3 <- ggplot(base_map_zoom_3_data) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_sf(data = mexico) +
  coord_sf(xlim = c(-120, -112), ylim = c(24, 35))+
  geom_text_repel(data = base_map_zoom_3_data, aes(long, lat), 
                  label = base_map_zoom_3_data$map_number, 
                  color = "red",
                  #nudge_x = -126, 
                  size = 3,
                  min.segment.length = 0,
                  point.padding = NA,
                  box.padding = 0,
                  segment.color = "black"
                 )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) #changing angle of x-axis)

base_map_zoom_4 <- ggplot(base_map_zoom_4_data) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_sf(data = mexico) +
  coord_sf(xlim = c(-120, -112), ylim = c(31, 35))+
  geom_text_repel(data = base_map_zoom_4_data, aes(long, lat), 
                  label = base_map_zoom_4_data$map_number, 
                  color = "red",
                  #nudge_x = -126, 
                  size = 3,
                  min.segment.length = 0,
                  point.padding = NA,
                  box.padding = 0,
                  segment.color = "black"
                 )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) #changing angle of x-axis)

base_map_zoom_5 <- ggplot(base_map_zoom_5_data) +
  geom_sf(data = us_boundaries()) +
  geom_sf(data = canada) + 
  geom_sf(data = mexico) +
  coord_sf(xlim = c(-120, -112), ylim = c(24, 31))+
  geom_text_repel(data = base_map_zoom_5_data, aes(long, lat), 
                  label = base_map_zoom_5_data$map_number, 
                  color = "red",
                  #nudge_x = -126, 
                  size = 3,
                  min.segment.length = 0,
                  point.padding = NA,
                  box.padding = 0,
                  segment.color = "black"
                 )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) #changing angle of x-axis)

barchart_supplement_v1 <- (all_sites_wrap) + (base_map_zoom_1) + (base_map_zoom_2) + (base_map_zoom_3)+
  plot_layout(nrow = 1, ncol = 4, widths = c(4, 1, 1, 1))

barchart_supplement_v2 <- (all_sites_wrap) + (base_map_zoom_1) + (base_map_zoom_2 / base_map_zoom_4) + (base_map_zoom_5)+
  plot_layout(nrow = 1, ncol = 5, widths = c(4, 1, 1, 1, 1))

prelim_barchart_supplement <- (all_sites_wrap) + (base_map_zoom_1) + (base_map_zoom_2 / base_map_zoom_3)

```

```{r}
ggsave("figures/barchart_supplement_v1.png", barchart_supplement_v1, width = 12, height = 6, dpi = 300)
ggsave("figures/barchart_supplement_v2.png", barchart_supplement_v2, width = 12, height = 6, dpi = 300)
barchart_supplement_v1
barchart_supplement_v2
```