#In this file, we do the coding work to produce "heat maps" like those in the manuscript that show predicted probability of occupancy as a function of predictors at different points inside lakes. Also, we do the work needed to produce maps of the presence/absence data we also show in the manuscript. 

#Load required libraries.
library(tidyverse)
library(sf)
library(viridis)

#First, load the raw data, as we also do in the modeling scripts.
ssw_data = readRDS("ssw_data_final")
#Also load in the psi samples from the posterior from year 8 (the most recent year)
psi.samplesYear8 = readRDS("psi_samplesYR8.rds")

#Take the mean Psi (probability of SSW occupancy) for each point using colMeans. This is to collapse down to an "average" psi for each point across all the posterior samples. We also similiarly get a median Psi. 
meanPsis = colMeans(psi.samplesYear8)
medianPsis = apply(psi.samplesYear8, 2, function(x) {quantile(x, p=c(0.5))})

#Also get the upper and lower bounds opf the 95% Credible Interval. 
lowerboundPsi = apply(psi.samplesYear8, 2, function(x) {quantile(x, p=c(0.025))})
upperboundPsi = apply(psi.samplesYear8, 2, function(x) {quantile(x, p=c(0.975))})

#For graphing, we will want just a single record from each point id, so we will group_by and filter to get just the last record chronologically from each point id. 
ssw_data2 = ssw_data %>% 
  arrange(pt.id, survey_dateStart) %>% 
  group_by(pt.id) %>% 
  filter(survey_dateStart == max(survey_dateStart)) %>% 
  arrange(pt.id)

#Then, we attach to this the mean and median psis as well as the upper and lower bounds from earlier.
ssw_data2$meanPsi = meanPsis
ssw_data2$medianPsi = medianPsis
ssw_data2$lowerboundPsi = lowerboundPsi
ssw_data2$upperboundPsi = upperboundPsi

#Take the difference between the lower and upper bound to represent the uncertainty in the predictions. 
ssw_data2$boundDiff = ssw_data2$upperboundPsi-ssw_data2$lowerboundPsi

#For each lake, we will also calculate the mean occupancy probability in the most recent year of sampling given the psi samples we have from all the points in that lake and year. 

#Create storage objects first. 
lakeMeanPsis = rep(NA, length(unique(ssw_data2$lake)))
lakeUpperPsi = rep(NA, length(unique(ssw_data2$lake)))
lakeLowerPsi = rep(NA, length(unique(ssw_data2$lake)))

#Then, go one by one through the lakes...
lakeNames = unique(ssw_data2$lake)

for(l in 1:length(unique(ssw_data2$lake))) {
  
  #Figure out which point ids are this lake's
  this_lakes_pts = unique(ssw_data2$pt.id[ssw_data2$lake == unique(ssw_data2$lake)[l]])
  
  #Calculate mean psis for each point id first. 
  meanPsis = rowMeans(psi.samplesYear8[,this_lakes_pts])
  
  #Then take the mean of those for the lake as well as the 95% Credible Interval. 
  lakeMeanPsis[l] = mean(meanPsis)
  lakeLowerPsi[l] = unname(quantile(meanPsis, p = 0.025))
  lakeUpperPsi[l] = unname(quantile(meanPsis, p = 0.975))
  
}

#Store the data in a data frame and round the values, then take the difference between the upper and lower credible bounds. 
lakePsiTyping = data.frame(lakeNames, 
                           meanPsi = round(lakeMeanPsis,3), 
                           lowerBound = round(lakeLowerPsi,3), 
                           upperBound = round(lakeUpperPsi,3)) %>% 
  mutate(boundDiff = upperBound-lowerBound)

#Fix depth data to be in meters
ssw_data2 = ssw_data2 %>% mutate(depth_m = depth_ft / 3.281)

#Ok, now to bring in lake polygons and access locations. 
MN_lakes = readRDS("MN_lakes_local") 
WI_lakes = readRDS("WI_lakes_local") 
MN_accesses_local = readRDS("MN_accesses_local") 
WI_accesses_local = readRDS("WI_accesses_local") 

#Put all these data into the most appropriate CRSs. Hard to do for WI's data for some reason. 
ssw_data3MN = st_as_sf(ssw_data2) %>% 
  st_transform(crs=26915) %>% 
  filter(state=="MN")
ssw_data3WI = ssw_data2 %>% 
  filter(state=="WI") %>% 
  st_drop_geometry()
ssw_data3WI = st_as_sf(ssw_data3WI, 
                       coords = c("long", "lat"), 
                       crs=4326) %>% 
  st_transform(crs=26916)

