In Burke-Lauenroth-Lab/rSOILWAT2: An Ecohydrological Ecosystem-Scale Water Balance Simulation Model

\pagebreak

Setup a weather database

Example data set with three sites

  # Load rSOILWAT2
  library("rSOILWAT2")

  # Load data for an example site: see ?sw_exampleData
  sw_in0 <- rSOILWAT2::sw_exampleData
  tmp <- swSite_IntrinsicSiteParams(sw_in0)

  sites1 <- data.frame(
    Longitude = rep(tmp[["Longitude"]], 3),
    Latitude = rep(tmp[["Latitude"]], 3),
    Label = paste0("ExampleSite_", seq_len(3))
  )

  scenario1 <- "Reference"

Create a new weather database

  my_dbW <- "dbWeatherData.sqlite3"

  is_dbW_success <- dbW_createDatabase(
    dbFilePath = my_dbW,
    site_data = sites1,
    Scenarios = scenario1,
    scen_ambient = scenario1
  )

  # Check that database was successfully created
  stopifnot(is_dbW_success)

We will see later how to add more sites and scenarios

Connect to the newly created weather database and check if it is valid

  dbW_setConnection(dbFilePath = my_dbW)

  # Check that connection to weather database was successfully established
  stopifnot(dbW_IsValid())

Query site and scenario information

  # Read entire site table and compare to our inputs
  dbW_sites <- dbW_getSiteTable()
  stopifnot(all.equal(dbW_sites[, colnames(sites1)], sites1))

  # Read entire scenario table and compare to our inputs
  dbW_scenarios <- dbW_getScenariosTable()
  stopifnot(all.equal(dbW_scenarios[, "Scenario"], scenario1))

Put together daily meteorological data and add to weather database

Here, we create three replicates of our example data with the weather generator

  # Weather from our example data
  wdata <- get_WeatherHistory(sw_in0)
  years <- get_years_from_weatherData(wdata)

  # Estimate weather generator coefficients
  wgen_coeffs <- dbW_estimate_WGen_coefs(
    weatherData = wdata,
    imputation_type = "mean",
    imputation_span = 5
  )

  x_empty <- weatherHistory()

  # Loop over our three sites, generate weather, and add it to the database
  for (k in seq_len(nrow(sites1))) {
    # Generate weather
    tmp_wdata <- dbW_generateWeather(
      weatherData = x_empty,
      years = years,
      wgen_coeffs = wgen_coeffs,
      seed = k
    )

    # Add to weather database
    tmp_add <- dbW_addWeatherData(
      Label = sites1[k, "Label"],
      weatherData = tmp_wdata,
      Scenario = scenario1
    )

    # Check that weather data was successfully added to the weather database
    stopifnot(tmp_add)
  }

Work with the weather database

Query database for weather data

  x <- dbW_getWeatherData(
    Label = sites1[2, "Label"],
    Scenario = scenario1[1]
  )

  # Weather data is organized in a list
  # where each element, representing daily data for one Gregorian calendar year,
  # is of class "swWeatherData"
  stopifnot(
    inherits(x, "list"),
    sapply(x, inherits, what = "swWeatherData")
  )

  # Run a simulation with that weather object
  sw_out <- sw_exec(inputData = sw_in0, weatherList = x)

Add two additional sites

  # Information for additional sites
  sites2 <- data.frame(
    Longitude = rep(sites1[1, "Longitude"], 2),
    Latitude = rep(sites1[1, "Latitude"], 2),
    Label = paste0("ExampleSite_", nrow(sites1) + seq_len(2))
  )

  # Add sites to existing weather database
  tmp_add <- dbW_addSites(site_data = sites2)
  stopifnot(tmp_add)

  # There is no weather data associated with the new sites yet!
  try(
    dbW_getWeatherData(Label = sites2[1, "Label"], Scenario = scenario1[1]),
    silent = TRUE
  )

Add weather data for the additional sites

  # Loop over sites, generate weather, and add it to the database
  for (k in seq_len(nrow(sites2))) {
    # Generate weather
    tmp_wdata <- dbW_generateWeather(
      weatherData = x_empty,
      years = years,
      wgen_coeffs = wgen_coeffs,
      seed = k
    )

    # Add to weather database
    tmp_add <- dbW_addWeatherData(
      Label = sites2[k, "Label"],
      weatherData = tmp_wdata,
      Scenario = scenario1
    )

    # Check that weather data was successfully added to the weather database
    stopifnot(tmp_add)
  }

  # Now, we can query the weather at the new sites
  x <- dbW_getWeatherData(
    Label = sites2[1, "Label"],
    Scenario = scenario1[1]
  )
  stopifnot(
    inherits(x, "list"),
    sapply(x, inherits, what = "swWeatherData")
  )

