Merge branch 'evi_helpers' of https://github.com/abbylewis/pecan into evi_helpers

Author: abbylewis

book_source/03_topical_pages/06_data/01_meteorology.Rmd

@@ -153,8 +153,8 @@ Availability: 1950-present
 
 Notes: 
 
-It's important to know that the raw ERA5 tiles needs to be downloaded and registered in the database first. Inside the `inst` folder in the data.atmosphere package there are R files for downloading and registering files in the BETY. However, it assumes that you have registered and setup your API requirements. Check out how to setup your API [here] (https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5). 
-In the `inst` folder you can find two files (`ERA5_db_register.R` and `ERA5_USA_download.R`). If you setup your `ecmwf` account as it's explained in the link above, `ERA5_USA_download.R` will help you to download all the tiles with all the variables required for pecan `extract.nc.ERA5` function to generate pecan standard met files. Besides installing the required packages for this file, it should work from top to bottom with no problem. After downloading the tiles, there is simple script in `ERA5_db_register.R` which helps you register your tiles in the bety. `met.process` later on uses that entery to find the required tiles for extracting met data for your sites. There are important points about this file. 1- Make sure you don't change the site id in the script (which is the same the `ParentSite` in ERA5 registeration xml file). 2- Make sure the start and end date in that script matches the downloaded tiles. Set your `ERA5.files.path` to where you downloaded the tiles and then the rest of the script should be working fine.
+It's important to know that the raw ERA5 tiles needs to be downloaded and registered in the database first. Inside the `inst` folder in the data.atmosphere package there are R files for downloading and registering files in the BETY. However, it assumes that you have registered and setup your API requirements. The download uses the `ecmwfr` package, so you need the UID and API key from your [CDS profile page](https://cds.climate.copernicus.eu/profile), which you pass as the `user` and `key` arguments. Note that `ecmwfr` stores the key using `keyring`; on a headless machine there is no OS keyring and the default backend will prompt for a password, so set `Sys.setenv(R_KEYRING_BACKEND = "env")` first.
+In the `inst` folder you can find two files (`ERA5_db_register.R` and `ERA5_NA_download.R`). Once your CDS credentials are set up, `ERA5_NA_download.R` will help you to download all the tiles with all the variables required for pecan `extract.nc.ERA5` function to generate pecan standard met files. Besides installing the required packages for this file, it should work from top to bottom with no problem. After downloading the tiles, there is simple script in `ERA5_db_register.R` which helps you register your tiles in the bety. `met.process` later on uses that entery to find the required tiles for extracting met data for your sites. There are important points about this file. 1- Make sure you don't change the site id in the script (which is the same the `ParentSite` in ERA5 registeration xml file). 2- Make sure the start and end date in that script matches the downloaded tiles. Set your `ERA5.files.path` to where you downloaded the tiles and then the rest of the script should be working fine.
 
 ## ICOS Drought 2018
 

models/sipnet/NEWS.md

@@ -1,5 +1,16 @@
 # PEcAn.SIPNET 1.10.0.9000
 
+* Improvements to the job.sh written by `write_segmented_configs`:
+  - Now places README.txt, segments.csv, and the full log files from each
+    segment, into the outdir (as was already done for one-segment runs).
+  - Now respects `settings$model$delete.raw` by deleting sipnet.out from segment
+    dirs as well as the from the job outdir
+  - No longer calls model2netcdf.SIPNET() inside each segment dir (it was
+    redundant with the whole-job netcdf output).
+  - Fixed incorrect run paths in the settings$host block passed to single segments
+* `model2netcdf.SIPNET` takes `LAI` from sipnet.out if present (which requires
+    Sipnet > v2.2). If it is not present, LAI is calculated as
+    `plantLeafC / leafCSpWt` as previously.
 * `split_inputs.SIPNET` now avoids internal time format conversions, giving a
   substantial speedup and reduced memory use when processing multi-year files.
 * `model2netcdf.SIPNET` now detects the number of timesteps per day by taking the maximum count across all days in the first simulation year, rather than reading only from day 1. This prevents a factor-of-N error in flux unit conversions when the first day of output is partial (fewer timesteps than a complete day) (#3624, #3989).

models/sipnet/R/model2netcdf.SIPNET.R

@@ -150,16 +150,21 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     )
 
 
-  # calculate LAI for standard output
-  # LAI = plantLeafC / leafCSpWt
-  # both operands are in carbon units (gC/m2 and gC/m2_leaf),
-  # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
-  param <- utils::read.table(file.path(gsub(pattern = "/out/",
-                                            replacement = "/run/", x = outdir),
-                                       "sipnet.param"), stringsAsFactors = FALSE)
-  leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
-  SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
-
+  if (!("LAI" %in% colnames(sipnet_output))) {
+    # When LAI not reported, calculate it from leaf C and leaf specific weight,
+    # with the latter read from the parameter file.
+    # Note the hardcoded + undocumented assumption that ../run/sipnet.param
+    # exists and contains the values that were used to generate this output.
+    # LAI = plantLeafC [gC/m2] / leafCSpWt [gC/m2 leaf] = m2 leaf / m2
+    # both operands are in carbon units (gC/m2 and gC/m2_leaf),
+    # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
+    param <- utils::read.table(file.path(gsub(pattern = "/out/",
+                                              replacement = "/run/", x = outdir),
+                                         "sipnet.param"), stringsAsFactors = FALSE)
+    leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
+    SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
+    sipnet_output$LAI = sipnet_output$plantLeafC * SLA
+  }
 
   ### Loop over years in SIPNET output to create separate netCDF outputs
   for (y in year_seq) {
@@ -208,7 +213,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       "litter_carbon_content" = sub.sipnet.output$litter,
       "fine_root_carbon_content" = sub.sipnet.output$fineRootC,
       "coarse_root_carbon_content" = sub.sipnet.output$coarseRootC,
-      "LAI" = sub.sipnet.output$plantLeafC * SLA,
+      "LAI" = sub.sipnet.output$LAI,
       "TotLivBiom" = sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC +
                        sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC,
       "TotSoilCarb" = sub.sipnet.output$soil + sub.sipnet.output$litter,

models/sipnet/R/write.configs.SIPNET.R

@@ -121,6 +121,16 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     has_microbeInit = rev_str == "v1"
   )
 
