Add script for creating extreme vars

data_processing/Extreme_Vars.qmd

@@ -0,0 +1,137 @@
+---
+title: "Extreme Climate Variables"
+author: "Kiri Daust"
+format: html
+editor: visual
+---
+
+## Load Data
+
+```{r}
+library(climr)
+library(data.table)
+library(ggplot2)
+
+pnts <- data.table(lon = c(-121.20225,-126.39689,-117.97568,-127.29956,-127.12704,-118.7898,-123.45617,-123.37194), lat = c(50.77239,54.73596,50.28127,50.28127,54.83288,50.24118,48.79616,52.49457), elev = c(588,985,1067,55,563,799,306,1103), id = c("BGxm1","SBSmc2","ICHmw2","CWHvm1","SBSdk","ICHxm1","CDFmm","SBPSdc"))
+
+```
+
+## Downscale
+
+```{r}
+res <- downscale(pnts, which_refmap = "refmap_climr", 
+                 obs_years = 1961:1990, 
+                 obs_ts_dataset = "cru.gpcc", 
+                 vars = list_vars()) #vars = list_vars("Monthly")[grep("CMD|PPT|Tmax",list_vars("Monthly"))]
+
+res <- res[!is.na(DATASET),]
+
+```
+
+## Statistics for Bioclim Vars
+
+```{r}
+bio_vars <- c("MAT","MWMT","MCMT","TD","MAP","MSP","AHM","SHM","DDsub0","DD5","DDsub18","DD18",
+              "NFFD","FFP","bFFP","eFFP","PAS","EMT","EXT","Eref","CMD","RH","CMI",
+              "Tmax_at", "Tmax_sm", "Tmax_sp", "Tmax_wt", 
+              "Tmin_at", "Tmin_sm", "Tmin_sp", "Tmin_wt",
+              "PPT_at", "PPT_sm", "PPT_sp", "PPT_wt")
+
+bio_stats <- res[,lapply(.SD, function(x) quantile(x, c(0.01,0.5,0.99))),
+                 .SDcols = bio_vars, 
+                 by = .(id)]
+bio_stats[,Quantile := rep(c("01","50","99"),nrow(pnts))]
+temp <- melt(bio_stats, id.vars = c("id","Quantile"))
+bio_stats <- dcast(temp, id ~ variable + Quantile)
+
+```
+
+## CMD Variables
+
+```{r}
+calc_cmd <- function(vals){
+  res <- numeric(length(vals))
+  res[1] <- vals[1]
+  for(i in 2:length(vals)){
+    prev <- res[i-1]
+    if(vals[i] < 0){
+      if(prev > 500) prev <- 500
+    }
+    tmp <- prev + vals[i]
+    if(tmp < 0) tmp <- 0
+    res[i] <- tmp
+  }
+  return(res)
+}
+```
+
+```{r}
+res <- downscale(pnts, which_refmap = "refmap_climr", 
+                 obs_years = 1961:1990, 
+                 obs_ts_dataset = "cru.gpcc", 
+                 vars = list_vars("Monthly")[grep("CMD|PPT|Tmax",list_vars("Monthly"))]) #
+res <- res[!is.na(DATASET),]
+res2 <- melt(res, id.vars = c("id","DATASET","PERIOD"))
+res2[,c("Var","Month") := tstrsplit(variable,"_")]
+res2[,Month := as.numeric(Month)]
+res_precip <- res2[Var == "PPT" & Month %in% c(10,11,12,1,2,3),.(id,PERIOD,Month,value)]
+setorder(res_precip, id, PERIOD, Month)
+precip_recharge <- res_precip[,.(precip_sum = sum(value)), by = .(id,PERIOD)]
+precip_recharge[,Month := 3]
+
+cmd_res <- res2[Var == "CMD",.(id,PERIOD,Var,Month,value)]
+setorder(cmd_res, id, PERIOD, Month)
+cmd_res[precip_recharge, precip_sum := precip_sum, on = c("id","PERIOD","Month")]
+cmd_res[!is.na(precip_sum),value := value - precip_sum]
+cmd_res[,cmd_sum := calc_cmd(value), by = .(id)]
+
+cmd_res[,Date := as.Date(paste0(PERIOD,"-",Month,"-01"))]
+
+```
+
+## Plot CMD Sum
+
+```{r}
+t_res <- res2[Var == "Tmax",.(id,PERIOD,Var,Month,value)]
+setorder(t_res, id, PERIOD, Month)
+cmd_res[t_res, tmax := i.value, on = c("id","PERIOD","Month")]
+
+temp <- cmd_res[cmd_sum > 600 & tmax > 30,]
+
+ggplot(cmd_res, aes(x = Date, y = cmd_sum)) +
+  geom_line() +
+  geom_point(data = temp)+
+  facet_wrap(~id, ncol = 1)
+```
+
+## Normal Period Variables
+
+```{r}
+##yearly CMD_sum max
+cmd_max_ann <- cmd_res[,.(CMD_Max = max(cmd_sum)), by = .(id, PERIOD)]
+
+##three year rolling sum of cmd max
+cmd_max_ann[,CMD_3yr := frollsum(CMD_Max, n = 3), by = .(id)]
+
+## number of months with drought
+cmd_drought <- cmd_res[cmd_sum > 600, .(CMD_Drought = .N), by = .(id)]
+
+##number of months with heat drought
+hdr_25 <- cmd_res[cmd_sum > 600 & tmax > 25, .(HDR_25 = .N), by = .(id)]
+hdr_30 <- cmd_res[cmd_sum > 600 & tmax > 30, .(HDR_30 = .N), by = .(id)]
+hdr_35 <- cmd_res[cmd_sum > 600 & tmax > 35, .(HDR_35 = .N), by = .(id)]
+
+##normal period CMD_sum statistics
+cmd_max_nrm <- cmd_max_ann[,.(CMD_Max_01 = quantile(CMD_Max,0.01),
+                              CMD_Max_5 = quantile(CMD_Max,0.5),
+                              CMD_Max_99 = quantile(CMD_Max,0.99)),
+                           by = .(id)]
+all_normal <- merge(bio_stats, cmd_max_nrm, by = "id", all = T)
+all_normal <- merge(all_normal, hdr_25, by = "id",all = T)
+all_normal <- merge(all_normal, hdr_30, by = "id", all = T)
+all_normal <- merge(all_normal, hdr_35, by = "id", all = T)
+all_normal <- merge(all_normal, cmd_drought, by = "id", all = T)
+setnafill(all_normal, fill = 0, cols = names(all_normal[-1]))
+
+all_normal[is.na(all_normal)] <- 0
+```