Add a new climate scenario to the weather database

  # Information for additional scenario
  scenario2 <- "+2C"

  # Add scenario to weather database
  tmp_add <- dbW_addScenarios(scenario2)
  stopifnot(tmp_add)

  # All scenarios
  scenarios <- c(scenario1, scenario2)

Add weather data to all sites for the new scenario

Here, we use again the weather generator to generate weather that is 2 C warmer

  # Add 2 C to our weather data
  tmp <- dbW_weatherData_to_dataframe(wdata)
  var_temp <- c("Tmax_C", "Tmin_C")
  tmp[, var_temp] <- tmp[, var_temp] + 2
  wdata_2C <- dbW_dataframe_to_weatherData(tmp)

  # Estimate weather generator coefficients
  wgen_coeffs_2C <- dbW_estimate_WGen_coefs(
    weatherData = wdata_2C,
    imputation_type = "mean",
    imputation_span = 5
  )

  # Loop over our three sites, generate weather, and add it to the database
  sites <- rbind(sites1, sites2)

  for (k in seq_len(nrow(sites))) {
    # Generate weather
    tmp_wdata <- dbW_generateWeather(
      weatherData = x_empty,
      years = years,
      wgen_coeffs = wgen_coeffs_2C,
      seed = k
    )

    # Add to weather database
    tmp_add <- dbW_addWeatherData(
      Label = sites[k, "Label"],
      weatherData = tmp_wdata,
      Scenario = scenario2
    )

    # Check that weather data was successfully added to the weather database
    stopifnot(tmp_add)
  }

Plot daily climate means for all sites and scenarios

  vars <- c("Tmax_C", "Tmin_C", "PPT_cm")

  x <- array(
    data = NA,
    dim = c(365, length(vars), nrow(sites), length(scenarios)),
    dimnames = list(NULL, vars, sites[, "Label"], scenarios)
  )

  used_doys <- seq_len(365)

  # Query weather data for all sites and scenarios and calculate means
  for (k1 in seq_len(nrow(sites))) for (k2 in seq_along(scenarios)) {
    tmp <- dbW_getWeatherData(
      Label = sites[k1, "Label"],
      Scenario = scenarios[k2]
    )

    tmp_df <- dbW_weatherData_to_dataframe(tmp)
    tmp_clim <- aggregate(
      x = tmp_df[, vars],
      by = list(DOY = tmp_df[, "DOY"]),
      FUN = mean
    )

    x[, vars, k1, k2] <- as.matrix(tmp_clim[used_doys, vars])
  }


  # Calculate overall means per variable and scenario
  round(apply(x, c(2, 4), mean, na.rm = TRUE), 2)


  # Plot mean climate variables
  colsc <- c("orange", "purple")
  colsc2 <- adjustcolor(colsc, alpha.f = 0.4)

  par_prev <- par(
    mfrow = c(2, 2),
    mar = c(2.5, 2.5, 0.5, 0.5), mgp = c(1, 0, 0),
    tcl = 0.3
  )

  for (k1 in seq_along(vars)) {
    for (k2 in seq_along(scenarios)) {
      if (k2 == 1) {
        matplot(
          x[, k1, , k2],
          type = "l", lty = 1, lwd = 0.5, col = colsc2[k2],
          ylim = range(x[, k1, , ], na.rm = TRUE),
          xlab = "Time (day of year)", ylab = vars[k1]
        )
      } else {
        matlines(x[, k1, , k2], lty = 1, lwd = 0.5, col = colsc2[k2])
      }
    }

    if (k1 == 1) {
      legend(
        "bottom",
        bty = "n",
        legend = scenarios,
        col = colsc,
        lwd = 2
      )
    }
  }

  par(par_prev)

Disconnect and clean up

  dbW_disconnectConnection()
  unlink(my_dbW)

Burke-Lauenroth-Lab/rSOILWAT2 documentation built on Dec. 9, 2023, 1:46 a.m.