fix multiplicative factor when imposing climate changes

Deltares-research · Jan 23, 2024 · 4242a52 · 4242a52
1 parent 80ca637
commit 4242a52
Show file tree

Hide file tree

Showing 3 changed files with 7,750 additions and 20,876 deletions.
diff --git a/R/imposeClimateChanges.R b/R/imposeClimateChanges.R
@@ -33,44 +33,55 @@ imposeClimateChanges <- function(
 
  {
 
+  # Number of grid cells
   ngrids <- length(climate.data)
 
+  # Year and month indices
   year_vec <- as.numeric(format(sim.dates,"%Y"))
   year_ind <- year_vec - min(year_vec) + 1
   month_ind <- as.numeric(format(sim.dates,"%m"))
 
-  # Define daily temperature change factors
+
+  ### Define Daily temperature change factors
+
+  # Transient changes (gradual)
   if(isTRUE(transient.temp.change)) {
       tempf1 <- sapply(1:12, function(x)
-          seq(0, change.factor.temp.mean[x], length.out = max(year_ind)))
-  } else {
+          seq(0, change.factor.temp.mean[x] * 2, length.out = max(year_ind)))
+
+      # Step changes
+      } else {
       tempf1 <- sapply(1:12, function(x)
-          rep(change.factor.temp.mean[x]/2, length.out = max(year_ind)))
+          rep(change.factor.temp.mean[x], length.out = max(year_ind)))
   }
 
-  tempf2 <- sapply(1:length(sim.dates), function(x) tempf1[year_ind[x], month_ind[x]])
+  tempf2 <- sapply(1:length(sim.dates),
+                   function(x) tempf1[year_ind[x], month_ind[x]])
 
-  # Define daily precipitation change factors
+  ### Define Daily precipitation change factors
+
+  # Transient changes (gradual)
   if(isTRUE(transient.precip.change)) {
 
     precip_meanf <- sapply(1:12,
-      function(m) seq(1, change.factor.precip.mean[m], length.out = max(year_ind)))
+      function(m) seq(1, change.factor.precip.mean[m] * 2, length.out = max(year_ind)))
     precip_varf  <- sapply(1:12,
-      function(m) seq(1, change.factor.precip.variance[m], length.out = max(year_ind)))
+      function(m) seq(1, change.factor.precip.variance[m] * 2, length.out = max(year_ind)))
 
