Specify type of soil density inputs

- close #209 (Integrate SOILWAT2's soil density descriptor) - see DrylandEcology/SOILWAT2#324 (SOILWAT2 submodule points to this devel branch for now) - new behavior: user specifies whether input soil density represents matric or bulk density; code then determines both matric and bulk density values from inputs - Class `swSite` gains new slot `"SoilDensityInputType"` and associated methods `swSite_SoilDensityInputType()` (#209). "SoilDensityInputType" encodes if soil density inputs represent matric soil or bulk soil values.
DrylandEcology · Aug 2, 2022 · 4179476 · 4179476
1 parent f8c1e44
commit 4179476
Show file tree

Hide file tree

Showing 20 changed files with 154 additions and 10 deletions.
diff --git a/.gitmodules b/.gitmodules
@@ -1,4 +1,4 @@
 [submodule "src/SOILWAT2"]
 	path = src/SOILWAT2
 	url = https://github.com/DrylandEcology/SOILWAT2
-	branch = master
+	branch = feature_soildensity_type
diff --git a/NAMESPACE b/NAMESPACE
@@ -211,6 +211,7 @@ exportMethods("swSite_ModelCoefficients<-")
 exportMethods("swSite_ModelFlags<-")
 exportMethods("swSite_SWClimits<-")
 exportMethods("swSite_SnowSimulationParams<-")
+exportMethods("swSite_SoilDensityInputType<-")
 exportMethods("swSite_SoilTemperatureConsts<-")
 exportMethods("swSite_SoilTemperatureFlag<-")
 exportMethods("swSite_TranspCoefficients<-")
@@ -307,6 +308,7 @@ exportMethods(swSite_ModelCoefficients)
 exportMethods(swSite_ModelFlags)
 exportMethods(swSite_SWClimits)
 exportMethods(swSite_SnowSimulationParams)
+exportMethods(swSite_SoilDensityInputType)
 exportMethods(swSite_SoilTemperatureConsts)
 exportMethods(swSite_SoilTemperatureFlag)
 exportMethods(swSite_TranspCoefficients)