#Rename some of the Wisconsin lake names to look more proper.
WI_lakenames = unique(as.character(ssw_data3WI$lake))
WI_lakenames[7] = "Little Cedar"
WI_lakenames[8] = "Little Muskego"
WI_lakenames[1] = "Big Muskego"
WI_lakenames[10] = 'Lower Nemahbin'

WI_lakes$WBDY_NAME[5] = "Big Muskego Lake"

### The first set of plots will just show median predicted psi values for all points in all lakes, starting with those in Minnesota, based on the calculations done above. However, point size will be mapped to the size of the 95% Credible Interval so that points of greater uncertainty stand out. 

for(lak in 1:length(unique(ssw_data3MN$lake))) { #For each lake...
  
  current.dat = ssw_data3MN %>% 
    filter(lake == unique(ssw_data3MN$lake)[lak]) #Get the data from the current lake we're graphing
  
  #Get the polygon for that lake
  current.lake.poly = MN_lakes %>% 
    filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  
  #Get the accesses on that lake.
  current.accesses = MN_accesses_local %>% 
    filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  
  #Pull the known infested points to plot them distinctly as squares.
  only.1s = filter(current.dat, W == 1)
  
  #Designate the two access types as distinctive for proper plotting.
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Generate the plot.
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") + #Demarkate the polygon edge and fill the lake interior distinctly from the background.
    geom_sf(data = current.dat,  #Add the median psis as colors and map the size of the points to uncertainty.
            aes(geometry = geometry,
                color = medianPsi,
                size = boundDiff), 
            pch = 16) +
    scale_size_continuous("95% CI width", range=c(0.3, 1.3),  #Set custom labels, size ranges, and breaks. 
                          breaks = c(0.1, 0.3, 0.5, 0.7),
                          labels = c("10%", "30%", "50%", "70%")) +
    scale_color_viridis(paste0("\n\n", #Create a programmatic, custom label, pick a specific viridis palette, and set specific limits and breaks. 
                               current.dat$lake[1], 
                               " Lake (",
                               current.dat$state[1], 
                               ")\n", 
                               "\nMedian occupancy\nprobability\n"), 
                        discrete = FALSE,
                        option = "H", 
                        limits = c(0,1), 
                        breaks = c(0,0.2, 0.4, 0.6, 0.8, 1), 
                        labels = c("0%", "20%", "40%", "60%", "80%", "100%")) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) + #Apply red squares in known-occupied sites.
    geom_sf(data = only.1s, size = 1.25, pch = 22, color = "transparent", 
             aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) + #Plot public and private accesses as yellow diamonds.
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray70"), #Custom theming.
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(fill = guide_legend(order = 2, override.aes = list(size = 1.6)), #Customize the legends to all display in a certain manner and order.
           color = guide_colorbar(order = 1, barheight = unit(2, "cm"), override.aes = list(size = 1.6)),
           size = guide_legend(order = 3))
  ggsave(p, file = paste0(current.dat$lake[1], "psi.png"), #Save the graph.
         width = 1400, height = 800, units = "px")
  
}

#Same but now for Wisconsin. Annotations here on out will be minimal. 
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak])
  
  if(lak == 1) { current.dat$lake = "Big Muskego"} #Fix some naming issues unique to WI. 
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  
  current.accesses = WI_accesses_local %>% 
    filter(grepl(unique(ssw_data3WI$lake)[lak], lake) | grepl(WI_lakenames[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")

  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = boundDiff), 
            pch = 16) +
    scale_size_continuous("95% CI width", range=c(0.3, 1.3), 
                          breaks = c(0.1, 0.3, 0.5, 0.7),
                          labels = c("10%", "30%", "50%", "70%")) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), discrete = FALSE,
                        option = "H", limits = c(0,1), breaks = c(0,0.2, 0.4, 0.6, 0.8, 1), labels = c("0%", "20%", "40%", "60%", "80%", "100%")) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.25, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray70"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(fill = guide_legend(order = 2, override.aes = list(size = 1.6)),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm"), override.aes = list(size = 1.6)),
           size = guide_legend(order = 3))
  ggsave(p, file = paste0(current.dat$lake[1], "psi.png"), 
         width = 1400, height = 800, units = "px")
  
}


### The next set of plots is very similar except that we map predictor data to the size of the circles we are graphing for each sample point. 

