OTN and DFO Activity Summary in Eastern Canada: DRAFT

proj_start <- ymd("20190101")
proj_end <- ymd("20230909")
proj_long_upp <- -40.00
proj_long_low <- -70.00
proj_lat_upp <- 60.00
proj_lat_low <- 40.00

url <- "https://gisp.dfo-mpo.gc.ca/arcgis/rest/services/FGP/Oceans_Act_Marine_Protected_Areas/MapServer/0"
MPAs <- get_spatial_layer(url) %>% 
  st_make_valid() %>% 
  st_crop(xmin=proj_long_low,
          ymin=proj_lat_low,
          xmax=proj_long_upp,
          ymax=proj_lat_upp)

geoserver_receivers <- readr::read_csv('https://members.oceantrack.org/geoserver/otn/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=otn:stations_receivers&outputFormat=csv', guess_max = 13579)

geoserver_tag_releases <- readr::read_csv('https://members.oceantrack.org/geoserver/otn/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=otn:animals&outputFormat=csv', guess_max = 13579)

geoserver_tag_releases <- geoserver_tag_releases %>%
  filter(yearcollected > 2018)

geoserver_projects <- readr::read_csv('https://members.oceantrack.org/geoserver/otn/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=otn:otn_resources_metadata_points&outputFormat=csv', guess_max = 13579)

# Filter the geoserver_receivers dataset to select relevant stations 

otn_stations <- geoserver_receivers %>%
  # Remove rows with missing deploy_date values
  filter(!is.na(deploy_date)) %>%
  
  # Select stations that fall within the project timeframe defined above
  filter((deploy_date > proj_start & deploy_date < proj_end) |
           (recovery_date < proj_end & recovery_date > proj_start) |
           (deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(18, 'months')) |
           
           # Select specific models within certain date ranges
           (grepl('VR3', model) & deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(4, 'years')) |
           (grepl('VR4', model) & deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(6, 'years'))) %>%
  
  # Filter stations based on latitude and longitude bounds
  filter(stn_lat >= proj_lat_low & stn_lat <= proj_lat_upp &
           stn_long >= proj_long_low & stn_long <= proj_long_upp)

# Define the columns to keep
columns_to_keep <- c("FID", "institutioncode", "datacenter_reference", "scientificname", "vernacularname","longitude", "latitude", "basisofrecord", "yearcollected", "collector", "classname")

# Filter and select the desired columns
otn_animals <- geoserver_tag_releases %>%
  filter(longitude >= proj_long_low & longitude <= proj_long_upp &
           latitude >= proj_lat_low & latitude <= proj_lat_upp) %>%
  select(all_of(columns_to_keep))

# Define the columns to keep
columns_to_keep_projects <- c("FID", "resource_full_name", "ocean", "seriescode", "status", "collaborationtype")

# Filter and select the desired columns
otn_projects <- geoserver_projects %>%
  filter(seriescode == "DFOCanada") %>%
  select(all_of(columns_to_keep_projects))

# Create a Leaflet map for receiver stations
map <- leaflet(otn_stations) %>%
  # Add default basemap
  addTiles() %>%
  # Add MPA polygons
  addPolygons(data = MPAs, 
              color = "yellow", 
              weight = 0, 
              fillOpacity = 1, 
              label = ~NAME_E) %>%
  
  # Add circle markers for each station
  addCircleMarkers(
    lng = ~stn_long,
    lat = ~stn_lat,
    label = ~paste(seriescode, instrumenttype, sep = " - "),
    radius = 0.5,  # Adjust the circle marker size 
    color = ~ifelse(seriescode == "DFOCanada", "red", "blue"),  # Marker color
    fillOpacity = 0.8  # Adjust fill opacity
  ) %>%
  
  # Add bounding box based on filter
  addRectangles(
    lng1 = proj_long_low,      # Left longitude
    lat1 = proj_lat_low,       # Lower latitude
    lng2 = proj_long_upp,      # Right longitude
    lat2 = proj_lat_upp,       # Upper latitude
    weight = 2.5,              # Border width
    color = "darkorange",      # Border color
    fill = FALSE               # Don't fill the rectangle
  ) %>%
  
  # Add legend for all elements of the plot
  addLegend(
    position = "bottomleft",
    colors = c("blue", "red", "darkorange", "yellow"),
    labels = c("All OTN", "OTN-DFO", "Data filter", "Marine Protected Areas"),
    title = "Legend"
  ) %>%
  
  # Add filter control displaying the range date
  addControl(
    html = sprintf("<strong>Filter:</strong> %s - %s", proj_start, proj_end)
  )

# Calculate the total number of receivers for each seriescode (OTN vs OTN-DFO).
summary_table <- otn_stations %>%
  group_by(seriescode) %>%
  summarize(
    Total_Receivers = n()
  ) %>%
  # Calculate the percentage relative to the total and round it.
  mutate(
    Percentage_Relative_To_Total = paste(round((Total_Receivers / sum(Total_Receivers)) * 100, 0), "%"))

# Rename columns
colnames(summary_table)[colnames(summary_table) == "Total_Receivers"] <- "Total # of Receivers"
colnames(summary_table)[colnames(summary_table) == "Percentage_Relative_To_Total"] <- "% Receivers relative to total"

# Print the summary table as a Markdown table.
table_receivers <- knitr::kable(
  summary_table, 
  format = "markdown",
  align = "l", 
  width = "100%"
)