+
+  # find out where to write run/ouput
+  rundir <- file.path(settings$host$rundir, as.character(run.id))
+  outdir <- file.path(settings$host$outdir, as.character(run.id))
+  if (is.null(settings$host$qsub) && (settings$host$name == "localhost")) {
+    rundir <- file.path(settings$rundir, as.character(run.id))
+    outdir <- file.path(settings$modeloutdir, as.character(run.id))
+  }
+
+
   ### WRITE sipnet.in
   template.in <- system.file(
     paste0("sipnet.in_", rev_str),
@@ -140,7 +150,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   config.text <- update_flag_lines(config.text, user_flags)
 
-  writeLines(config.text, con = file.path(settings$rundir, run.id, "sipnet.in"))
+  writeLines(config.text, con = file.path(rundir, "sipnet.in"))
 
   ### WRITE *.clim
   template.clim <- settings$run$inputs$met$path  ## read from settings
@@ -152,14 +162,6 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   PEcAn.logger::logger.info(paste0("Writing SIPNET configs with input ", template.clim))
 
-  # find out where to write run/ouput
-  rundir <- file.path(settings$host$rundir, as.character(run.id))
-  outdir <- file.path(settings$host$outdir, as.character(run.id))
-  if (is.null(settings$host$qsub) && (settings$host$name == "localhost")) {
-    rundir <- file.path(settings$rundir, as.character(run.id))
-    outdir <- file.path(settings$modeloutdir, as.character(run.id))
-  }
-
   # create launch script (which will create symlink)
   if (!is.null(settings$model$jobtemplate) && file.exists(settings$model$jobtemplate)) {
     jobsh <- readLines(con = settings$model$jobtemplate, n = -1)
@@ -251,8 +253,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   }
   jobsh <- gsub("@DELETE.RAW@", settings$model$delete.raw, jobsh)
 
-  writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh"))
-  Sys.chmod(file.path(settings$rundir, run.id, "job.sh"))
+  writeLines(jobsh, con = file.path(rundir, "job.sh"))
+  Sys.chmod(file.path(rundir, "job.sh"))
 
 
   ### Copy event file
@@ -268,7 +270,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   ### WRITE *.param-spatial
   if (caps$has_param_spatial) {
     template.paramSpatial <- system.file("template.param-spatial", package = "PEcAn.SIPNET")
-    file.copy(template.paramSpatial, file.path(settings$rundir, run.id, "sipnet.param-spatial"))
+    file.copy(template.paramSpatial, file.path(rundir, "sipnet.param-spatial"))
   }
 
   ### WRITE *.param
@@ -1028,12 +1030,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     }
 
   }
