adding habitat data

habitats/air_tmp_vulnerability.Rmd

@@ -0,0 +1,119 @@
+---
+title: "Untitled"
+output: html_document
+date: "2024-09-12"
+---
+
+Going to use current pressure data and habitat distributions to determine how vulnerable some key habitats are to the air temperature pressure.
+
+```{r setup, include=FALSE}
+library(tidyverse)
+library(terra)
+library(here)
+library(countrycode)
+
+```
+
+```{r}
+heat <- rast("/home/shares/ohi/stressors_2021/combining_pressures/rescaled_pressures/air-heat-index_ssp245_current.tif")
+seagrass <- rast(here("habitats/data/seagrass.tif"))
+
+stack <- c(heat, seagrass)
+stack_df <- data.frame(stack)
+
+seagrass <- stack_df %>%
+  filter(!is.na(mean)) %>%
+  mutate(seagrass = ifelse(is.na(seagrass), 0, seagrass))
+
+Hmisc::wtd.quantile(seagrass$mean, weights=seagrass$seagrass, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+
+tmp <-  seagrass %>%
+  group_by(seagrass) %>%
+  mutate(mean = ifelse(mean >0.21, 1, 0)) %>%
+  mutate(seagrass = ifelse(seagrass > 0, 1, 0)) 
+table(tmp)
+
+# > 9855/205281
+#[1] 0.04800737
+#> 254*0.048
+#[1] 12.192
+#> (1-(9855/205281))^254
+#[1] 3.740515e-06
+
+ggplot(seagrass, aes(x=mean, y=seagrass)) +
+  geom_jitter(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("seagrass cover") +
+  xlab("heat-index pressure")
+
+
+
+mangrove <- rast(here("habitats/data/mangroves.tif"))
+
+stack <- c(heat, mangrove)
+stack_df <- data.frame(stack)
+
+mangrove <- stack_df %>%
+  filter(!is.na(mean)) %>%
+  mutate(mangroves = ifelse(is.na(mangroves), 0, mangroves))
+
+
+ggplot(mangrove, aes(x=mean, y=mangroves)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("mangrove cover") +
+  xlab("heat-index pressure")
+
+
+ice <- rast(here("habitats/data/ice.tif"))
+
+stack <- c(heat, ice)
+stack_df <- data.frame(stack)
+
+ice <- stack_df %>%
+  filter(!is.na(mean)) %>%
+  mutate(ice = ifelse(is.na(ice), 0, ice))
+
+Hmisc::wtd.quantile(ice$mean, weights=ice$ice, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+mean(ice$mean[ice$ice>0])
+ggplot(ice, aes(x=mean, y=ice)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red")+ 
+  ylab("seaice cover") +
+  xlab("heat-index pressure")
+
+
+marsh <- rast(here("habitats/data/salt-marsh.tif"))
+
+stack <- c(heat, marsh)
+stack_df <- data.frame(stack)
+
+marsh <- stack_df %>%
+  filter(!is.na(mean)) %>%
+  mutate(salt_marsh = ifelse(is.na(salt_marsh), 0, salt_marsh))
+
+Hmisc::wtd.quantile(marsh$mean, weights=marsh$salt_marsh, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+
+tmp <-  marsh %>%
+  group_by(salt_marsh) %>%
+  mutate(mean = ifelse(mean >0.20, 1, 0)) %>%
+  mutate(salt_marsh = ifelse(salt_marsh > 0, 1, 0)) 
+table(tmp)
+#> 4997/210048
+#[1] 0.0237898
+#> 345*0.0237898
+#[1] 8.207481
+#> (1 - 4997/210048)^345
+#[1] 0.0002468611
+ 
+ggplot(marsh, aes(x=mean, y=salt_marsh)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("saltmarsh cover") +
+  xlab("heat-index pressure")
+
+```

habitats/sst_vulnerability.Rmd

@@ -0,0 +1,96 @@
+---
+title: "Untitled"
+output: html_document
+date: "2024-09-12"
+---
+
+Going to use current pressure data and habitat distributions to determine how vulnerable some key habitats are to the air temperature pressure.
+
+```{r setup, include=FALSE}
+library(tidyverse)
+library(terra)
+library(here)
+library(countrycode)
+
+```
+
+```{r}
+sst <- rast("/home/shares/ohi/stressors_2021/_dataprep/SST/five_year_average_final/sst_avg_ssp245_2015_2019.tif")
+sst <- sst-273.15
+
+kelp <- rast(here("habitats/data/kelp.tif"))
+
+stack <- c(sst, kelp)
+stack_df <- data.frame(stack)
+
+kelp <- stack_df %>%
+  filter(!is.na(focal_mean)) %>%
+  mutate(kelp = ifelse(is.na(kelp), 0, kelp))
+
+Hmisc::wtd.quantile(kelp$focal_mean, weights=kelp$kelp,probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+
+ggplot(kelp, aes(x=focal_mean, y=kelp)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("kelp") +
+  xlab("SST")
+
+
+coral <- rast(here("habitats/data/coral-reef.tif"))
+
+stack <- c(sst, coral)
+stack_df <- data.frame(stack)
+
+coral <- stack_df %>%
+  filter(!is.na(focal_mean)) %>%
+  mutate(coral_reef = ifelse(is.na(coral_reef), 0, coral_reef))
+
+Hmisc::wtd.quantile(coral$focal_mean, weights=coral$coral_reef, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.95, 0.99, 0.999))
+
+ggplot(coral, aes(x=focal_mean, y=coral_reef)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("coral") +
+  xlab("SST")
+
+
+marsh <- rast(here("habitats/data/salt-marsh.tif")) 
+
+stack <- c(sst, marsh) 
+
+stack_df <- data.frame(stack)
+
+marsh <- stack_df %>%
+  filter(!is.na(focal_mean)) %>%
+  mutate(salt_marsh = ifelse(is.na(salt_marsh), 0, salt_marsh))
+
+Hmisc::wtd.quantile(marsh$focal_mean, weights=marsh$salt_marsh, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+
+ggplot(marsh, aes(x=focal_mean, y=salt_marsh)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("marsh") +
+  xlab("SST")
+
+
+seagrass <- rast(here("habitats/data/seagrass.tif")) 
+
+stack <- c(sst, seagrass) 
+
+stack_df <- data.frame(stack)
+
+seagrass <- stack_df %>%
+  filter(!is.na(focal_mean)) %>%
+  mutate(seagrass = ifelse(is.na(seagrass), 0, seagrass))
+
+Hmisc::wtd.quantile(seagrass$focal_mean, weights=seagrass$seagrass, probs= c(0.001, 0.01, 0.1, 0.5, 0.90, 0.99, 0.999))
+
+ggplot(seagrass, aes(x=focal_mean, y=seagrass)) +
+  geom_point(size=0.2, alpha=0.5, height=0.1) +
+  geom_smooth(method = "gam", se = TRUE) + 
+  geom_hline(yintercept=0, color="red") + 
+  ylab("seagrass") +
+  xlab("SST")