# Display the map and table
map

table_receivers

# Create a summary table with counts of animals tagged by species.
summary_otn_animals <- otn_animals %>%
  group_by(vernacularname) %>%
  summarise(
    Count = n(),
    Count_DFO = sum(grepl("DFO", institutioncode))
  )

# Reorder the levels of vernacularname based on Count 
summary_otn_animals <- summary_otn_animals %>%
  arrange(desc(Count)) %>%
  mutate(vernacularname = factor(vernacularname, levels = vernacularname))

# Create a bar plot using ggplot2.
plot_animals <- ggplot(summary_otn_animals, aes(x = vernacularname)) +
  geom_bar(aes(y = Count, fill = "OTN"), stat = "identity", position = "dodge") +
  geom_bar(aes(y = Count_DFO, fill = "OTN-DFO"), stat = "identity", position = "dodge") +
  labs(
    title = "Number of Animals Tagged by Species in Eastern Canada",
    x = "Species",
    y = "Number of Animals Tagged"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  coord_flip() +
  scale_fill_manual(values = c("OTN" = "blue", "OTN-DFO" = "red"))

# Convert the ggplot plot to an interactive plot using plotly.
interactive_plot <- ggplotly(plot_animals) %>%
  layout(legend = list(x = 0.5, y = 1))

# Create a Leaflet map to visualize were animals were originally tagged
map_animals <- leaflet(otn_animals) %>%
  addTiles() %>%
  
  # Add MPA polygons
  addPolygons(data = MPAs, 
              color = "yellow", 
              weight = 0, 
              fillOpacity = 1, 
              label = ~NAME_E) %>%

  # Add circle in the map for each animal tagged.
  addCircleMarkers(
    lng = ~longitude,
    lat = ~latitude,
    label = ~paste(institutioncode, scientificname, vernacularname, yearcollected, collector, sep = " - "),
    radius = 0.5,
    color = ~ifelse(grepl("DFO", institutioncode), "red", "blue"),
    fillOpacity = 0.8
  ) %>%

  # Add bounding box based on filter defined at the beginning of code.
  addRectangles(
    lng1 = proj_long_low,
    lat1 = proj_lat_low,
    lng2 = proj_long_upp,
    lat2 = proj_lat_upp,
    weight = 2.5,
    color = "darkorange",
    fill = FALSE
  ) %>%

  # Add legend for different elements of the map.
  addLegend(
    position = "bottomleft",
    colors = c("blue", "red", "darkorange", "yellow"),
    labels = c("All OTN", "OTN-DFO", "Data filter", "Marine Protected Areas"),
    title = "Legend"
  )

# Calculate the minimum and maximum years in the dataset.
min_year <- min(otn_animals$yearcollected)
max_year <- max(otn_animals$yearcollected)

# Add a text box to display the year collected range.
map_animals <- map_animals %>%
  addControl(
    html = sprintf("<strong>Year Collected Range:</strong> %d - %d", min_year, max_year),
    position = "topleft"
  )

interactive_plot

map_animals

# Filter the data frame to include only rows with "DFO" in the institutioncode
# Group the data by collector and institutioncode, and select the first row of each group
filtered_data <- otn_animals %>%
  filter(grepl("DFO", institutioncode)) %>%
  group_by(collector, institutioncode) %>%
  slice(1) %>%
  ungroup()

# Create a table of contributors and tagged species
knitr::kable(filtered_data[, c("collector", "institutioncode", "vernacularname")], caption = "DFO Researchers and Tagged Species (also plotted in the map above in blue)")

# Calculate the percentage of OTN-DFO tagging for each species.
summary_otn_animals$Percentage = ifelse(summary_otn_animals$Count_DFO == 0, "0%", paste0(round(summary_otn_animals$Count / summary_otn_animals$Count_DFO * 100), "%"))

# Create a table with species, OTN total count, and % OTN-DFO tagging.
knitr::kable(
  summary_otn_animals[, c("vernacularname", "Count", "Percentage")],
  format = "markdown", 
  col.names = c("Species", "OTN Total", "% OTN-DFO Tagging"),
  caption = "Summary of Animals Tagged by Species in Eastern Canada"
)

# Create a stacked bar plot of collaboration types by status using ggplot2.
stacked_bar_plot <- ggplot(otn_projects, aes(x = collaborationtype, fill = status)) +
  geom_bar() +  # Create a bar plot.
  labs(title = "Collaboration Type Distribution") +  # Set the plot title.
  theme_minimal() +  # Use a minimal theme.
  scale_fill_manual(values = c("ongoing" = "green", "completed" = "blue", "proposed" = "orange")) +  # Define fill colors.
  coord_flip()  # Flip the coordinates for horizontal bars.

# Convert the ggplot plot to a plotly object.
interactive_stacked_bar_plot <- ggplotly(stacked_bar_plot)

# Customize the layout of the interactive plot.
interactive_stacked_bar_plot %>%
  layout(legend = list(x = 0.5, y = 1))  # Adjust the position of the legend.

# Create an HTML table using knitr::kable to display a summary of DFO-OTN projects.
knitr::kable(
  otn_projects[, c("resource_full_name", "ocean", "status", "collaborationtype")],
  caption = "Summary of DFO-OTN Projects",
  format = "html",
  col.names = c("Project Title", "Ocean", "Status", "Collaboration Type")
)