-  if (file.exists(file.path(settings$rundir, run.id, "sipnet.param"))) {
+  if (file.exists(file.path(rundir, "sipnet.param"))) {
     file.rename(
-      file.path(settings$rundir, run.id, "sipnet.param"),
+      file.path(rundir, "sipnet.param"),
       file.path(
-        settings$rundir,
-        run.id,
+        rundir,
         paste0("sipnet_", lubridate::year(settings$run$start.date), "_", lubridate::year(settings$run$end.date), ".param")
       )
     )
@@ -1042,7 +1043,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
 
   utils::write.table(
     param,
-    file.path(settings$rundir, run.id, "sipnet.param"),
+    file.path(rundir, "sipnet.param"),
     row.names = FALSE,
     col.names = FALSE,
     quote = FALSE

models/sipnet/R/write_segmented_configs.R

@@ -97,6 +97,15 @@ write_segmented_configs.SIPNET <- function(settings, input_design = NULL, ...) {
         by = "ens_num",
         relationship = "many-to-one")
   }
+  # Use a custom job.sh template with no model2netcdf step;
+  # we convert all segments together at the end of the timeseries instead.
+  # ...unless user provided their own template, in which case use it as-is
+  if (is.null(settings$model$jobtemplate) || !file.exists(settings$model$jobtemplate)) {
+    settings$model$jobtemplate <- system.file(
+      "template_singlesegment.job",
+      package = "PEcAn.SIPNET"
+    )
+  }
 
   new_jobfiles <- character()
 
@@ -218,6 +227,7 @@ write_segment_configs <- function(
   )
 
   jobsh_files <- character()
+  job_logfiles <- character()
 
   for (isegment in seq_len(nrow(segments))) {
     segment <- segments[isegment, ]
@@ -259,6 +269,18 @@ write_segment_configs <- function(
     segment_settings[["outdir"]] <- segment_outdir
     segment_settings[["modeloutdir"]] <- segment_outdir
     segment_settings[["rundir"]] <- segment_rundir
+    segment_settings[[c("host", "rundir")]] <-  sub(
+      pattern = settings$rundir,
+      replacement = settings$host$rundir,
+      x = segment_rundir,
+      fixed = TRUE
+    )
+    segment_settings[[c("host", "outdir")]] <-  sub(
+      pattern = settings$modeloutdir,
+      replacement = settings$host$outdir,
+      x = segment_outdir,
+      fixed = TRUE
+    )
     segment_settings[[c("run", "start.date")]] <- dstart
     segment_settings[[c("run", "end.date")]] <- dend
     segment_settings[[c("run", "inputs")]] <- segment_inputs
@@ -291,12 +313,19 @@ write_segment_configs <- function(
       run.id = runid_dummy
     )
 
+    segment_log <- file.path(segment_settings$modeloutdir,
+                             runid_dummy,
+                             "logfile.txt")
     segment_jobsh <- file.path(segment_settings$rundir, runid_dummy, "job.sh")
     stopifnot(file.exists(segment_jobsh))
+    job_logfiles <- c(job_logfiles, segment_log)
     jobsh_files <- c(jobsh_files, segment_jobsh)
   }
 
   # Now, get the run's jobsh file
+  # NB this doesn't do any path expansion or @TEMPLATE@ string replacement --
+  # those are done in the per-segment job.sh.
+  # Note especially that means no @HOST_SETUP@ / @HOST_TEARDOWN@ in this script
   run_jobsh <- file.path(run_dir, "job.sh")
   target_sipnet_out <- file.path(run_modeloutdir, "sipnet.out")
   segmented_jobsh_file <- file.path(run_dir, "job_segmented.sh")
@@ -326,10 +355,31 @@ write_segment_configs <- function(
         sprintf("sitelon = %s", as.character(settings$run$site$lon)),
         sprintf("start_date = %s", shQuote(settings$run$start.date)),
         sprintf("end_date = %s", shQuote(settings$run$end.date)),
+        sprintf("delete.raw = %s", shQuote(settings$model$delete.raw)),
         sprintf("revision = %s", shQuote(settings$model$revision)),
         sep = ", "
       )
-    )
+    ),
+    "",
+    "# copy readme with specs to output",
+    paste("cp", file.path(run_dir, "README.txt"), file.path(run_modeloutdir, "README.txt")),
+    paste("cp", file.path(run_dir, "segments.csv"), file.path(run_modeloutdir, "segments.csv")),
+    "",
+    "# Concatenate segment log files",
+    paste("echo \"\n--> contents of", job_logfiles, ":\" && cat", job_logfiles),
+    "",
+    if (isTRUE(as.logical(settings$model$delete.raw))) {
+      # Removing all outputs that have been safely copied to the outdir
+      c(
+        "# Remove per-segment outputs & logs after concatenating to job outdir",
+        sprintf("find %s -name sipnet.out -delete", segment_rootdir),
+        sprintf("find %s -name logfile.txt -delete", segment_rootdir),
+        sprintf("rm %s/README.txt", run_dir),
+        sprintf("rm %s/segments.csv", run_dir)
+      )
+    },
+    "",
+    "echo -e \"MODEL FINISHED\nLogfile is located at '${OUTDIR}/logfile.txt'\" >&3"
   )
   writeLines(segmented_jobsh_lines, segmented_jobsh_file)
   if (replace_and_link) {

models/sipnet/inst/template_singlesegment.job

@@ -0,0 +1,53 @@
+#!/bin/bash
+
+# redirect output
+exec 3>&1
+exec &> "@OUTDIR@/logfile.txt"
+
+# host specific setup
+@HOST_SETUP@
+
+# cdo setup
+@CDO_SETUP@
+
+# create output folder
+mkdir -p "@OUTDIR@"
+
+# Convert any relative paths to absolute
+# (otherwise we'd lose track of them when cd'ing into rundir)
+OUTDIR=$(cd "@OUTDIR@" && pwd -P)
+RUNDIR=$(cd "@RUNDIR@" && pwd -P)
+SITE_MET=$(cd $(dirname "@SITE_MET@") && pwd -P)/$(basename "@SITE_MET@")
+BINARY=$(cd $(dirname "@BINARY@") && pwd -P)/$(basename "@BINARY@")
+
+# see if application needs running
+if [ ! -e "${OUTDIR}/sipnet.out" ]; then
+  cd "$RUNDIR"
+  ln -s "${SITE_MET}" sipnet.clim
+
+  "${BINARY}"
+  STATUS=$?
+  
+  # copy output
+  mv "${RUNDIR}/sipnet.out" "$OUTDIR"
+
+  # check the status
+  if [ $STATUS -ne 0 ]; then
+    echo -e "ERROR IN MODEL RUN\nLogfile is located at '${OUTDIR}/logfile.txt'" >&3
+    exit $STATUS
+  fi
+
+  # conversion to MsTMIP not done in this single-segment step.
+  # The whole-job script will handle conversion after all segments have run.
+fi
+
+# copy readme with specs to output
+cp  "${RUNDIR}/README.txt" "${OUTDIR}/README.txt"
+
+# run getdata to extract right variables
+
+# host specific teardown
+@HOST_TEARDOWN@
+
+# all done
+echo -e "MODEL FINISHED\nLogfile is located at '${OUTDIR}/logfile.txt'" >&3

models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R

@@ -291,3 +291,22 @@ test_that("out_day is correct when the first day of output is partial", {
   ts_correct <- 86400 / 2
   expect_equal(gpp, rep(base_row$gpp * 1e-3 / ts_correct, 5), tolerance = 1e-10)
 })
+
+test_that("LAI taken from leafCSpWt if not present in output", {
+  dat <- make_v2_sipnet(2)
+  # setup_sipnet_test sets leafCSpWt = 32 => expect LAI of 1 and 10
+  dat$plantLeafC = c(32, 320)
+
+  # First without LAI in dat => looks for sipnet.param in run dir
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc1 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc1), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc1, "LAI")), c(1, 10))
+
+  # Now with LAI specified => uses it as-is
+  dat$LAI = c(2, 3)
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc2 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc2), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc2, "LAI")), c(2, 3))
+})