### FIRST - DEPTH AS THE PREDICTOR, MAPPED TO POINT SIZE. 
for(lak in 1:length(unique(ssw_data3MN$lake))) {
  
  current.dat = ssw_data3MN %>% filter(lake == unique(ssw_data3MN$lake)[lak]) 
  
  current.lake.poly = MN_lakes %>% filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  
  current.accesses = MN_accesses_local %>% filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x depth plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = depth_m), 
            pch = 16) +
    scale_size("Depth (m)",
               breaks = seq(from = round(min(current.dat$depth_m)), 
                            to = round(max(current.dat$depth_m)), length=5),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black")),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm"),),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdepth.png"), 
         width = 1400, height = 800, units = "px")
  
}

#Same but for WI
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak]) 
  
  if(lak == 1) { current.dat$lake = "Big Muskego"} 
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  
  current.accesses = WI_accesses_local %>% filter(grepl(unique(ssw_data3WI$lake)[lak], lake) |
                                                    grepl(WI_lakenames[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x depth plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = depth_m), 
            pch = 16) +
    scale_size("Depth (m)",
               breaks = seq(from = round(min(current.dat$depth_m)), 
                            to = round(max(current.dat$depth_m)), length=5),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black")),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdepth.png"), 
         width = 1400, height = 800, units = "px")
  
}


### FETCH AS THE PREDICTOR, MAPPED TO POINT SIZE. 
for(lak in 1:length(unique(ssw_data3MN$lake))) {
  
  current.dat = ssw_data3MN %>% filter(lake == unique(ssw_data3MN$lake)[lak]) 
  
  current.lake.poly = MN_lakes %>% filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  
  current.accesses = MN_accesses_local %>% filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x fetch plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = max.fetch), 
            pch = 16) +
    scale_size("Fetch (m)",
               breaks = seq(from = round(min(current.dat$max.fetch)), 
                            to = round(max(current.dat$max.fetch)), length=5),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black")),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm"),),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixfetch.png"), 
         width = 1400, height = 800, units = "px")
  
}

#Same but for WI
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak]) 
  
  if(lak == 1) { current.dat$lake = "Big Muskego"} 
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  
  current.accesses = WI_accesses_local %>% filter(grepl(unique(ssw_data3WI$lake)[lak], lake) |
                                                    grepl(WI_lakenames[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x fetch plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = max.fetch), 
            pch = 16) +
    scale_size("Fetch (m)",
               breaks = seq(from = round(min(current.dat$max.fetch)), 
                            to = round(max(current.dat$max.fetch)), length=5),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black")),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixfetch.png"), 
         width = 1400, height = 800, units = "px")
  
}


### NEXT--DEPTHxFETCH AS PREDICTOR MAPPED TO POINT SIZE. 

for(lak in 1:length(unique(ssw_data3MN$lake))) {
  
  current.dat = ssw_data3MN %>% filter(lake == unique(ssw_data3MN$lake)[lak]) 
  
  current.lake.poly = MN_lakes %>% filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  
  current.accesses = MN_accesses_local %>% filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  current.dat$depthfetch = ((current.dat$depth_m - mean(ssw_data2$depth_m))/sd(ssw_data2$depth_m)) *
    ((current.dat$max.fetch - mean(ssw_data2$max.fetch))/sd(ssw_data2$max.fetch))
  
  #Psi x depthfetch plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = depthfetch), 
            pch = 16) +
    scale_size("Depth/Fetch", trans = "reverse",
               breaks = seq(from = min(current.dat$depthfetch), 
                           to = max(current.dat$depthfetch), length=5),
               labels = c("Shallower/Less turbulent", "", "", "","Deeper/more turbulent"),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black", 
                                                   size = c(0.25, 0.525, 0.8, 1.075, 1.35))),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm"),),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdepthfetch.png"), 
         width = 1400, height = 800, units = "px")
  
}

#Same but for WI
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak]) 
  
  if(lak == 1) { current.dat$lake = "Big Muskego"}
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  
  current.accesses = WI_accesses_local %>% filter(grepl(unique(ssw_data3WI$lake)[lak], lake) |
                                                    grepl(WI_lakenames[lak], lake))
  
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  current.dat$depthfetch = ((current.dat$depth_m - mean(ssw_data2$depth_m))/sd(ssw_data2$depth_m)) *
    ((current.dat$max.fetch - mean(ssw_data2$max.fetch))/sd(ssw_data2$max.fetch))
  
  #Psi x depthfetch plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = depthfetch), 
            pch = 16) +
    scale_size("Depth/Fetch", trans = "reverse",
               breaks = seq(from = min(current.dat$depthfetch), 
                            to = max(current.dat$depthfetch), length=5),
               labels = c("Shallower/Less turbulent", "", "", "","Deeper/more turbulent"),
               range = c(0.3, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6), discrete=FALSE) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.2, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3, 
                               override.aes = list(color = "black", 
                                                   size = c(0.25, 0.525, 0.8, 1.075, 1.35))),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdepthfetch.png"), 
         width = 1400, height = 800, units = "px")
  
}


