Updates on plot code

analysis_paf.R

@@ -0,0 +1,3 @@
+### Proportion of affected species
+
+# prop_critical ~ sites + (1|method) + (1|richness)

codes/analysis_occurrence_data.R

@@ -1,22 +1,21 @@
 #################### eDNA expeditions - scientific analysis ####################
 ########################## Environmental data download #########################
 # January of 2024
-# Author: Mike Burrows
-# Contact: [email protected]
+# Author: Silas Principe, Mike Burrows
+# Contact: [email protected]; [email protected]
 #
-########## Read and analyse eDNA data from UNESCO World Heritage Sites #########
+############## Plot eDNA data from UNESCO World Heritage Sites #################
 
 # Load packages
 library(dplyr)
 library(stringr)
 library(purrr)
 library(arrow)
 library(data.table)
+library(ggplot2)
 fs::dir_create("figures")
 
 
-####
-
 speciestherm<-read_parquet("results/species_tsummaries.parquet")
 
 speciesthermdt<-data.table(speciestherm)
@@ -30,28 +29,18 @@ speciesthermsite<-read_parquet("results/tsummaries_aggregated.parquet")
 
 speciesthermsitedt<-data.table(speciesthermsite)
 
-### CTI by higher geography areas #####
-
-names(speciesthermsitedt)
-# [1] "species"            "AphiaID"            "source_obis"        "source_gbif"        "source_dna"        
-# [6] "higherGeography"    "q_0.5"              "q_0.75"             "q_0.9"              "q_0.95"            
-# [11] "q_0.99"             "q_1"                "mean"               "sd"                 "q_0"               
-# [16] "q_0.01"             "q_0.05"             "q_0.1"              "q_0.25"             "site_current"      
-# [21] "site_ssp126_dec100" "site_ssp126_dec50"  "site_ssp245_dec100" "site_ssp245_dec50"  "site_ssp370_dec100"
-# [26] "site_ssp370_dec50"  "site_ssp460_dec100" "site_ssp460_dec50"  "site_ssp585_dec100" "site_ssp585_dec50" 
-# [31] "records"            "where"
+speciesthermsitedt$where <- ifelse(speciesthermsitedt$where == "Both" | 
+                                     speciesthermsitedt$where == "OBIS/GBIF", "Databases", "eDNA")
 
+### CTI by higher geography areas #####
 ctihighergeog<-speciesthermsitedt[,list(ctiavg=mean(q_0.5),
                                         sdcti=sd(q_0.5),
                                         str=mean(q_0.9-q_0.1),
                                         sstavg=mean(site_current),
                                         nspp=.N),by=c("higherGeography","where")]
 
-# Save for later use
-write_parquet(ctihighergeog, "results/cti_highergeo.parquet")
 
 # Save figures
-pdf("figures/AllSpeciesThermalMetrics.pdf",height=3.5,width=8)
 par(mfrow=c(1,3))
 plot(data=ctihighergeog,ctiavg~sstavg,pch=16,col=unclass(as.factor(where)),cex=1.5,main="All species",
      xlab="Sea surface temperature (°C)", ylab="Community Temperature Index (°C)")