modules/data.atmosphere/R/ERA5_download.R

@@ -9,6 +9,12 @@
 #' To get a Copernicus CDS API key, register at \url{https://cds.climate.copernicus.eu/profile}.
 #' You must provide both \code{user} (UID) and \code{key} parameters from your CDS profile.
 #'
+#' Credentials are stored by \code{ecmwfr::wf_set_key}, which uses the \code{keyring} package.
+#' On a headless machine (cluster node, server, container) there is no OS keyring available and
+#' the default backend will prompt for a keyring password. To avoid this, set
+#' \code{Sys.setenv(R_KEYRING_BACKEND = "env")} before calling this function, which keeps the
+#' key in an environment variable for the session instead.
+#'
 #' You can check the "CC-BY" license under the \href{https://cds.climate.copernicus.eu/profile?tab=licences}{'licences' tab of your profile page}.
 #' @param outfolder Character. Directory where downloaded NetCDF files will be saved.
 #' @param start_date character: the start date of the data to be downloaded. Format is YYYY-MM-DD (will only use the year part of the date)

modules/data.atmosphere/inst/ERA5/ERA5_NA_download.R

@@ -2,6 +2,12 @@ outfolder <- "/projectnb/dietzelab/dongchen/anchorSites/ERA5/"
 start_date <- "2012-01-01"
 end_date <- "2021-12-31"
 extent <- c(-179, -20, 7, 85)
+# cds credentials, from https://cds.climate.copernicus.eu/profile
+user <- "your_cds_uid"
+key <- "your_cds_api_key"
+# ecmwfr keeps the key in the keyring. a headless machine has no os keyring and
+# the default backend prompts for a password, so keep the key in the environment
+Sys.setenv(R_KEYRING_BACKEND = "env")
 variables <- c("2m_temperature",
                "surface_pressure",
                "2m_dewpoint_temperature",
@@ -15,4 +21,5 @@ results <- PEcAn.data.atmosphere::download.ERA5_cds(outfolder = outfolder,
                                                     end_date = end_date,
                                                     extent = extent,
                                                     variables = variables,
-                                                    auto.create.key = T)
+                                                    user = user,
+                                                    key = key)