### NEXT--DIST TO NEAREST ACCESS AS PREDICTOR MAPPED TO POINT SIZE.
#FIRST, FOR MN

for(lak in 1:length(unique(ssw_data3MN$lake))) {
  
  current.dat = ssw_data3MN %>% filter(lake == unique(ssw_data3MN$lake)[lak])
  current.lake.poly = MN_lakes %>% filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  current.accesses = MN_accesses_local %>% filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  dists.range = round(seq(from = min(ssw_data2$min.dist), 
                    to = max(ssw_data2$min.dist), length=5))
  
  #Psi x DIST plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = min.dist), 
            pch = 16) +
    scale_size_continuous("Distance to\nnearest access (m)",
               limits = c(0, dists.range[5]),
               range = c(0.2, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6)) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.3, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdist.png"), 
         width = 1400, height = 800, units = "px")
  
}

#Now for WI
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak]) 
  if(lak == 1) { current.dat$lake = "Big Muskego"}
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  current.accesses = WI_accesses_local %>% filter(grepl(unique(ssw_data3WI$lake)[lak], lake) |
                                                grepl(WI_lakenames[lak], lake))
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  dists.range = round(seq(from = min(ssw_data2$min.dist), 
                          to = max(ssw_data2$min.dist), length=5))
  
  #Psi x dist plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = min.dist), 
            pch = 16) +
    scale_size_continuous("Distance to\nnearest access (m)",
                      limits = c(0, dists.range[5]),
                      range = c(0.2, 1.3)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6)) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.3, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixdist.png"), 
         width = 1400, height = 800, units = "px")
  
}


### NEXT--NEAREST ACCESS TYPE AS PREDICTOR MAPPED TO POINT SIZE.
#FIRST, FOR MN

for(lak in 1:length(unique(ssw_data3MN$lake))) {
  
  current.dat = ssw_data3MN %>% filter(lake == unique(ssw_data3MN$lake)[lak])
  current.lake.poly = MN_lakes %>% filter(grepl(unique(ssw_data3MN$lake)[lak], map_label))
  current.accesses = MN_accesses_local %>% filter(grepl(unique(ssw_data3MN$lake)[lak], lake))
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x TYPE plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = type), 
            pch = 16) +
    scale_size_discrete("Nearest access type",range = c(0.5, 1)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6)) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.3, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixtype.png"), 
         width = 1400, height = 800, units = "px")
  
}

#Now for WI
for(lak in 1:length(unique(ssw_data3WI$lake))) {
  
  current.dat = ssw_data3WI %>% filter(lake == unique(ssw_data3WI$lake)[lak]) 
  if(lak == 1) { current.dat$lake = "Big Muskego"} 
  if(lak == 8) { current.dat$lake = "Little Muskego"} 
  if(lak == 10) { current.dat$lake = "Lower Nemahbin"}
  if(lak == 7) { current.dat$lake = "Little Cedar"} 
  current.lake.poly = WI_lakes %>% filter(grepl(WI_lakenames[lak], WBDY_NAME))
  current.accesses = WI_accesses_local %>% filter(grepl(unique(ssw_data3WI$lake)[lak], lake) |
                                                    grepl(WI_lakenames[lak], lake))
  only.1s = filter(current.dat, W == 1)
  public.accesses = filter(current.accesses, type == "PUBLIC")
  private.accesses = filter(current.accesses, type == "PRIVATE")
  
  #Psi x type plots
  p = ggplot() +
    geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
    geom_sf(data = current.dat, 
            aes(geometry = geometry,
                color = medianPsi,
                size = type), 
            pch = 16) +
    scale_size_discrete("Nearest access type",range = c(0.5, 1)) +
    scale_color_viridis(paste0("\n\n", current.dat$lake[1], " Lake (",current.dat$state[1], ")\n", "\nMedian occupancy\nprobability\n"), 
                        option = "H", limits = c(0,1), breaks = seq(from=0, to=1, length=6)) +
    scale_fill_manual("", values = "darkred", labels = c("Known occupied", "")) +
    geom_sf(data = only.1s, size = 1.3, pch = 22, color = "transparent", 
            aes(fill = factor(state))) +
    geom_sf(data = public.accesses, shape = 23,  fill = "#F0B523", size = 2) +
    geom_sf(data = public.accesses, size = 1, shape = "B") +
    geom_sf(data = private.accesses, size = 2, pch = 23, fill = "#F0B523") +
    geom_sf(data = private.accesses, size = 1, shape = "V") +
    theme(panel.grid = element_line(color="gray60"), 
          panel.background = element_rect(fill="white"),
          axis.text = element_text(size=5),
          axis.text.x = element_text(angle=90), 
          legend.text = element_text(size = 6),
          legend.title = element_text(size=6),
          axis.ticks = element_line(color="gray60"),
          legend.key.size = unit(0.2, "cm"),
          legend.key = element_rect(fill="white"),
          legend.spacing = unit(0.001, "cm"), 
          legend.box.just = "left",
          legend.background = element_rect(fill="white")) +
    guides(size = guide_legend(order = 3),
           color = guide_colorbar(order = 1, barheight = unit(2, "cm")),
           fill = guide_legend(order = 2))
  ggsave(p, file = paste0(current.dat$lake[1], "psixtype.png"), 
         width = 1400, height = 800, units = "px")
  
}



