Diversity of fish and vulnerable species in Marine World Heritage Sites based on OBIS data

get_redlist_species <- function() {
  redlist <- tibble()
  page <- 0
  while (TRUE) {
    res <- rl_sp(page, key = "a936c4f78881e79a326e73c4f97f34a6e7d8f9f9e84342bff73c3ceda14992b9")$result
    if (length(res) == 0) {
      break
    }
    redlist <- bind_rows(redlist, res)
    page <- page + 1
  }
  redlist <- redlist %>%
    as_tibble() %>%
    filter(category %in% c("CR", "EN", "VU", "EW", "EX")) %>%
    mutate(category = factor(category, levels = c("EX", "EW", "CR", "EN", "VU"))) %>%
    group_by(scientific_name) %>%
    filter(row_number() == 1) %>%
    ungroup()
  return(redlist)  
}

redlist <- get_redlist_species()

worms <- worms %>%
  left_join(redlist %>% select(scientific_name, category), by = c("scientificName" = "scientific_name"))

openxlsx::write.xlsx(worms, "output/worms_species.xlsx", rowNames = FALSE)
write.csv(worms, "output/worms_species.csv", row.names = FALSE)

# fixes issue with bold package
# invisible(Sys.setlocale("LC_ALL", "C"))

if (!file.exists("data/bold.Rdata")) {
  fish_species <- worms %>% filter(is_fish) %>% pull(scientificName) %>% unique()
  bold_list <- sapply(fish_species, function(x) NULL)
  for (i in 1:length(bold_list)) {
    message(i, " ", names(bold_list)[i])
    if (is.null(bold_list[[i]])) {
      bold_list[[i]] <- tryCatch({
        caspr::bold_statistics(names(bold_list)[i])
      }, warning = function(warning_condition) {
      }, error = function(error_condition) {
      }, finally = {
      })
    }
  }
  save(bold_list, file = "data/bold.Rdata")
} else {
  load("data/bold.Rdata")
}

sequence_numbers <- unlist(map(bold_list, nrow))
fish_sequences <- tibble(species = names(sequence_numbers), sequences = sequence_numbers) %>%
  filter(sequences > 0)

worms <- worms %>%
  left_join(fish_sequences, by = c("scientificName" = "species"))

worms %>%
  summarize(
    species = n(),
    fish = sum(is_fish),
    vulnerable = length(na.omit(category)),
    vulnerable_fish = sum(is_fish * !is.na(category)),
    barcode_fish = sum(sequences > 0, na.rm = T)
  ) %>%
  kable(format.args = list(big.mark = ","))

ggplot(worms %>% filter(!is.na(category))) +
  geom_bar(aes(x = phylum, fill = category)) +
  coord_flip() +
  scale_fill_viridis_d(direction = -1)

if (!file.exists("data/marine_world_heritage.gpkg")) {
  download.file("https://github.com/iobis/mwhs-shapes/blob/master/output/marine_world_heritage.gpkg?raw=true", "data/marine_world_heritage.gpkg")
}

shapes <- st_read("data/marine_world_heritage.gpkg")

shapes_processed <- shapes %>%
  group_by(name) %>%
  summarize()

occ_for_geom <- function(geom) {
  wkt <- st_as_text(st_as_sfc(st_bbox(geom)), digits = 6)
  message(wkt)
  occ <- occurrence(geometry = wkt, fields = c("decimalLongitude", "decimalLatitude", "date_year", "scientificName", "aphiaID", "dataset_id")) %>%
    st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
  occ_filtered <- occ %>%
    filter(st_intersects(geometry, geom, sparse = FALSE)) %>%
    as_tibble() %>%
    select(-geometry)
  return(occ_filtered)
}

sf_use_s2(FALSE)

if (!file.exists("data/occurrence.Rdata")) {
  occs <- map(shapes_processed$geom, occ_for_geom)
  for (i in 1:nrow(shapes_processed)) {
    occs[[i]]$name <- shapes_processed$name[i]
  }
  occ <- bind_rows(occs)
  save(occ, file = "data/occurrence.Rdata")
} else {
  load("data/occurrence.Rdata")
}