diff --git a/NEWS.md b/NEWS.md
@@ -1,6 +1,9 @@
 # rSOILWAT2 devel
 * `r-lib` Github Actions updated to `v2`;
   separate workflows for `R-CMD-check` and `test-coverage` (#202).
+* Class `swSite` gains new slot `"SoilDensityInputType"` and associated
+  methods `swSite_SoilDensityInputType()` (#209).
+  This encodes if soil density inputs represent matric soil or bulk soil values.
 
 
 # rSOILWAT2 v5.3.1

diff --git a/R/A_swGenericMethods.R b/R/A_swGenericMethods.R
@@ -1483,6 +1483,16 @@ setGeneric(
   function(object) standardGeneric("swSite_SoilTemperatureConsts")
 )
 
+#' \code{swSite_SoilTemperatureFlag}
+#' @param object An object of class \code{\linkS4class{swSite}} or
+#'   \code{\linkS4class{swInputData}}.
+#' @seealso \code{\linkS4class{swSite}} and \code{\linkS4class{swInputData}}
+setGeneric(
+  "swSite_SoilDensityInputType",
+  function(object) standardGeneric("swSite_SoilDensityInputType")
+)
+
+
 #' \code{swSite_TranspirationRegions}
 #' @param object An object of class \code{\linkS4class{swSite}} or
 #'   \code{\linkS4class{swInputData}}.
@@ -1608,6 +1618,16 @@ setGeneric(
   function(object, value) standardGeneric("swSite_SoilTemperatureConsts<-")
 )
 
+#' \code{swSite_SoilDensityInputType<-}
+#' @param object An object of class \code{\linkS4class{swSite}} or
+#'   \code{\linkS4class{swInputData}}.
+#' @param value A value to assign to a specific slot of the \code{object}.
+#' @seealso \code{\linkS4class{swSite}} and \code{\linkS4class{swInputData}}
+setGeneric(
+  "swSite_SoilDensityInputType<-",
+  function(object, value) standardGeneric("swSite_SoilDensityInputType<-")
+)
+
 #' \code{swSite_TranspirationRegions<-}
 #' @param object An object of class \code{\linkS4class{swSite}} or
 #'   \code{\linkS4class{swInputData}}.

diff --git a/R/F_swSite.R b/R/F_swSite.R
@@ -60,6 +60,7 @@ setClass(
     IntrinsicSiteParams = "numeric",
     SoilTemperatureFlag = "logical",
     SoilTemperatureConstants = "numeric",
+    SoilDensityInputType = "integer",
     TranspirationRegions = "matrix"
   )
 )
@@ -128,6 +129,14 @@ setValidity("swSite", function(object) {
     msg <- "@SoilTemperatureConstants length != 10."
     val <- if (isTRUE(val)) msg else c(val, msg)
   }
+  if (typeof(object@SoilDensityInputType) != "integer") {
+    msg <- "@SoilDensityInputType is of integer type."
+    val <- if (isTRUE(val)) msg else c(val, msg)
+  }
+  if (length(object@SoilDensityInputType) != 1L) {
+    msg <- "@SoilDensityInputType length != 1."
+    val <- if (isTRUE(val)) msg else c(val, msg)
+  }
   if (typeof(object@TranspirationRegions) != "integer") {
     msg <- "@TranspirationRegions is of integer type."
     val <- if (isTRUE(val)) msg else c(val, msg)
@@ -239,6 +248,14 @@ setMethod("swSite_SoilTemperatureFlag", "swSite",
 setMethod("swSite_SoilTemperatureConsts", "swSite",
   function(object) slot(object, "SoilTemperatureConstants"))
 
+#' @rdname swSite-class
+#' @export
+setMethod(
+  "swSite_SoilDensityInputType",
+  "swSite",
+  function(object) slot(object, "SoilDensityInputType")
+)
+
 #' @rdname swSite-class
 #' @export
 setMethod("swSite_TranspirationRegions", "swSite",
@@ -343,6 +360,17 @@ setReplaceMethod("swSite_SoilTemperatureConsts", signature = "swSite",
     object
 })
 
+#' @rdname swSite-class
+#' @export
+setReplaceMethod(
+  "swSite_SoilDensityInputType",
+  signature = "swSite",
+  definition = function(object, value) {
+    object@SoilDensityInputType <- as.integer(value[1L])
+    validObject(object)
+    object
+})
+
 #' @rdname swSite-class
 #' @export
 setReplaceMethod("swSite_TranspirationRegions", signature = "swSite",

diff --git a/R/K_swContainer.R b/R/K_swContainer.R
@@ -1503,6 +1503,14 @@ setMethod(
   function(object) swSite_SoilTemperatureConsts(object@site)
 )
 
+#' @rdname swInputData-class
+#' @export
+setMethod(
+  "swSite_SoilDensityInputType",
+  signature = "swInputData",
+  function(object) swSite_SoilDensityInputType(object@site)
+)
+
 #' @rdname swInputData-class
 #' @export
 setMethod(
@@ -1633,6 +1641,17 @@ setReplaceMethod(
   }
 )
 
+#' @rdname swInputData-class
+#' @export
+setReplaceMethod(
+  "swSite_SoilDensityInputType",
+  signature = "swInputData",
+  function(object, value) {
+    swSite_SoilDensityInputType(object@site) <- value
+    object
+  }
+)
+
 #' @rdname swInputData-class
 #' @export
 setReplaceMethod(

diff --git a/data/sw_exampleData.rda b/data/sw_exampleData.rda
diff --git a/inst/extdata/example1/Input/siteparam.in b/inst/extdata/example1/Input/siteparam.in
@@ -83,6 +83,10 @@ NAN			# aspect = surface azimuth angle (degrees): S=0, E=-90, N=±180, W=90;
 # Name of CO2 scenario: see input file `carbon.in`
 RCP85
 
+# --- Soil characterization ---
+# Are inputs of density representing bulk soil (type 1) or the matric component (type 0)?
+0
+
 #---- Transpiration regions
 # ndx  : 1=shallow, 2=medium, 3=deep, 4=very deep
 # layer: deepest soil layer number of the region.

diff --git a/man/swInputData-class.Rd b/man/swInputData-class.Rd
diff --git a/man/swSite-class.Rd b/man/swSite-class.Rd
diff --git a/man/swSite_SoilDensityInputType-set.Rd b/man/swSite_SoilDensityInputType-set.Rd
diff --git a/man/swSite_SoilDensityInputType.Rd b/man/swSite_SoilDensityInputType.Rd
diff --git a/src/SOILWAT2 b/src/SOILWAT2
diff --git a/src/rSW_Site.c b/src/rSW_Site.c
@@ -49,7 +49,9 @@ static char *cSW_SIT[] = {
   "SWClimits", "ModelFlags", "ModelCoefficients",
   "SnowSimulationParameters", "DrainageCoefficient", "EvaporationCoefficients",
   "TranspirationCoefficients", "IntrinsicSiteParams", "SoilTemperatureFlag",
-  "SoilTemperatureConstants", "TranspirationRegions"
+  "SoilTemperatureConstants",
+  "SoilDensityInputType",
+  "TranspirationRegions"
 };
 
 static char *cLayers[] = {
@@ -76,7 +78,7 @@ SEXP onGet_SW_LYR() {
 	p_Layers = REAL(Layers);
 	for (i = 0; i < (v->n_layers); i++) {
 		p_Layers[i + (v->n_layers) * 0] = dmax = v->lyr[i]->width + dmax;
-		p_Layers[i + (v->n_layers) * 1] = v->lyr[i]->soilMatric_density;
+		p_Layers[i + (v->n_layers) * 1] = v->lyr[i]->soilDensityInput;
 		p_Layers[i + (v->n_layers) * 2] = v->lyr[i]->fractionVolBulk_gravel;
 		p_Layers[i + (v->n_layers) * 3] = v->lyr[i]->evap_coeff;
 		p_Layers[i + (v->n_layers) * 4] = v->lyr[i]->transp_coeff[SW_GRASS];
@@ -107,7 +109,7 @@ void onSet_SW_LYR(SEXP SW_SOILS) {
 	SW_SITE *v = &SW_Site;
 	LyrIndex lyrno;
 	int i, j, k, columns;
-	RealF dmin = 0.0, dmax, evco, trco_veg[NVEGTYPES], psand, pclay, matricd, imperm, soiltemp, f_gravel;
+	RealF dmin = 0.0, dmax, evco, trco_veg[NVEGTYPES], psand, pclay, soildensity, imperm, soiltemp, f_gravel;
 	RealD *p_Layers;
 	SEXP SW_LYR;
 
@@ -135,7 +137,7 @@ void onSet_SW_LYR(SEXP SW_SOILS) {
 		lyrno = _newlayer();
 
 		dmax = p_Layers[i + j * 0];
-		matricd = p_Layers[i + j * 1];
+		soildensity = p_Layers[i + j * 1];
 		f_gravel = p_Layers[i + j * 2];
 		evco = p_Layers[i + j * 3];
 		trco_veg[SW_GRASS] = p_Layers[i + j * 4];
@@ -149,7 +151,7 @@ void onSet_SW_LYR(SEXP SW_SOILS) {
 
 		v->lyr[lyrno]->width = dmax - dmin;
 		dmin = dmax;
-		v->lyr[lyrno]->soilMatric_density = matricd;
+		v->lyr[lyrno]->soilDensityInput = soildensity;
 		v->lyr[lyrno]->fractionVolBulk_gravel = f_gravel;
 		v->lyr[lyrno]->evap_coeff = evco;
 		ForEachVegType(k) {
@@ -207,6 +209,8 @@ SEXP onGet_SW_SIT() {
 	char *cSoilTempValues[] = { "BiomassLimiter_g/m^2", "T1constant_a", "T1constant_b", "T1constant_c", "cs_constant_SoilThermCondct", "cs_constant", "sh_constant_SpecificHeatCapacity",
 			"ConstMeanAirTemp", "deltaX_Param", "MaxDepth" };
 
+	SEXP SoilDensityInputType;
+
 	SEXP TranspirationRegions, TranspirationRegions_names, TranspirationRegions_names_y;
 	char *cTranspirationRegions[] = { "ndx", "layer" };
 	int *p_transp; // ideally `LyrIndex` so that same type as `_TranspRgnBounds`, but R API INTEGER() is signed
@@ -310,6 +314,8 @@ SEXP onGet_SW_SIT() {
 		SET_STRING_ELT(SoilTemperatureConstants_names, i, mkChar(cSoilTempValues[i]));
 	setAttrib(SoilTemperatureConstants, R_NamesSymbol, SoilTemperatureConstants_names);
 
+	PROTECT(SoilDensityInputType = ScalarInteger(v->type_soilDensityInput));
+
 	PROTECT(TranspirationRegions = allocMatrix(INTSXP,(v->n_transp_rgn),2));
 	p_transp = INTEGER(TranspirationRegions);
 	for (i = 0; i < (v->n_transp_rgn); i++) {
@@ -333,9 +339,10 @@ SEXP onGet_SW_SIT() {
 	SET_SLOT(SW_SIT, install(cSW_SIT[7]), IntrinsicSiteParams);
 	SET_SLOT(SW_SIT, install(cSW_SIT[8]), SoilTemperatureConstants_use);
 	SET_SLOT(SW_SIT, install(cSW_SIT[9]), SoilTemperatureConstants);
-	SET_SLOT(SW_SIT, install(cSW_SIT[10]), TranspirationRegions);
+	SET_SLOT(SW_SIT, install(cSW_SIT[10]), SoilDensityInputType);
+	SET_SLOT(SW_SIT, install(cSW_SIT[11]), TranspirationRegions);
 
-	UNPROTECT(24);
+	UNPROTECT(25);
 	return SW_SIT;
 }
 
@@ -353,7 +360,9 @@ void onSet_SW_SIT(SEXP SW_SIT) {
 	SEXP IntrinsicSiteParams;
 	SEXP SoilTemperatureConstants_use;
 	SEXP SoilTemperatureConstants;
+	SEXP SoilDensityInputType;
 	SEXP TranspirationRegions;
+
 	int *p_transp; // ideally `LyrIndex` so that same type as `_TranspRgnBounds`, but R API INTEGER() is signed
 
   #ifdef RSWDEBUG
@@ -459,6 +468,12 @@ void onSet_SW_SIT(SEXP SW_SIT) {
 	if (debug) swprintf(" > 'soiltemp-constants'");
 	#endif
 
+	PROTECT(SoilDensityInputType = GET_SLOT(SW_SIT, install("SoilDensityInputType")));
+	v->type_soilDensityInput = INTEGER(SoilDensityInputType)[0];
+	#ifdef RSWDEBUG
+	if (debug) swprintf(" > 'density-type'");
+	#endif
+
 	PROTECT(TranspirationRegions = GET_SLOT(SW_SIT, install("TranspirationRegions")));
 	p_transp = INTEGER(TranspirationRegions);
 	v->n_transp_rgn = nrows(TranspirationRegions);
@@ -488,5 +503,5 @@ void onSet_SW_SIT(SEXP SW_SIT) {
 	if (debug) swprintf(" ... done. \n");
 	#endif
 
-	UNPROTECT(11);
+	UNPROTECT(12);
 }
diff --git a/tests/test_data/Ex1_input.rds b/tests/test_data/Ex1_input.rds
diff --git a/tests/test_data/Ex2_input.rds b/tests/test_data/Ex2_input.rds
diff --git a/tests/test_data/Ex3_input.rds b/tests/test_data/Ex3_input.rds
diff --git a/tests/test_data/Ex4_input.rds b/tests/test_data/Ex4_input.rds
diff --git a/tests/test_data/Ex5_input.rds b/tests/test_data/Ex5_input.rds
diff --git a/vignettes/rSOILWAT2_demo.Rmd b/vignettes/rSOILWAT2_demo.Rmd
@@ -463,6 +463,9 @@ You may organize weather data in a variety of ways:
 
 ## Assign new soil data to rSOILWAT2 input object
 ```{r, soils_final1}
+  # Specify that our soil density values are bulk density
+  rSOILWAT2::swSite_SoilDensityInputType(sw_in) <- 1L
+
   # This fails because rooting profile values are still missing
   try(rSOILWAT2::swSoils_Layers(sw_in) <- data.matrix(soil_new))
 ```