Merge pull request #115 from Burke-Lauenroth-Lab/new_CO2

new_CO2: simulation experiments that can use the CO2-concentration effects on biomass and water-use efficiency offered by SOILWAT2

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: rSFSW2
 Title: Simulation Framework for SOILWAT2
-Version: 2.4.1
-Date: 2017-12-12
+Version: 2.5.1
+Date: 2017-12-15
 Authors@R: c(person("Daniel", "Schlaepfer", email = "daniel.schlaepfer@yale.edu", role = c("aut", "cre")),
       person("Caitlin", "Andrews", role = "ctb"),
       person("Zach", "Kramer", role = "ctb"),
@@ -11,7 +11,7 @@ Description: Setting up, carrying out, and analyzing ecosystem water balance
 Depends:
     R (>= 3.3.0)
 Imports:
-    rSOILWAT2 (>= 1.6.0),
+    rSOILWAT2 (>= 1.7.0),
     RSQLite (>= 2.0),
     DBI (>= 0.7),
     Rcpp (>= 0.12.12),

NAMESPACE

@@ -30,6 +30,7 @@ export(calc_ExtendSoilDatafileToRequestedSoilLayers)
 export(calc_SiteClimate)
 export(check_cltool)
 export(check_lock_content)
+export(check_monotonic_increase)
 export(check_outputDB_completeness)
 export(check_rSFSW2_project_input_data)
 export(check_weatherDB)

R/IO_datafiles.R

@@ -365,7 +365,7 @@ complete_with_defaultpaths <- function(project_paths, fnames_in) {
   }
 
   # full names of files located in 'dir_in_treat'
-  ftemp <- c("LookupClimatePPTScenarios", "LookupClimateTempScenarios",
+  ftemp <- c("LookupCarbonScenarios", "LookupClimatePPTScenarios", "LookupClimateTempScenarios",
     "LookupShiftedPPTScenarios", "LookupEvapCoeffFromTable", "LookupTranspCoeffFromTable",
     "LookupTranspRegionsFromTable", "LookupSnowDensityFromTable",
     "LookupVegetationComposition")
@@ -465,6 +465,8 @@ process_inputs <- function(project_paths, fnames_in, use_preprocin = TRUE, verbo
     sw_input_soillayers <- tryCatch(SFSW2_read_csv(fnames_in[["fslayers"]],
       nrowsClasses = nrowsClasses), error = print)
     sw_input_soillayers <- fix_rowlabels(sw_input_soillayers, SWRunInformation)
+    sw_input_soillayers[, -(1:2)] <- check_monotonic_increase(data.matrix(sw_input_soillayers[, -(1:2)]),
+      strictly = TRUE, fail = TRUE, na.rm = TRUE)
 
     temp <- tryCatch(SFSW2_read_inputfile(fnames_in[["ftreatDesign"]],
       nrowsClasses = nrowsClasses), error = print)
@@ -554,13 +556,16 @@ process_inputs <- function(project_paths, fnames_in, use_preprocin = TRUE, verbo
     tr_weather <- load_Rsw_treatment_templates(project_paths, create_treatments, "weathersetupin", "swWeather")
     tr_cloud <- load_Rsw_treatment_templates(project_paths, create_treatments, "cloudin", "swCloud")
 
-    tr_input_climPPT <- tr_input_climTemp <- tr_input_shiftedPPT <- list()
+    tr_input_CarbonScenario <- tr_input_climPPT <- tr_input_climTemp <- tr_input_shiftedPPT <- list()
     tr_input_EvapCoeff <- tr_input_TranspCoeff_Code <- tr_input_TranspCoeff <- list()
     tr_input_TranspRegions <- tr_input_SnowD <- tr_VegetationComposition <- list()
 
     if (any(create_treatments == "LookupClimatePPTScenarios"))
       tr_input_climPPT <- SFSW2_read_csv(fnames_in[["LookupClimatePPTScenarios"]])
 
+    if (any(create_treatments == "LookupCarbonScenarios"))
+      tr_input_CarbonScenario <- SFSW2_read_csv(fnames_in[["LookupCarbonScenarios"]])
+
     if (any(create_treatments == "LookupClimateTempScenarios"))
       tr_input_climTemp <- SFSW2_read_csv(fnames_in[["LookupClimateTempScenarios"]])
 
@@ -635,7 +640,7 @@ process_inputs <- function(project_paths, fnames_in, use_preprocin = TRUE, verbo
       sw_input_climscen_use = sw_input_climscen_use, sw_input_climscen = sw_input_climscen,
       sw_input_climscen_values_use = sw_input_climscen_values_use, sw_input_climscen_values = sw_input_climscen_values,
       tr_files = tr_files, tr_prod = tr_prod, tr_site = tr_site, tr_soil = tr_soil,
-      tr_weather = tr_weather, tr_cloud = tr_cloud, tr_input_climPPT = tr_input_climPPT,
+      tr_weather = tr_weather, tr_cloud = tr_cloud, tr_input_CarbonScenario = tr_input_CarbonScenario, tr_input_climPPT = tr_input_climPPT,
       tr_input_climTemp = tr_input_climTemp, tr_input_shiftedPPT = tr_input_shiftedPPT,
       tr_input_EvapCoeff = tr_input_EvapCoeff, tr_input_TranspCoeff_Code = tr_input_TranspCoeff_Code,
       tr_input_TranspCoeff = tr_input_TranspCoeff, tr_input_TranspRegions = tr_input_TranspRegions,

R/Mathematical_Functions.R

@@ -304,3 +304,53 @@ do_compare <- function(x1, x2) {
     }
   }
 }
+
+
+#' Check that values in matrix-like object are (strictly) monotonically increasing/decreasing
+#'
+#' @param x A numeric matrix like object.
+#' @param MARGIN An integer value giving the subscripts over which the monotonicity will
+#'   be checked; 1 indicates rows, 2 indicates columns.
+#' @param increase A logical value. If \code{TRUE}, check monotic increase; if
+#'   \code{FALSE}, check monotonic decrease.
+#' @param strictly A logical value. If \code{TRUE}, check for a strict monotic pattern.
+#' @param fail A logical value. If \code{TRUE}, throw error if monotic check fails.
+#' @param replacement A value that replaces non-(strictly) monotically increasing/decreasing
+#'   values if \code{fail} is \code{FALSE}.
+#' @param na.rm A logical value. If \code{TRUE}, then ignore \code{NA}s; if \code{FALSE},
+#'   then fail if \code{strictly} or replace with \code{replacement}.
+#' @return The updated \code{x}.
+#' @export
+check_monotonic_increase <- function(x, MARGIN = 1, increase = TRUE, strictly = FALSE,
+  fail = FALSE, replacement = NA, na.rm = FALSE) {
+
+  stopifnot(MARGIN %in% c(1, 2), length(dim(x)) == 2)
+
+  x <- as.matrix(x)
+  if (MARGIN == 2) {
+    x <- t(x)
+  }
+
+  mfun <- if (increase) {
+      if (strictly) '>' else '>='
+    } else {
+      if (strictly) '<' else '=<'
+    }
+
+  ord <- !match.fun(mfun)(x[, -1, drop = FALSE], x[, -ncol(x), drop = FALSE])
+  has_na <- is.na(x)
+
+  if (any(ord, na.rm = TRUE) || (has_na && !na.rm && strictly)) {
+    if (fail) {
+      stop(paste0("'check_monotonic_increase': data are not ", if (strictly) "strictly ",
+        "monotically ", if (increase) "increasing " else "decreasing ",
+        if (MARGIN == 1) "in rows." else "in columns."))
+
+    } else {
+      x[, -1][is.na(ord) | ord] <- replacement
+      x[is.na(x[, 1]), 1] <- replacement
+    }
+  }
+
+  if (MARGIN == 1) x else t(x)
+}

R/OutputDatabase.R

@@ -1550,6 +1550,11 @@ dbOutput_create_Design <- function(con_dbOut, SFSW2_prj_meta, SFSW2_prj_inputs)
     if (aon$input_VegetationBiomassMonthly) {
       temp <- c(temp, paste0(c(rep("Grass", 36), rep("Shrub", 36), rep("Tree", 36), rep("Forb", 36)), "_", c(rep("Litter", 12), rep("TotalBiomass", 12), rep("LiveBiomass", 12)), "_m", SFSW2_glovars[["st_mo"]], "_gPERm2"))
     }
+  #2b
+    if (aon$input_VegetationBiomassTrends) {
+      temp <- c(temp, paste0(rep(c("Grass", "Shrub", "Tree", "Forb", "Total"), 2), "_",
+        rep(c("Total", "Live"), each = 5), "Biomass_gPERm2_mean"))
+    }
   #3.
     if (aon$input_VegetationPeak) {
       temp <- c(temp, paste0("SWinput.PeakLiveBiomass_", c("month_mean", "months_duration")))
@@ -1593,6 +1598,11 @@ dbOutput_create_Design <- function(con_dbOut, SFSW2_prj_meta, SFSW2_prj_inputs)
     if (aon$input_ClimatePerturbations) {
       temp <- c(temp, paste0(rep(paste0("SWinput.ClimatePerturbations.", c("PrcpMultiplier.m", "TmaxAddand.m", "TminAddand.m")), each = 12), SFSW2_glovars[["st_mo"]], rep(c("_none", "_C", "_C"), each = 12), "_const"))
     }
+  #6b
+    if (aon$input_CO2Effects) {
+      temp <- c(temp, paste0(rep(c("Grass", "Shrub", "Tree", "Forb"), 2), "_",
+        rep(c("Biomass", "WUE"), each = 4), "_CO2multiplier_fraction_mean"))
+    }
 
     ##############################################################---Aggregation: Climate and weather---##############################################################
 
@@ -2131,7 +2141,8 @@ dbOutput_create_DailyAggregationTable <- function(con_dbOut, req_aggs) {
         EvaporationTotal = 1, VWCbulk = 2, VWCmatric = 2, SWCbulk = 2, SWPmatric = 2,
         SWAbulk = 2, Snowpack = 1, Rain = 1, Snowfall = 1, Snowmelt = 1, SnowLoss = 1,
         Infiltration = 1, DeepDrainage = 1, PET = 1, TotalPrecipitation = 1,
-        TemperatureMin = 1, TemperatureMax = 1, SoilTemperature = 2, Runoff = 1)
+        TemperatureMin = 1, TemperatureMax = 1, SoilTemperature = 2, Runoff = 1,
+        Runon = 1)
       tableName <- paste0("aggregation_doy_", req_aggs[["tag"]][doi])
 
       if (agg.analysis == 1) {
@@ -2204,7 +2215,8 @@ dbOutput_create_EnsembleTables <- function(con_dbOut, dbOutput, prj_todos, sim_s
             EvaporationTotal = 1, VWCbulk = 2, VWCmatric = 2, SWCbulk = 2, SWPmatric = 2,
             SWAbulk = 2, Snowpack = 1, Rain = 1, Snowfall = 1, Snowmelt = 1, SnowLoss = 1,
             Infiltration = 1, DeepDrainage = 1, PET = 1, TotalPrecipitation = 1,
-            TemperatureMin = 1, TemperatureMax = 1, SoilTemperature = 2, Runoff = 1)
+            TemperatureMin = 1, TemperatureMax = 1, SoilTemperature = 2, Runoff = 1,
+            Runon = 1)
 
           if (agg.analysis == 1) {
             sql1 <- paste0("CREATE TABLE \"", EnsembleFamilyLevelTables[1], "\" (",

R/Simulation_Project.R

@@ -850,7 +850,7 @@ simulate_SOILWAT2_experiment <- function(actions, SFSW2_prj_meta, SFSW2_prj_inpu
   if (SFSW2_prj_meta[["prj_todos"]][["use_SOILWAT2"]] &&
     SFSW2_prj_meta[["sim_size"]][["runsN_todo"]] > 0) {
 
-    swof <- sw_out_flags()
+    swof <- rSOILWAT2::sw_out_flags()
     swDataFromFiles <- read_SOILWAT2_FileDefaults(SFSW2_prj_meta[["project_paths"]][["dir_in_sw"]])
     args_do_OneSite <- gather_args_do_OneSite(SFSW2_prj_meta, SFSW2_prj_inputs)