occ <- occ %>%
  inner_join(worms %>% select(phylum, class, order, family, genus, is_fish, is_coral, aphiaID, category, sequences), by = "aphiaID") %>%
  mutate(
    aphia_fish = ifelse(is_fish == TRUE, aphiaID, NA),
    aphia_vulnerable = ifelse(!is.na(category), aphiaID, NA),
    aphia_vulnerable_fish = ifelse(is_fish == TRUE & !is.na(category), aphiaID, NA),
    aphia_barcode_fish = ifelse(is_fish == TRUE & sequences > 0, aphiaID, NA),
  ) 

occ %>%
  summarize(
    records = n(),
    species = length(unique(aphiaID)),
    fish = length(unique(aphia_fish)),
    vulnerable = length(unique(aphia_vulnerable)),
    vulnerable_fish = length(unique(aphia_vulnerable_fish)),
    barcode_fish =  length(unique(aphia_barcode_fish))
  ) %>%
  kable(format.args = list(big.mark = ","))

site_stats <- occ %>%
  group_by(name) %>%
  summarize(
    records = n(),
    species = length(unique(aphiaID)),
    fish = length(unique(aphia_fish)),
    vulnerable = length(unique(aphia_vulnerable)),
    vulnerable_fish = length(unique(aphia_vulnerable_fish)),
    barcode_fish =  length(unique(aphia_barcode_fish))
  )

xlsx::write.xlsx(site_stats, file = "output/sites.xlsx")

site_stats %>%
  kable(format.args = list(big.mark = ","))

subset <- site_stats %>%
  filter(name %in% c("Great Barrier Reef", "Belize Barrier Reef Reserve System", "Ha Long Bay", "The Wadden Sea", "Heard and McDonald Islands", "Aldabra Atoll")) %>%
  arrange(name)

subset$x <- c(600, 100, 800, 500, 10, 60)
subset$y <- c(140, 80, 210, 3, 20, 40)

ggplot() +
  geom_segment(subset, mapping = aes(x = x, y = y, xend = fish, yend = vulnerable)) +
  geom_label(subset, mapping = aes(x = x, y = y, label = name), size = 3.5, label.size = 0) +
  geom_point(site_stats, mapping = aes(x = fish, y = vulnerable, size = species, color = species), shape = 21, stroke = 1.2, fill = "white") +
  scale_x_continuous(trans = "log10") +
  scale_y_continuous(trans = "log10") +
  scale_size(range = c(1, 8), breaks = c(1000, 5000, 10000)) +
  scale_color_viridis_c(end = 0.8, trans = "log10", breaks = c(1000, 5000, 10000)) +
  guides(color = guide_legend(), size = guide_legend()) +
  xlab("fish species") + ylab("vulnerable species") +
  ft_theme() +
  ggtitle("Fish and vulnerable species diversity at marine World Heritate sites")

ggsave("output/sites.png", width = 10, height = 6, dpi = 600)

species_lists <- occ %>%
  group_by(name, scientificName, aphiaID, phylum, class, order, family, genus, is_fish, is_coral, category, sequences) %>%
  summarize() %>%
  select(site = name, scientificName, aphiaID, phylum, class, order, family, genus, is_fish, is_coral, iucn_category = category, bold_sequences = sequences) %>%
  arrange(site, phylum, class, order, family, scientificName) %>%
  ungroup()

write.csv(species_lists, file = "output/species_lists.csv", row.names = FALSE, na = "")
xlsx::write.xlsx(as.data.frame(species_lists), file = "output/species_lists.xlsx", showNA = FALSE, row.names = FALSE)
xlsx::write.xlsx(as.data.frame(species_lists %>% filter(is_fish)), file = "output/species_lists_fish.xlsx", showNA = FALSE, row.names = FALSE)