### Next, we plot our starry stonewort presence/absence data in map form, such that we can see where SSW is present, where it is absent, and where it is present/absent in the year we are currently focused on but absent/present in previous years.

ssw_dataM = ssw_data #Make a copy of our data object to work with. 

#Break these objects into MN and WI so as to make this a bit easier workflow-wise. 
ssw_dataM_MN = ssw_dataM %>% filter(state == "MN")
ssw_dataM_WI = ssw_dataM %>% 
  filter(state=="WI") %>% 
  st_drop_geometry()
ssw_dataM_WI = st_as_sf(ssw_dataM_WI, 
                       coords = c("long", "lat"), 
                       crs=4326) %>% 
  st_transform(crs=26916)

#Now, for each state, do a for loop that goes lake by lake and then year by year and renders a map like those above. However, this time, at each point, we are marking presence or absence of starry stonewort at each spot, based on whether its been detected there previously or not. 

MN_lakes2for = unique(ssw_dataM_MN$lake) #The lakes to for loop over.

for(f in 1:length(MN_lakes2for)) {

  lake1 = MN_lakes2for[f] #Get current lake
  print(lake1) #Progress bar of sorts
  
  dat1 = filter(ssw_dataM_MN, lake == lake1) %>% 
    arrange(SSW) #Pull all observational data for that lake
  
  years2for = sort(unique(dat1$survey_year)) #Get the unique years that lake was surveyed
  
  for(y in 1:length(years2for)) { #For each year the lake was surveyed
    years1 = years2for[y] #Find current year
    print(years1) #Another progress bar
    
    #Arrange the data so that all SSW == 1s go to the bottom to accommodate multiple surveys at the same lake in a given year (which only occurred a few times). This will put known detections always on top and visible in the graphs we make here. 
    dat2 = filter(dat1, survey_year == years1) %>% 
      select(SSW, geometry, survey_year) %>% 
      arrange(SSW) 
    
    #What we do next will depend a lot on whether this is year 1 for this lake or not...if it's the first year a lake was surveyed, there's no previous data to summarize and incorporate--otherwise, there will be.
    if(y == 1) {
      all.years.dat = dat2
    } else {
      all.years.dat = bind_rows(all.years.dat, dat2)
    }
    dat2$INF.NOW = 0 #Default case

    #Get all previous years' data and find any spots that were infested previously, but did not yield detections in the current survey year.
    prev.years = all.years.dat %>% 
      filter(survey_year != years1, SSW == 1)
    
    if(nrow(prev.years) > 0) { #So long as there were any spots with detections in previous years...
      dat2$INF.NOW[dat2$geometry %in% prev.years$geometry] = 1 #Mark those spots as previously infested but not currently infested thusly. 
    }
    
    dat2$INF.NOW[dat2$SSW == 1] = 2 #All infested spots in the current year, possibly overriding the 1s we just put there (we aren't interested in marking sites that have repeatedly produced detections in our maps).

    #Find any points that were only surveyed in previous years and not surveyed in the current year. We'll use their presence/absence data as our data to plot as white and gray as "old news" in our maps
    old.pts = all.years.dat %>% filter(
      survey_year != years1, #Find points not surveyed in this year but in previous years. 
      !.$geometry %in% dat2$geometry 
      )

    #Create some summary columns in our data structure to reference in our ggplot2 calls. 
    truepos = dat2 %>%  filter(INF.NOW == 2) #Currently infested points
    trueneg = dat2 %>%  filter(INF.NOW == 0) #Currently not infested
    falseneg = dat2 %>% filter(INF.NOW == 1) #Not detected now, but was detected in the past
    oldpos = old.pts %>% filter(SSW == 1) #Not surveyed now, but detected in the past
    oldneg = old.pts %>% filter(SSW == 0) #Not surveyed now, and no past detections.
    
    #Turn the data to the appropriate spatial data in the proper CRS to be plotted. 
    truepos = st_as_sf(truepos) %>% st_transform(crs=26915)
    trueneg = st_as_sf(trueneg) %>% st_transform(crs=26915)
    falseneg = st_as_sf(falseneg) %>% st_transform(crs=26915)
    oldpos = st_as_sf(oldpos) %>% st_transform(crs=26915)
    oldneg = st_as_sf(oldneg) %>% st_transform(crs=26915)
    
    #Create nice labels for our ggplot2 call later.
    if(nrow(truepos) > 0) { truepos$status1 = "Detected" }
    if(nrow(trueneg) > 0) { trueneg$status1 = "Not detected" }
    if(nrow(falseneg) > 0) { falseneg$status1 = "Not detected\nbut prior detection(s)" }
    if(nrow(oldpos) > 0) { oldpos$status2 = "Not surveyed,\nbut prior detection(s)" }
    if(nrow(oldneg) > 0) { oldneg$status2 = "Not surveyed,\nno prior detections" }
    
    #Stitch all these categorical status data into one object.
    alldat = bind_rows(truepos, trueneg, falseneg, oldpos, oldneg)
    #Categorize these variables as factors. 
    alldat$status1 = factor(alldat$status1, 
                           levels = c("Detected",
                                      "Not detected",
                                      "Not detected\nbut prior detection(s)"))
    if(is.null(alldat$status2)) { alldat$status2 = NA }
    alldat$status2 = factor(alldat$status2, 
                            levels = c("Not surveyed,\nbut prior detection(s)", 
                                       "Not surveyed,\nno prior detections"))

    #Get other needed graphic data, such as the lake polygons and the locations of public and private accesses on the lakes. 
    current.lake.poly = MN_lakes %>% filter(grepl(lake1, map_label))
    current.accesses = MN_accesses_local %>% filter(grepl(lake1, lake))
    public.accesses = filter(current.accesses, type == "PUBLIC")
    private.accesses = filter(current.accesses, type == "PRIVATE")

    #The graphics are finicky, especially thanks to the last two categories ("not surveyed now, but..."). So, we use an if/else to make sure these plots render as intended. 
    if(nrow(alldat %>% filter(!is.na(status2))) > 0) {
      
    #Render the ggplot
    p = ggplot() +
            geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
            geom_sf(data = alldat %>% 
                      filter(!is.na(status2)), aes(fill=status2), size = 1, pch = 22) + 
            geom_sf(data = alldat %>% 
                      filter(!is.na(status1)), aes(fill=status1), size = 1, pch = 22) +
            geom_sf(data = public.accesses, shape = 23,  fill = "goldenrod", size = 2) +
            geom_sf(data = public.accesses, size = 1, shape = "B") +
            geom_sf(data = private.accesses, size = 2, pch = 23, fill = "goldenrod") +
            geom_sf(data = private.accesses, size = 1, shape = "V") +
            theme(panel.background = element_rect(fill="white"),
                  axis.text = element_text(size=5),
                  axis.text.x = element_text(angle=90),
                  legend.key = element_rect(fill = "white"),
                  legend.text = element_text(size = 6),
                  legend.title = element_text(size = 6),
                  legend.background = element_rect(fill="white"),
                  axis.ticks = element_line(color="gray60"),
                  panel.grid = element_line(color="gray60"),
                  panel.grid.minor = element_blank()) +
      guides(fill = guide_legend(override.aes = list(size = 5))) +
      scale_fill_discrete(paste0(lake1, " Lake ", "(", years1, ")"), 
                          type = c("Not detected" = turbo(21)[4], 
                                   "Detected" = turbo(21)[21], 
                                   "Not detected\nbut prior detection(s)" = turbo(21)[17], 
                                   "Not surveyed,\nbut prior detection(s)" = "gray65", 
                                   "Not surveyed,\nno prior detections" = "white"),
                          na.value = "white") 
    } else {
      p = ggplot() +
        geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
        geom_sf(data = alldat %>% 
                  filter(!is.na(status1)), aes(fill=status1), size = 1, pch = 22) +
        geom_sf(data = public.accesses, shape = 23,  fill = "goldenrod", size = 2) +
        geom_sf(data = public.accesses, size = 1, shape = "B") +
        geom_sf(data = private.accesses, size = 2, pch = 23, fill = "goldenrod") +
        geom_sf(data = private.accesses, size = 1, shape = "V") +
        theme(panel.background = element_rect(fill="white"),
              axis.text = element_text(size=5),
              axis.text.x = element_text(angle=90),
              legend.key = element_rect(fill = "white"),
              legend.text = element_text(size = 6),
              legend.title = element_text(size = 6),
              legend.background = element_rect(fill="white"),
              axis.ticks = element_line(color="gray60"),
              panel.grid = element_line(color="gray60"),
              panel.grid.minor = element_blank()) +
        guides(fill = guide_legend(override.aes = list(size = 5))) +
        guides(fill = guide_legend(override.aes = list(size = 5))) +
        scale_fill_discrete(paste0(lake1, " Lake ", "(", years1, ")"), 
                            type = c("Not detected" = turbo(21)[4], 
                                     "Detected" = turbo(21)[21], 
                                     "Not detected\nbut prior detection(s)" = turbo(21)[17], 
                                     "Not surveyed,\nbut prior detection(s)" = "gray65", 
                                     "Not surveyed,\nno prior detections" = "white"),
                            na.value = "white") 
    }
    ggsave(p, file = paste0(lake1, years1, ".png"), width = 1400, height = 800, units = "px")
  }
}