@@ -68,6 +57,53 @@ plot(data=ctihighergeog,str~sstavg,pch=16,col=unclass(as.factor(where)),cex=1.5,
 #abline(a=0,b=1,lty=3)
 legend("bottomleft",legend=levels(as.factor(ctihighergeog$where)),inset=0.02,col=1:3,pt.cex=1.5,pch=16)
 
+
+(p1 <- ggplot(ctihighergeog) +
+  geom_abline(slope = 1, intercept = 0, color = "grey90") +
+  geom_point(aes(x = sstavg, y = ctiavg, colour = where, size = nspp), alpha = .5) +
+  scale_color_manual(values = c("#00c3c9",  "#f58000")) + #"#5151db",
+  scale_y_continuous(limits = c(5, 30)) +
+  scale_x_continuous(limits = c(5, 30)) +
+  scale_size(range = c(2,7)) +
+  guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3))) +
+  ylab("Community Thermal Index (°C)") + 
+  xlab("Sea Surface Temperature (°C)") +
+  theme_light() +
+  theme(panel.grid = element_blank(),
+        panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.3, linetype = 2),
+        legend.title = element_blank()))
+
+(p2 <- ggplot(ctihighergeog) +
+    geom_point(aes(x = sstavg, y = sdcti, colour = where, size = nspp), alpha = .5) +
+    scale_color_manual(values = c("#00c3c9",  "#f58000")) + #"#5151db",
+    # scale_y_continuous(limits = c(5, 30)) +
+    # scale_x_continuous(limits = c(5, 30)) +
+    scale_size(range = c(2,7)) +
+    guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3))) +
+    ylab("SD CTI (°C)") + 
+    xlab("Sea Surface Temperature (°C)") +
+    theme_light() +
+    theme(panel.grid = element_blank(),
+          panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.3, linetype = 2),
+          legend.title = element_blank()))
+
+(p3 <- ggplot(ctihighergeog) +
+    geom_point(aes(x = sstavg, y = str, colour = where, size = nspp), alpha = .5) +
+    scale_color_manual(values = c("#00c3c9",  "#f58000")) + #"#5151db",
+    # scale_y_continuous(limits = c(5, 30)) +
+    # scale_x_continuous(limits = c(5, 30)) +
+    scale_size(range = c(2,7)) +
+    guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3))) +
+    ylab("Species thermal range (average, °C)") + 
+    xlab("Sea Surface Temperature (°C)") +
+    theme_light() +
+    theme(panel.grid = element_blank(),
+          panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.3, linetype = 2),
+          legend.title = element_blank()))
+
+library(patchwork)
+p1 + p2 + p3 + patchwork::plot_layout(guides = "collect")
+
 ### Add taxon information
 
 specieslist<-read_parquet("results/species_list.parquet")

codes/get_sst_marinespecies.R

@@ -8,10 +8,13 @@
 
 # Load packages ----
 library(DBI)
+library(furrr)
+library(storr)
 library(dplyr)
 library(tidyr)
 library(arrow)
 library(terra)
+library(polars)
 fs::dir_create("results")
 
 # Load files/settings ----
@@ -65,15 +68,66 @@ names(env) <- env_vars
 # Get H3 ids from the environmental layer
 biooracle_h3 <- open_dataset("data/biooracle_h3id.parquet")
 
-bio_h3_sel <- biooracle_h3 %>%
-  filter(h3 %in% unique(db_species$h3)) %>%
+# Get unique h3 ids
+un_h3s <- unique(db_species$h3)
+
+# Get cells for those available
+bio_h3_presel <- biooracle_h3 %>%
+  filter(h3 %in% un_h3s) %>%
   select(cell, h3) %>%
   collect()
 
+# Get not available
+not_available <- un_h3s[!un_h3s %in% unique(bio_h3_presel$h3)]
+
+# Get disks
+# Create a temporary folder to hold disks
+fs::dir_create("t_disks")
+
+plan(multisession, workers = 8)
+not_available_h3_b <- split(not_available, as.integer((seq_along(not_available) - 1) / 2000))
+
+un_h3s_disk <- future_map(not_available_h3_b, function(x){
+
+  bio_h3 <- open_dataset("data/biooracle_h3id.parquet")
+
+  disk_cells <- lapply(x, function(z){
+    dc <- unlist(h3jsr::get_disk(z))
+    data.frame(h3_original = z, h3 = dc)
+  })
+  disk_cells <- bind_rows(disk_cells)
+
+  pc <- bio_h3 %>%
+    filter(h3 %in% unique(disk_cells$h3)) %>%
+    select(cell, h3) %>%
+    collect()
+
+  disk_cells <- left_join(disk_cells, pc) %>%
+    filter(!is.na(cell))
+
+  write_parquet(disk_cells,
+                paste0("t_disks/", x[1], ".parquet"))
+  NULL
+}, .progress = T)
+
+# Aggregate expanded disks data
+disk_f <- open_dataset("t_disks/")
+disk_f <- disk_f %>%
+  collect()
+
+# Merge with previous one
+bio_h3_presel <- bio_h3_presel %>%
+  mutate(h3_original = h3) %>%
+  select(h3_original, h3, cell)
+
+bio_h3_sel <- bind_rows(disk_f, bio_h3_presel)
+
+
 # Extract env data
 env_ext <- terra::extract(env, y = as.vector(bio_h3_sel$cell)) %>%
   bind_cols(bio_h3_sel) %>%
-  select(-cell) %>%
+  select(-cell, -h3) %>%
+  rename(h3 = h3_original) %>%
   group_by(h3) %>%
   summarise(across(everything(), ~ mean(.x, na.rm = TRUE)))