-  } else {
+    # Step changes
+    } else {
     precip_meanf <- sapply(1:12,
-      function(m) rep(change.factor.precip.mean[m]/2, max(year_ind)))
+      function(m) rep(change.factor.precip.mean[m], max(year_ind)))
     precip_varf  <- sapply(1:12,
-      function(m) rep(change.factor.precip.variance[m]/2, max(year_ind)))
+      function(m) rep(change.factor.precip.variance[m], max(year_ind)))
   }
 
   ##############################################################################
 
   # Apply climate changes (per grid)
   for (x in 1:ngrids) {
 
-    # Perturb daily precipitation using quantile mapping
+    # Perturb daily precipitation (by quantile mapping)
     climate.data[[x]]$precip <- weathergenr::quantileMapping(
             value = climate.data[[x]]$precip,
             mon.ts = month_ind,
@@ -79,74 +90,77 @@ imposeClimateChanges <- function(
             var.change = precip_varf,
             fit.method = "mme")
 
-    # Perturb temp, temp_min, and temp_max by delta factors
+    # Perturb temp, temp_min, and temp_max (by additive delta factors)
     climate.data[[x]]$temp <- climate.data[[x]]$temp + tempf2
     climate.data[[x]]$temp_min <- climate.data[[x]]$temp_min + tempf2
     climate.data[[x]]$temp_max <- climate.data[[x]]$temp_max + tempf2
 
+    # Calculate adjusted PET based on temperature variables
     if(isTRUE(calculate.pet)) {
-      climate.data[[x]]$pet <- with(climate.data[[x]], hargreavesPet(
-          months = month_ind, temp = temp, tdiff = temp_max - temp_min,
-          lat = climate.grid$y[x]))
+      climate.data[[x]]$pet <- with(climate.data[[x]],
+                hargreavesPet(months = month_ind,
+                              temp = temp,
+                              tdiff = temp_max - temp_min,
+                              lat = climate.grid$y[x]))
     }
 
   }
 
-  # Replace possible infinite/NA values with zero
+  # Replace possible infinite/NA values with zeros
   climate.data <- lapply(1:length(climate.data), function(y)
-    do.call(tibble, lapply(climate.data[[y]], function(x) replace(x, is.infinite(x), 0))))
-
-
-  # # Set number of cores for parallel computing
-  # if(compute.parallel == TRUE) {
-  #
-  #   if(is.null(num.cores)) num.cores <- parallel::detectCores()-1
-  #   cl <- parallel::makeCluster(num.cores)
-  #   doParallel::registerDoParallel(cl)
-  #   `%d%` <- foreach::`%dopar%`
-  #
-  # } else {
-  #
-  #   `%d%` <- foreach::`%do%`
-  # }
-
-  ##############################################################################
-  ##############################################################################
-
-  # precip <- foreach::foreach(x=seq_len(ngrids)) %d% {
-  #
-  #     weathergenr::quantileMapping(
-  #           value = climate.data[[x]]$precip,
-  #           mon.ts = month_ind,
-  #           year.ts = year_ind,
-  #           mean.change = precip_meanf,
-  #           var.change = precip_varf,
-  #           fit.method = fit.method)
-  # }
-  #
-  # if(compute.parallel == TRUE) parallel::stopCluster(cl)
-  #
-  # for (x in 1:ngrids) {
-  #
-  #   # Perturb temp, temp_min, and temp_max by delta factors
-  #   climate.data[[x]]$precip <- precip[[x]]
-  #   climate.data[[x]]$temp   <- climate.data[[x]]$temp + tempf2
-  #   climate.data[[x]]$temp_min <- climate.data[[x]]$temp_min + tempf2
-  #   climate.data[[x]]$temp_max <- climate.data[[x]]$temp_max + tempf2
-  #
-  #   if(isTRUE(calculate.pet)) {
-  #     climate.data[[x]]$pet <- with(climate.data[[x]], hargreavesPet(
-  #         months = month_ind, temp = temp, tdiff = temp_max - temp_min,
-  #         lat = climate.grid$y[x]))
-  #   }
-  #
-  # }
-  #
-  # # Replace possible infinite/NA values with zero
-  # climate.data <- lapply(1:length(climate.data), function(y)
-  #       do.call(tibble::tibble, lapply(climate.data[[y]],
-  #         function(x) replace(x, is.infinite(x) | is.na(x), 0))))
+    do.call(tibble, lapply(climate.data[[y]],
+                           function(x) replace(x, is.infinite(x), 0))))
 
   return(climate.data)
 
 }
+
+# # Set number of cores for parallel computing
+# if(compute.parallel == TRUE) {
+#
+#   if(is.null(num.cores)) num.cores <- parallel::detectCores()-1
+#   cl <- parallel::makeCluster(num.cores)
+#   doParallel::registerDoParallel(cl)
+#   `%d%` <- foreach::`%dopar%`
+#
+# } else {
+#
+#   `%d%` <- foreach::`%do%`
+# }
+
+##############################################################################
+##############################################################################
+
+# precip <- foreach::foreach(x=seq_len(ngrids)) %d% {
+#
+#     weathergenr::quantileMapping(
+#           value = climate.data[[x]]$precip,
+#           mon.ts = month_ind,
+#           year.ts = year_ind,
+#           mean.change = precip_meanf,
+#           var.change = precip_varf,
+#           fit.method = fit.method)
+# }
+#
+# if(compute.parallel == TRUE) parallel::stopCluster(cl)
+#
+# for (x in 1:ngrids) {
+#
+#   # Perturb temp, temp_min, and temp_max by delta factors
+#   climate.data[[x]]$precip <- precip[[x]]
+#   climate.data[[x]]$temp   <- climate.data[[x]]$temp + tempf2
+#   climate.data[[x]]$temp_min <- climate.data[[x]]$temp_min + tempf2
+#   climate.data[[x]]$temp_max <- climate.data[[x]]$temp_max + tempf2
+#
+#   if(isTRUE(calculate.pet)) {
+#     climate.data[[x]]$pet <- with(climate.data[[x]], hargreavesPet(
+#         months = month_ind, temp = temp, tdiff = temp_max - temp_min,
+#         lat = climate.grid$y[x]))
+#   }
+#
+# }
+#
+# # Replace possible infinite/NA values with zero
+# climate.data <- lapply(1:length(climate.data), function(y)
+#       do.call(tibble::tibble, lapply(climate.data[[y]],
+#         function(x) replace(x, is.infinite(x) | is.na(x), 0))))