#Now, the same series of graphs but for Wisconsin, with some babying of lake names for formal presentation. 
WI_lakes2forA = as.character(unique(ssw_dataM_WI$lake))
WI_lakes2forB = WI_lakes2forA
WI_lakes2forB[3] = "Little Cedar"
WI_lakes2forB[4] = "Little Muskego"
WI_lakes2forB[9] = "Big Muskego"
WI_lakes2forB[12] = "Lower Nemahbin"

for(f in 1:length(WI_lakes2forA)) {
  
  lake1 = WI_lakes2forA[f] 
  lake2 = WI_lakes2forB[f] 
  
  print(lake2)
  
  dat1 = filter(ssw_dataM_WI, lake %in% c(lake1, lake2)) %>% arrange(SSW) 
  
  years2for = sort(unique(dat1$survey_year))
  
  for(y in 1:length(years2for)) { 
    years1 = years2for[y] 
    print(years1)
    dat2 = filter(dat1, survey_year == years1) %>% 
      select(SSW, geometry, survey_year) %>% 
      arrange(SSW)
    
    if(y == 1) {
      all.years.dat = dat2
    } else {
      all.years.dat = bind_rows(all.years.dat, dat2)
    }
    dat2$INF.NOW = 0 
    
    prev.years = all.years.dat %>% 
      filter(survey_year != years1, SSW == 1)
    
    if(nrow(prev.years) > 0) { 
      dat2$INF.NOW[dat2$geometry %in% prev.years$geometry] = 1 
    }
    
    dat2$INF.NOW[dat2$SSW == 1] = 2 
   
    old.pts = all.years.dat %>% filter(
      survey_year != years1, 
      !.$geometry %in% dat2$geometry 
    )
    
    truepos = dat2 %>%  filter(INF.NOW == 2)
    trueneg = dat2 %>%  filter(INF.NOW == 0)
    falseneg = dat2 %>% filter(INF.NOW == 1)
    
    oldpos = old.pts %>% filter(SSW == 1)
    oldneg = old.pts %>% filter(SSW == 0)

    truepos = st_as_sf(truepos, coords = c("long", "lat"), crs = 26916)
    trueneg = st_as_sf(trueneg, coords = c("long", "lat"), crs = 26916)
    falseneg = st_as_sf(falseneg, coords = c("long", "lat"), crs = 26916)
    oldpos = st_as_sf(oldpos, coords = c("long", "lat"), crs = 26916)
    oldneg = st_as_sf(oldneg, coords = c("long", "lat"), crs = 26916)
    
    if(nrow(truepos) > 0) { truepos$status1 = "Detected" }
    if(nrow(trueneg) > 0) { trueneg$status1 = "Not detected" }
    if(nrow(falseneg) > 0) { falseneg$status1 = "Not detected\nbut prior detection(s)" }
    if(nrow(oldpos) > 0) { oldpos$status2 = "Not surveyed,\nbut prior detection(s)" }
    if(nrow(oldneg) > 0) { oldneg$status2 = "Not surveyed,\nno prior detections" }
    
    alldat = bind_rows(truepos, trueneg, falseneg, oldpos, oldneg)
    alldat$status1 = factor(alldat$status1, 
                            levels = c("Detected",
                                       "Not detected",
                                       "Not detected\nbut prior detection(s)"))
    if(is.null(alldat$status2)) { alldat$status2 = NA }
    alldat$status2 = factor(alldat$status2, 
                            levels = c("Not surveyed,\nbut prior detection(s)", 
                                       "Not surveyed,\nno prior detections"))
    
    if(lake2 == "Big Muskego") { 
      current.lake.poly = WI_lakes %>% filter(grepl("Muskego Lake", WBDY_NAME), 
                                              !grepl("Little", WBDY_NAME))
      } else {
    current.lake.poly = WI_lakes %>% filter(grepl(lake2, WBDY_NAME))
      }
    
    current.accesses = WI_accesses_local %>% filter(grepl(lake1, lake) | grepl(lake2, lake))
    
    public.accesses = filter(current.accesses, type == "PUBLIC")
    private.accesses = filter(current.accesses, type == "PRIVATE")
    
    if(nrow(alldat %>% filter(!is.na(status2))) > 0) {
      
      p = ggplot() +
        geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
        geom_sf(data = alldat %>% 
                  filter(!is.na(status2)), aes(fill=status2), size = 1, pch = 22) + 
        geom_sf(data = alldat %>% 
                  filter(!is.na(status1)), aes(fill=status1), size = 1, pch = 22) +
        geom_sf(data = public.accesses, shape = 23,  fill = "goldenrod", size = 2) +
        geom_sf(data = public.accesses, size = 1, shape = "B") +
        geom_sf(data = private.accesses, size = 2, pch = 23, fill = "goldenrod") +
        geom_sf(data = private.accesses, size = 1, shape = "V") +
        theme(panel.background = element_rect(fill="white"),
              axis.text = element_text(size=5),
              axis.text.x = element_text(angle=90),
              legend.key = element_rect(fill = "white"),
              legend.text = element_text(size = 6),
              legend.title = element_text(size = 6),
              legend.background = element_rect(fill="white"),
              axis.ticks = element_line(color="gray60"),
              panel.grid = element_line(color="gray60"),
              panel.grid.minor = element_blank()) +
        guides(fill = guide_legend(override.aes = list(size = 5))) +
        scale_fill_discrete(paste0(lake2, " Lake ", "(", years1, ")"), 
                            type = c("Not detected" = turbo(21)[4], 
                                     "Detected" = turbo(21)[21], 
                                     "Not detected\nbut prior detection(s)" = turbo(21)[17], 
                                     "Not surveyed,\nbut prior detection(s)" = "gray65", 
                                     "Not surveyed,\nno prior detections" = "white"),
                            na.value = "white") 
    } else {
      p = ggplot() +
        geom_sf(data = current.lake.poly, size = 0.6, color = "black", fill = "gray95") +
        geom_sf(data = alldat %>% 
                  filter(!is.na(status1)), aes(fill=status1), size = 1, pch = 22) +
        geom_sf(data = public.accesses, shape = 23,  fill = "goldenrod", size = 2) +
        geom_sf(data = public.accesses, size = 1, shape = "B") +
        geom_sf(data = private.accesses, size = 2, pch = 23, fill = "goldenrod") +
        geom_sf(data = private.accesses, size = 1, shape = "V") +
        theme(panel.background = element_rect(fill="white"),
              axis.text = element_text(size=5),
              axis.text.x = element_text(angle=90),
              legend.key = element_rect(fill = "white"),
              legend.text = element_text(size = 6),
              legend.title = element_text(size = 6),
              legend.background = element_rect(fill="white"),
              axis.ticks = element_line(color="gray60"),
              panel.grid = element_line(color="gray60"),
              panel.grid.minor = element_blank()) +
        guides(fill = guide_legend(override.aes = list(size = 5))) +
        guides(fill = guide_legend(override.aes = list(size = 5))) +
        scale_fill_discrete(paste0(lake2, " Lake ", "(", years1, ")"), 
                            type = c("Not detected" = turbo(21)[4], 
                                     "Detected" = turbo(21)[21], 
                                     "Not detected\nbut prior detection(s)" = turbo(21)[17], 
                                     "Not surveyed,\nbut prior detection(s)" = "gray65", 
                                     "Not surveyed,\nno prior detections" = "white"),
                            na.value = "white") 
    }
    ggsave(p, file = paste0(lake1, years1, ".png"), width = 1400, height = 800, units = "px")
  }
}