le functiones

bcgov · Nov 16, 2024 · fd50a7a · fd50a7a
1 parent 610bde0
commit fd50a7a
Show file tree

Hide file tree

Showing 2 changed files with 264 additions and 3 deletions.
diff --git a/data_processing/Extreme_Vars.qmd b/data_processing/Extreme_Vars.qmd
@@ -138,7 +138,7 @@ ggplot(cmd_res, aes(x = Date, y = cmd_sum)) +
 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)]
+cmd_max_ann[,CMD_3yr := frollsum(CMD_Max, n = 3), by = .(id)] ## maximum value
 
 ## number of months with drought
 cmd_drought <- cmd_res[cmd_sum > 600, .(CMD_Drought = .N), by = .(id)]
@@ -153,12 +153,182 @@ 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)]
+cmd_3yr <- cmd_max_ann[,.(CMD_Max3yr = max(CMD_3yr, na.rm = T)), 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 <- merge(all_normal, cmd_3yr, by = "id", all = T)
+#setnafill(all_normal, fill = 0, cols = names(all_normal[-1]))
 
 all_normal[is.na(all_normal)] <- 0
+
+##test climate future with 3 model runs
+##is the variability increasing (0.01 - 0.99)
+```
+
+```{r}
+
+make_extreme_vars <- function(pnts){
+  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")
+  cmd_vars <- list_vars("Monthly")[grep("CMD|PPT|Tmax",list_vars("Monthly"))]
+  res <- downscale(pnts, which_refmap = "refmap_climr", 
+                 obs_years = 1961:1990, 
+                 obs_ts_dataset = "cru.gpcc", 
+                 vars = c(bio_vars, cmd_vars))
+  res <- res[!is.na(DATASET),]
+  
+  ##quantiles
+  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 vars
+  vars_tmp <- c("id","DATASET","PERIOD",cmd_vars)
+  res2 <- res[,..vars_tmp]
+  res2 <- melt(res2, 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 := cmd_sum(value, cutoff = 500), by = .(id)] ##this line uses the C++ function
+  
+  cmd_res[,Date := as.Date(paste0(PERIOD,"-",Month,"-01"))]
+  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")]
+  
+  ##summarise variables
+  ##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)] ## maximum value
+  
+  ## 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)]
+  cmd_3yr <- cmd_max_ann[,.(CMD_Max3yr = max(CMD_3yr, na.rm = T)), 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)
+  all_normal <- merge(all_normal, cmd_3yr, by = "id", all = T)
+  #setnafill(all_normal, fill = 0, cols = names(all_normal[-1]))
+  
+  all_normal[is.na(all_normal)] <- 0
+  return(all_normal)
+}
+```
+
+```{r}
+make_extreme_vars_fut <- function(pnts, period = 2041:2060, gcm = list_gcms()[1]){
+  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")
+  cmd_vars <- list_vars("Monthly")[grep("CMD|PPT|Tmax",list_vars("Monthly"))]
+  res <- downscale(pnts, which_refmap = "refmap_climr", 
+                 gcms = gcm,
+                 ssps = list_ssps()[2],
+                 gcm_ssp_years = period,
+                 max_run = 3L,
+                 vars = c(bio_vars, cmd_vars))
+  res <- res[!is.na(GCM),][RUN != "ensembleMean",]
+  
+  ##quantiles
+  bio_stats <- res[,lapply(.SD, function(x) quantile(x, c(0.01,0.5,0.99))),
+                 .SDcols = bio_vars, 
+                 by = .(id, GCM, RUN)]
+  bio_stats[,Quantile := rep(c("01","50","99"),nrow(bio_stats)/3)]
+  temp <- melt(bio_stats, id.vars = c("id","GCM","RUN", "Quantile"))
+  bio_stats <- dcast(temp, id + GCM + RUN ~ variable + Quantile)
+  
+  ##CMD vars
+  vars_tmp <- c("id","GCM","RUN","PERIOD",cmd_vars)
+  res2 <- res[,..vars_tmp]
+  res2 <- melt(res2, id.vars = c("id","GCM","RUN","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,Month,GCM, RUN,PERIOD,value)]
+  setorder(res_precip, id,GCM,RUN,PERIOD, Month)
+  precip_recharge <- res_precip[,.(precip_sum = sum(value)), by = .(id,GCM,RUN,PERIOD)]
+  precip_recharge[,Month := 3]
+  
+  cmd_res <- res2[Var == "CMD",.(id,GCM,RUN,PERIOD,Var,Month,value)]
+  setorder(cmd_res, id, GCM,RUN,PERIOD, Month)
+  cmd_res[precip_recharge, precip_sum := precip_sum, on = c("id","GCM","RUN","PERIOD","Month")]
+  cmd_res[!is.na(precip_sum),value := value - precip_sum]
+  cmd_res[,cmd_sum := cmd_sum(value, cutoff = 500), by = .(id,GCM,RUN)] ##this line uses the C++ function
+  
+  cmd_res[,Date := as.Date(paste0(PERIOD,"-",Month,"-01"))]
+  t_res <- res2[Var == "Tmax",.(id,GCM,RUN,PERIOD,Var,Month,value)]
+  setorder(t_res, id,GCM,RUN, PERIOD, Month)
+  cmd_res[t_res, tmax := i.value, on = c("id","GCM","RUN","PERIOD","Month")]
+  
+  ##summarise variables
+  ##yearly CMD_sum max
+  cmd_max_ann <- cmd_res[,.(CMD_Max = max(cmd_sum)), by = .(id,GCM,RUN, PERIOD)]
+  
+  ##three year rolling sum of cmd max
+  cmd_max_ann[,CMD_3yr := frollsum(CMD_Max, n = 3), by = .(id,GCM,RUN)] ## maximum value
+  
+  ## number of months with drought
+  cmd_drought <- cmd_res[cmd_sum > 600, .(CMD_Drought = .N), by = .(id,GCM,RUN)]
+  
+  ##number of months with heat drought
+  hdr_25 <- cmd_res[cmd_sum > 600 & tmax > 25, .(HDR_25 = .N), by = .(id,GCM,RUN)]
+  hdr_30 <- cmd_res[cmd_sum > 600 & tmax > 30, .(HDR_30 = .N), by = .(id,GCM,RUN)]
+  hdr_35 <- cmd_res[cmd_sum > 600 & tmax > 35, .(HDR_35 = .N), by = .(id,GCM,RUN)]
+  
+  ##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,GCM,RUN)]
+  cmd_3yr <- cmd_max_ann[,.(CMD_Max3yr = max(CMD_3yr, na.rm = T)), by = .(id,GCM,RUN)]
+  all_normal <- merge(bio_stats, cmd_max_nrm, by = c("id","GCM","RUN"), all = T)
+  all_normal <- merge(all_normal, hdr_25, by = c("id","GCM","RUN"),all = T)
+  all_normal <- merge(all_normal, hdr_30, by = c("id","GCM","RUN"), all = T)
+  all_normal <- merge(all_normal, hdr_35, by = c("id","GCM","RUN"), all = T)
+  all_normal <- merge(all_normal, cmd_drought, by = c("id","GCM","RUN"), all = T)
+  all_normal <- merge(all_normal, cmd_3yr, by = c("id","GCM","RUN"), all = T)
+  #setnafill(all_normal, fill = 0, cols = names(all_normal[-1]))
+  
+  all_normal[is.na(all_normal)] <- 0
+  return(all_normal)
+}
+```
+
+## Test Functions
+
+```{r}
+dat_fut <- make_extreme_vars_fut(pnts, period = 2041:2060, gcm = list_gcms()[1])
 ```
diff --git a/data_processing/Test_Script.R b/data_processing/Test_Script.R
@@ -1,8 +1,95 @@
 library(data.table)
-
 library(terra)
 library(climr)
+library(RPostgres)
+
+dbCon <- data_connect()
+bbox <- c(55,51.5,-115,-128)
+
+cache_clear("reference")
+climna <- input_refmap(dbCon, bbox)
+plot(climna[[15]])
+
+pnts <- data.table(lon = c(-127.18, -123.46, -125), lat = c(54.89, 48.78, 48.3), elev = c(540,67, 102), id = 1:3)
+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"))
+name <- "normal_composite"
+bands <- 1:73
+
+dbGetTiles <- function(dbCon, name, pnts, rast = "rast", bands = 1:73){
+
+  pnts[,mls := paste0("(",lon," ",lat,")")]
+  wkt_str <- paste0("MULTIPOINT(", paste(pnts$mls, collapse = ", "),")")
+  qry <- paste0("select rid as id,
+            st_xmax(st_envelope(rast)) as xmx,
+            st_xmin(st_envelope(rast)) as xmn,
+            st_ymax(st_envelope(rast)) as ymx,
+            st_ymin(st_envelope(rast)) as ymn,
+            st_width(rast) as cols,
+            st_height(rast) as rows
+            from
+            ", name," 
+            WHERE ST_Intersects(
+              rast,ST_SetSRID(
+                ST_GeomFromText('",wkt_str,"'),
+                4326
+              )
+            )
+            ")
+
+  info <- dbGetQuery(dbCon,qry)
+
+  bandqs1 <- paste0("UNNEST(ST_Dumpvalues(rast, ", bands, ")) as vals_", bands)
+
+  qry2 <- paste0("select rid, ", 
+         paste(bandqs1, collapse = ", "),
+         " from ", name,
+         " where rid IN (",paste(info$id, collapse = ", "),")")
+  r_values <- dbGetQuery(dbCon,qry2)
+
+  rast_vals <- suppressWarnings(melt(r_values, id.vars = c("rid")))
+  out_list <- list()
+  for(tile in unique(info$id)){
+    info_tl <- info[info$id == tile,]
+    rout <- rast(
+      nrows = info_tl$rows, ncols = info_tl$cols, xmin = info_tl$xmn,
+      xmax = info_tl$xmx, ymin = info_tl$ymn, ymax = info_tl$ymx, nlyrs = length(bands),
+      crs = paste0("EPSG:", 4326), vals = rast_vals$value[rast_vals$rid == tile]
+    )
+    out_list[[as.character(tile)]] <- rout
+  }
+
+  rsrc <- sprc(out_list)
+  rast_res <- merge(rsrc)
+  return(rast_res)
+}
+
+out <- dbGetTiles(dbCon, "\"gcm_mpi-esm1-2-hr\"", pnts)
+plot(out[[15]])
+out <- dbGetTiles(dbCon, "normal_composite", pnts)
+plot(out[[15]])
+
 
+varsl = c("DD5", "DDsub0_at", "DDsub0_wt", "PPT_05", "PPT_06", "PPT_07", "PPT_08",
+          "PPT_09", "CMD", "PPT_at", "PPT_wt", "CMD_07", "SHM", "AHM", "NFFD", "PAS", 
+          "CMI", "Tmax_sm", "TD", "PPT_sm", "DD5_sp", "PPT_sp", "PPT_sm", "Tmax_at", "Tmax_wt", "Tmax_sp", "Tmax_sm",
+          "Tmin_at", "Tmin_wt", "Tmin_sp", "Tmin_sm")
+varsl = unique(varsl)
+
+cache_clear()
+load("my_grid.Rdata")
+setDT(my_grid)
+my_grid <- my_grid[lat < 60,]
+clim.bc <- downscale(
+  xyz = my_grid,   which_refmap = "auto",
+  obs_periods = "2001_2020", 
+  vars = varsl)
+
+
+
+bc_pnts <- c(59.9987500001111, 48.3126388889072, -114.057083333295, -138.998750000111)
 
 wlr <- rast("../Common_Files/composite_WNA_dem.tif")
 plot(wlr)
@@ -18,6 +105,10 @@ bbox <- get_bb(my_grid)
 bbox2 <- c(20,14.83,-80,-120)
 refmap <- input_refmap(db, bbox)
 
+db <- data_connect()
+refmap_na <- input_refmap(db, bbox = bc_pnts, reference = "refmap_climr")
+plot(refmap_na$Tmax_07)
+
 plot(refmap$Tmax_07)
 
 pts <- data.frame(lon = c(-124.11, -125.11), lat = rep(48.82, 2), elev = rep(25,2), id = 1:2)