\pagebreak
options(tinytex.verbose = TRUE)
v <- try(utils::packageVersion("rSOILWAT2"), silent = TRUE) if (inherits(v, "try-error") || v < "4.0.3") { message( "There is no or an outdated rSOILWAT2 version present. ", "Please update to a more recent version and build vignettes: see", "https://github.com/DrylandEcology/rSOILWAT2#installation" ) }
# Package documentation help(package = "rSOILWAT2") # Documentation of a specific function, e.g., ?rSOILWAT2::sw_exec # View code of a specific function, e.g., rSOILWAT2::sw_exec
?Extract #--- Subset a list x <- list(ID = 1, Lyr_1 = 1:10, Lyr_2 = 11:30) x[[1]] # extract element by position x[["Lyr_1"]] # extract element by name x[1:2] # subset by position #--- Subset a data.frame x <- as.data.frame(x) # extract column by column number x[[1]] x[, 1] # subset by column numbers x[1:2] x[, 1:2] # extract column by column name x[["Lyr_1"]] x[, "Lyr_1"] x$Lyr_1 #--- Subset a matrix x <- as.matrix(x) x[, 1] # extract element by position x[, "Lyr_1"] # extract element by name x[, 1:2] # subset by position
# Load data for an example site: see ?sw_exampleData sw_in0 <- rSOILWAT2::sw_exampleData # Run SOILWAT2 for the example site sw_out0 <- rSOILWAT2::sw_exec(inputData = sw_in0)
class(sw_out0) slotNames(sw_out0)
clim0 <- rSOILWAT2::calc_SiteClimate( weatherList = rSOILWAT2::get_WeatherHistory(sw_in0), do_C4vars = TRUE )
x <- slot(slot(sw_out0, "VWCBULK"), "Day") ids <- grep("Lyr", colnames(x), fixed = TRUE) vwc <- x[, ids] matplot( vwc, ylim = c(0, max(vwc)), xlab = "Time", ylab = "Volumetric water content (bulk; m3 / m3)", type = "l" ) legend("bottom", ncol = 4, legend = colnames(vwc), lty = 1, col = seq_along(ids) ) x <- slot(slot(sw_out0, "TRANSP"), "Month") ids <- grep("transp_total", colnames(x), fixed = TRUE) matplot( x[1:24, ids], type = "l", xlab = "Time (Months)", ylab = "Transpiration (cm)" ) x <- slot(slot(sw_out0, "EVAPSOIL"), "Year") ids <- grep("Lyr", colnames(x), fixed = TRUE) es <- rowSums(x[, ids]) plot( x = x[, "Year"], y = es, ylim = c(0, max(es)), xlab = "Time", ylab = "Bare-soil evaporation (cm)", type = "l" )
path_demo <- system.file("extdata", "example1", package = "rSOILWAT2")
# Task: Use function `sw_inputDataFromFiles` to read inputs from location # `path_demo` (e.g., this could be inputs for your site) # Read data from input files sw_in1 <- rSOILWAT2::sw_inputDataFromFiles(dir = path_demo) # Then run simulation sw_out1 <- rSOILWAT2::sw_exec(inputData = sw_in1)
## Quick option # Note: This approach is discouraged because it is very easy to miss to # set all relevant parameters and variables to site-specific values sw_in__bad <- rSOILWAT2::sw_exampleData ## Don't do this! ## Preferred option # All relevant site-specific parameters and variables are set to NA sw_in <- rSOILWAT2::swInputData()
# Set simulation period, e.g., to 1980-2019 rSOILWAT2::swYears_StartYear(sw_in) <- 0 rSOILWAT2::swYears_EndYear(sw_in) <- 2019 rSOILWAT2::swYears_StartYear(sw_in) <- 1980 # Specify geographic location of site rSOILWAT2::swSite_IntrinsicSiteParams(sw_in) <- c( Longitude = -105.58, Latitude = 39.59, Altitude = 1000, Slope = 0, Aspect = NA )
You may organize weather data in a variety of ways:
- rSOILWAT2
weather database (see vignette rSOILWAT2_WeatherDatabase
)
- text files as for SOILWAT2
- directly access external online data sets
wdata <- rSOILWAT2::getWeatherData_folders( LookupWeatherFolder = file.path(path_demo, "Input"), weatherDirName = "data_weather" )
DayMet
has_daymet <- FALSE if ( requireNamespace("daymetr") && requireNamespace("curl") && curl::has_internet() ) { mm_dm <- try( rSOILWAT2::sw_meteo_obtain_DayMet( x = c(longitude = -105.6833, latitude = 41.3167), # Laramie WY start_year = rSOILWAT2::swYears_StartYear(sw_in), end_year = rSOILWAT2::swYears_EndYear(sw_in) ) ) if (!inherits(mm_dm, "try-error")) { # Impute values for added leap days # Use weather generator for available variables, otherwise interpolate wdata <- rSOILWAT2::dbW_imputeWeather( weatherData = mm_dm[["weatherDF"]], use_wg = TRUE, method_after_wg = "interp", seed = 123L ) # Check that weather data is well-formed stopifnot(rSOILWAT2::dbW_check_weatherData(wdata)) # Set use flags sw_in@weather@desc_rsds <- mm_dm[["desc_rsds"]] sw_in@weather@use_cloudCoverMonthly <- mm_dm[["use_cloudCoverMonthly"]] sw_in@weather@use_windSpeedMonthly <- mm_dm[["use_windSpeedMonthly"]] sw_in@weather@use_humidityMonthly <- mm_dm[["use_humidityMonthly"]] sw_in@weather@dailyInputFlags <- mm_dm[["dailyInputFlags"]] has_daymet <- TRUE } } if (!has_daymet) { # We don't have live internet and # were not able to obtain weather data for requested years # --> instead, use data we have locally and adjust years years <- sapply(wdata, slot, name = "year") rSOILWAT2::swYears_StartYear(sw_in) <- 0 rSOILWAT2::swYears_EndYear(sw_in) <- max(years) rSOILWAT2::swYears_StartYear(sw_in) <- min(years) # Set use flags sw_in@weather@desc_rsds <- 0L sw_in@weather@use_cloudCoverMonthly <- TRUE sw_in@weather@use_windSpeedMonthly <- TRUE sw_in@weather@use_humidityMonthly <- TRUE sw_in@weather@dailyInputFlags <- c(rep(TRUE, 3L), rep(FALSE, 11L)) }
This may not be the case if there was no internet connection to download
data from DayMet
yrs <- rSOILWAT2::get_years_from_weatherData(wdata) if (rSOILWAT2::swYears_EndYear(sw_in) > max(yrs)) { message( "Insufficient weather data to end simulations in year ", rSOILWAT2::swYears_EndYear(sw_in), ":\nAdjust end of simulation period to year ", max(yrs), "." ) rSOILWAT2::swYears_EndYear(sw_in) <- max(yrs) } if (rSOILWAT2::swYears_StartYear(sw_in) < min(yrs)) { message( "Insufficient weather data to start simulations in year ", rSOILWAT2::swYears_StartYear(sw_in), ":\nadjust start of simulation period to year ", min(yrs), "." ) rSOILWAT2::swYears_StartYear(sw_in) <- min(yrs) }
clim <- rSOILWAT2::calc_SiteClimate(weatherList = wdata, do_C4vars = TRUE)
# Name of the CO2 concentration scenario co2_nametag <- "RCP85" # Obtain yearly values from look-up table co2_data <- rSOILWAT2::lookup_annual_CO2a( start = rSOILWAT2::swYears_StartYear(sw_in), end = rSOILWAT2::swYears_EndYear(sw_in), name_co2 = co2_nametag ) # Assign CO2 values to rSOILWAT2 input object rSOILWAT2::swCarbon_Scenario(sw_in) <- co2_nametag rSOILWAT2::swCarbon_CO2ppm(sw_in) <- co2_data
# Extract from external data set, or # Use fixed values describing an experiment treatment (as here) rH <- rep(50, 12) # relative humidity [%] ws <- rep(2, 12) # wind speed [m/s] sc <- rep(30, 12) # sky cover [%] # Assign monthly climate normals to rSOILWAT2 input object rSOILWAT2::swCloud_Humidity(sw_in) <- rH rSOILWAT2::swCloud_WindSpeed(sw_in) <- ws rSOILWAT2::swCloud_SkyCover(sw_in) <- sc
soils_fixed <- data.frame( depth = c(5, 10, 20, 50), bulkd = 1.3, gravel = 0.1, evco = NA, trco_grass = NA, trco_shrub = NA, trco_tree = NA, trco_forb = NA, sand = 0.65, clay = 0.05, impermeability = NA, soil_temp = NA ) soil_new <- data.frame(rSOILWAT2::swSoils_Layers(sw_in)[0, ]) soil_new[seq_len(nrow(soils_fixed)), ] <- soils_fixed
NRCS SSURGO/STATSGO via SDA
if ( requireNamespace("rSW2exter") && requireNamespace("soilDB") && requireNamespace("curl") && curl::has_internet() ) { tmp <- rSOILWAT2::swSite_IntrinsicSiteParams(sw_in) soils_sda <- try(rSW2exter::extract_soils_NRCS_SDA( x = matrix(tmp[c("Longitude", "Latitude")], nrow = 1), method = "SSURGO_then_STATSGO", remove_organic_horizons = "at_surface", replace_missing_fragvol_with_zero = "at_surface", estimate_missing_bulkdensity = TRUE, restrict_by_ec_or_ph = FALSE, impute = TRUE, progress_bar = FALSE, verbose = FALSE )) if (inherits(soils_sda, "try-error")) { warning("Vignette was not able to extract soils from NRCS SDA!") } else { # Add homogenized soil layers tmp <- grep( "depth_L", colnames(soils_sda[["table_depths"]]), fixed = TRUE ) soils_sda2 <- rSW2data::update_soil_profile( soil_layers = soils_sda[["table_depths"]][1, tmp, drop = FALSE], requested_soil_layers = c( 5, seq(10, 100, by = 10), seq(120, 200, by = 20) ), soil_data = soils_sda[["table_texture"]], variables = c("dbovendry_L", "sandtotal_L", "claytotal_L", "fragvol_L"), vars_exhaust = NULL, keep_prev_soildepth = TRUE, keep_prev_soillayers = FALSE ) # Copy wide-data to rSOILWAT2-style soil data table soil_new <- data.frame(rSOILWAT2::swSoils_Layers(sw_in)[0, ]) ids <- seq_len(ncol(soils_sda2[["soil_layers"]])) # Soil layer depths tmp <- grep( "depth_L", colnames(soils_sda2[["soil_layers"]]), fixed = TRUE ) soil_new[ids, "depth_cm"] <- soils_sda2[["soil_layers"]][1, tmp] # Soil density, texture, and gravel cns_sd <- colnames(soils_sda2[["soil_data"]]) tmp <- grep("dbovendry_L", cns_sd, fixed = TRUE) soil_new[ids, "bulkDensity_g.cm.3"] <- soils_sda2[["soil_data"]][1, tmp] tmp <- grep("sandtotal_L", cns_sd, fixed = TRUE) soil_new[ids, "sand_frac"] <- soils_sda2[["soil_data"]][1, tmp] tmp <- grep("claytotal_L", cns_sd, fixed = TRUE) soil_new[ids, "clay_frac"] <- soils_sda2[["soil_data"]][1, tmp] tmp <- grep("fragvol_L", cns_sd, fixed = TRUE) soil_new[ids, "gravel_content"] <- soils_sda2[["soil_data"]][1, tmp] } }
# Set impermeability to zero soil_new[, "impermeability_frac"] <- 0
if (requireNamespace("rSW2data")) { soil_new[, "EvapBareSoil_frac"] <- rSW2data::calc_BareSoilEvapCoefs( layers_depth = soil_new[, "depth_cm"], sand = soil_new[, "sand_frac"], clay = soil_new[, "clay_frac"] )[1, ] }
soil_new[, "soilTemp_c"] <- rSW2data::init_soiltemperature( layers_depth = soil_new[, "depth_cm"], # Estimated soil surface temperature on Jan 1 # (potentially underneath snow pack) Tsoil_upper = max(-1, mean(clim[["meanMonthlyTempC"]][c(1, 12)])), # Constant soil temperature (Celsius) at the lower boundary (max depth) # approximated by mean annual air temperature Tsoil_const = mean(clim[["meanMonthlyTempC"]]), depth_Tsoil_const = 990 )
# NRCS SSURGO `dbovendry` represents the soil density of the matric component rSOILWAT2::swSite_SoilDensityInputType(sw_in) <- 0L # This fails because rooting profile values are still missing try(rSOILWAT2::swSoils_Layers(sw_in) <- soil_new)
# Note: `estimate_PotNatVeg_composition()` does not estimate fractional cover # of trees, bare-ground, etc. and sets those to zero or any other user # defined value veg_cover <- rSOILWAT2::estimate_PotNatVeg_composition( MAP_mm = 10 * clim[["MAP_cm"]], MAT_C = clim[["MAT_C"]], mean_monthly_ppt_mm = 10 * clim[["meanMonthlyPPTcm"]], mean_monthly_Temp_C = clim[["meanMonthlyTempC"]], dailyC4vars = clim[["dailyC4vars"]] ) ids <- sapply( X = names(rSOILWAT2::swProd_Composition(sw_in)), FUN = function(x) { grep( pattern = substr(x, 1, 4), x = names(veg_cover[["Rel_Abundance_L1"]]), ignore.case = TRUE ) } ) # Assign fractional cover values to rSOILWAT2 input object rSOILWAT2::swProd_Composition(sw_in) <- veg_cover[["Rel_Abundance_L1"]][ids]
# Reference biomass values from Bradford et al. 2014 are used # Mean monthly reference temperature corresponding to default phenology values # for the median across 898 big sagebrush sites are used veg_biom <- rSOILWAT2::estimate_PotNatVeg_biomass( target_temp = clim[["meanMonthlyTempC"]], target_MAP_mm = 10 * clim[["MAP_cm"]], do_adjust_phenology = TRUE, do_adjust_biomass = TRUE, fgrass_c3c4ann = veg_cover[["Grasses"]] ) # Assign monthly biomass values to rSOILWAT2 input object # Note: monthly biomass values of forbs, trees, etc. need to be estimated v1 <- c("Litter", "Biomass", "Perc.Live") v2 <- c("Litter", "Biomass", "Live_pct") rSOILWAT2::swProd_MonProd_grass(sw_in)[, v2] <- veg_biom[["grass"]][, v1] rSOILWAT2::swProd_MonProd_shrub(sw_in)[, v2] <- veg_biom[["shrub"]][, v1]
# Select rooting profile types # Set those to "FILL" where cover == 0 (because of transpiration regions) trco_type_by_veg <- list( grass_C3 = if (veg_cover[["Rel_Abundance_L0"]][["Grasses_C3"]] > 0) { "SchenkJackson2003_PCdry_grasses" } else { "FILL" }, grass_C4 = if (veg_cover[["Rel_Abundance_L0"]][["Grasses_C4"]] > 0) { "SchenkJackson2003_PCdry_grasses" } else { "FILL" }, grass_annuals = if ( veg_cover[["Rel_Abundance_L0"]][["Grasses_Annuals"]] > 0 ) { "Jacksonetal1996_crops" } else { "FILL" }, shrub = if (veg_cover[["Rel_Abundance_L0"]][["Shrubs"]] > 0) { "SchenkJackson2003_PCdry_shrubs" } else { "FILL" }, forb = if (veg_cover[["Rel_Abundance_L0"]][["Forbs"]] > 0) { "SchenkJackson2003_PCdry_forbs" } else { "FILL" }, tree = if (veg_cover[["Rel_Abundance_L0"]][["Trees"]] > 0) { "Bradfordetal2014_LodgepolePine" } else { "FILL" } ) veg_roots <- rSOILWAT2::estimate_PotNatVeg_roots( layers_depth = soil_new[, "depth_cm"], trco_type_by_veg = trco_type_by_veg, fgrass_c3c4ann = veg_cover[["Grasses"]] ) # Add rooting profile to soil v1 <- c("Grass", "Shrub", "Tree", "Forb") v2 <- paste0("transp", v1, "_frac") soil_new[, v2] <- veg_roots[, v1]
rSOILWAT2::swSoils_Layers(sw_in) <- soil_new
Note: This step is usually not necessarily;
instead, the default setup is that parameters are estimated automatically
by SOILWAT2
with the selected pedotransfer function.
tmp <- rSOILWAT2::swSite_SWRCflags(sw_in) swrcp <- rSOILWAT2::ptf_estimate( sand = soil_new[, "sand_frac"], clay = soil_new[, "clay_frac"], fcoarse = soil_new[, "gravel_content"], bdensity = soil_new[, "bulkDensity_g.cm.3"], swrc_name = tmp["swrc_name"], ptf_name = tmp["ptf_name"] ) # Set estimated `SWRC` parameter values rSOILWAT2::swSoils_SWRCp(sw_in) <- swrcp # Declare that parameter values are already estimated rSOILWAT2::swSite_hasSWRCp(sw_in) <- TRUE # Alternatively (and equivalently), # set soil properties and SWRC parameters with one assignment rSOILWAT2::set_swSoils(sw_in) <- list(Layers = soil_new, SWRCp = swrcp)
# Prepare transpiration regions based on soil layers # Adjust to your specific soil depths, # e.g., if extracted from NRCS SDA without updating soil layers tr_example <- c(1, 1, 2, 3, 3) # Values corresponding to `CONUSSOIL_BSE_EVERY10cm` of "rSFSW2": tr_lyrs_10cm <- c(1, 1, 1, 2, 2, 2, 2, 3, 3, 3) tr <- rSOILWAT2::prepare_TranspirationRegions(tr_lyrs = tr_lyrs_10cm) rSOILWAT2::swSite_TranspirationRegions(sw_in) <- tr # Make necessary adjustments based on soil depth and rooting profiles rSOILWAT2::swSite_TranspirationRegions(sw_in) <- rSOILWAT2::adjust_TranspirationRegions(sw_in) # Check that all is fine stopifnot(rSOILWAT2::check_TranspirationRegions(sw_in))
sw_in2 <- sw_in # Example: 50% reduction in June-Aug precipitation rSOILWAT2::swWeather_MonScalingParams(sw_in2)[6:8, "PPT"] <- 0.5 # Example: 2 C warming rSOILWAT2::swWeather_MonScalingParams(sw_in2)[, c("MaxT", "MinT")] <- 2
Note: Knitting vignette crashes with C-level messages -> run quietly
sw_out <- rSOILWAT2::sw_exec( inputData = sw_in, weatherList = wdata, quiet = TRUE ) sw_out2 <- rSOILWAT2::sw_exec( inputData = sw_in2, weatherList = wdata, quiet = TRUE )
res_transp <- rSOILWAT2::get_transpiration(sw_out, "Year") res_transp2 <- rSOILWAT2::get_transpiration(sw_out2, "Year") plot( res_transp, res_transp2, xlab = "Annual transpiration [mm]", ylab = "Warm & summer-dry: annual transpiration [mm]", xlim = c(0, 200), ylim = c(0, 200) ) abline(0, 1, col = "red")
# Task: use function `getWeatherData_folders` to read weather data with # missing values from folder `data_weather_missing` wdata_gaps <- rSOILWAT2::getWeatherData_folders( LookupWeatherFolder = file.path(path_demo, "Input"), weatherDirName = "data_weather_missing" ) # This will fail: missing data sw_out <- try(rSOILWAT2::sw_exec( inputData = sw_in0, weatherList = wdata_gaps ))
(i) Estimate weather generator coefficients with
function dbW_estimate_WGen_coefs()
,
(ii) generate weather with function dbW_generateWeather()
, and
(iii) run SOILWAT2 with filled-in data
wgen_coeffs <- rSOILWAT2::dbW_estimate_WGen_coefs( weatherData = wdata, imputation_type = "locf" ) wdata_filled <- rSOILWAT2::dbW_generateWeather( weatherData = wdata_gaps, years = seq( from = rSOILWAT2::swYears_StartYear(sw_in0), to = rSOILWAT2::swYears_EndYear(sw_in0) ), wgen_coeffs = wgen_coeffs, seed = 123 ) sw_out <- rSOILWAT2::sw_exec( inputData = sw_in0, weatherList = wdata_filled )
(i) estimate weather generator coefficients and (ii) run SOILWAT2 with activated weather generator
rSOILWAT2::swWeather_UseMarkov(sw_in0) <- TRUE rSOILWAT2::swMarkov_Prob(sw_in0) <- wgen_coeffs[["mkv_doy"]] rSOILWAT2::swMarkov_Conv(sw_in0) <- wgen_coeffs[["mkv_woy"]] sw_out <- rSOILWAT2::sw_exec( inputData = sw_in0, weatherList = wdata_gaps )
# Task: plot a comparison between `wdata` and `wdata_filled` for # mean monthly precipitation and mean monthly temperature using the output # of function `calc_SiteClimate` clim1 <- rSOILWAT2::calc_SiteClimate(wdata) clim2 <- rSOILWAT2::calc_SiteClimate(wdata_filled) plot(clim1[["meanMonthlyPPTcm"]], clim2[["meanMonthlyPPTcm"]]) abline(0, 1, col = "red") plot(clim1[["meanMonthlyTempC"]], clim2[["meanMonthlyTempC"]]) abline(0, 1, col = "red")
# Task: compare weather data `wdata` with `wdata_filled` using function # `compare_weather` rSOILWAT2::compare_weather( ref_weather = rSOILWAT2::dbW_weatherData_to_dataframe(wdata), weather = rSOILWAT2::dbW_weatherData_to_dataframe(wdata_filled), N = 1, path = tempdir(), tag = "test" )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.