From 582f9b5959513aa13d58fb4b96d8e79a71522425 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 3 Jan 2025 21:55:03 -0800 Subject: [PATCH 1/4] move duplicate code out of both legs of if/else, line wrapping --- base/workflow/R/runModule.get.trait.data.R | 57 ++++++------------- base/workflow/R/runModule.run.write.configs.R | 23 +++++--- 2 files changed, 30 insertions(+), 50 deletions(-) diff --git a/base/workflow/R/runModule.get.trait.data.R b/base/workflow/R/runModule.get.trait.data.R index 3b312dcf330..ce084567714 100644 --- a/base/workflow/R/runModule.get.trait.data.R +++ b/base/workflow/R/runModule.get.trait.data.R @@ -18,54 +18,29 @@ runModule.get.trait.data <- function(settings) { pfts <- c(pfts, pfts.i[ind]) pft.names <- sapply(pfts, function(x) x$name) } - - PEcAn.logger::logger.info(paste0( - "Getting trait data for all PFTs listed by any Settings object in the list: ", + PEcAn.logger::logger.info( + "Getting trait data for all PFTs listed by any Settings object in the list:", paste(pft.names, collapse = ", ") - )) - - modeltype <- settings$model$type - dbfiles <- settings$database$dbfiles - database <- settings$database$bety - forceupdate <- - ifelse(is.null(settings$meta.analysis$update), - FALSE, - settings$meta.analysis$update) - write <- settings$database$bety$write - settings$pfts <- - PEcAn.DB::get.trait.data( - pfts = pfts, - modeltype = modeltype, - dbfiles = dbfiles, - database = database, - forceupdate = forceupdate, - write = write - ) - return(settings) + ) } else if (PEcAn.settings::is.Settings(settings)) { pfts <- settings$pfts if (!is.list(pfts)) { PEcAn.logger::logger.severe("settings$pfts is not a list") } - modeltype <- settings$model$type - dbfiles <- settings$database$dbfiles - database <- settings$database$bety - forceupdate <- - ifelse(is.null(settings$meta.analysis$update), - FALSE, - settings$meta.analysis$update) - write <- settings$database$bety$write - settings$pfts <- - PEcAn.DB::get.trait.data( - pfts = pfts, - modeltype = modeltype, - dbfiles = dbfiles, - database = database, - forceupdate = forceupdate, - write = write - ) - return(settings) } else { stop("runModule.get.trait.data only works with Settings or MultiSettings") } + + settings$pfts <- PEcAn.DB::get.trait.data( + pfts = pfts, + modeltype = settings$model$type, + dbfiles = settings$database$dbfiles, + database = settings$database$bety, + forceupdate = ifelse(is.null(settings$meta.analysis$update), + FALSE, + settings$meta.analysis$update), + write = settings$database$bety$write + ) + + return(settings) } diff --git a/base/workflow/R/runModule.run.write.configs.R b/base/workflow/R/runModule.run.write.configs.R index 265da4f6c1c..61b7d0bed28 100644 --- a/base/workflow/R/runModule.run.write.configs.R +++ b/base/workflow/R/runModule.run.write.configs.R @@ -15,19 +15,24 @@ runModule.run.write.configs <- function(settings, overwrite = TRUE) { } else if (PEcAn.settings::is.Settings(settings)) { write <- settings$database$bety$write # double check making sure we have method for parameter sampling - if (is.null(settings$ensemble$samplingspace$parameters$method)) settings$ensemble$samplingspace$parameters$method <- "uniform" - ens.sample.method <- settings$ensemble$samplingspace$parameters$method - - - #check to see if there are posterior.files tags under pft - posterior.files <-settings$pfts %>% - purrr::map(purrr::possibly('posterior.files', NA_character_)) %>% + if (is.null(settings$ensemble$samplingspace$parameters$method)) { + settings$ensemble$samplingspace$parameters$method <- "uniform" + } + + # check to see if there are posterior.files tags under pft + posterior.files <- settings$pfts %>% + purrr::map(purrr::possibly("posterior.files", NA_character_)) %>% purrr::modify_depth(1, function(x) { ifelse(is.null(x), NA_character_, x) }) %>% unlist() - - return(PEcAn.workflow::run.write.configs(settings, write, ens.sample.method, posterior.files = posterior.files, overwrite = overwrite)) + + return(PEcAn.workflow::run.write.configs( + settings = settings, + write = settings$database$bety$write, + ens.sample.method = settings$ensemble$samplingspace$parameters$method, + posterior.files = posterior.files, + overwrite = overwrite)) } else { stop("runModule.run.write.configs only works with Settings or MultiSettings") } From 578c20c069a2f10665828afe2fba97b2cff814f0 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 3 Jan 2025 22:12:31 -0800 Subject: [PATCH 2/4] ran styler on ERA5 fns --- modules/data.atmosphere/R/ERA5_met_process.R | 217 ++++++++------- modules/data.atmosphere/R/extract_ERA5.R | 252 +++++++++--------- .../data.atmosphere/man/ERA5_met_process.Rd | 3 +- .../data.atmosphere/man/extract.nc.ERA5.Rd | 6 +- 4 files changed, 244 insertions(+), 234 deletions(-) diff --git a/modules/data.atmosphere/R/ERA5_met_process.R b/modules/data.atmosphere/R/ERA5_met_process.R index d483be8f1b1..4fd8f035272 100644 --- a/modules/data.atmosphere/R/ERA5_met_process.R +++ b/modules/data.atmosphere/R/ERA5_met_process.R @@ -6,62 +6,69 @@ #' @param write.db if write into Bety database #' @param write if write the settings into pecan.xml file in the outdir of settings. #' -#' @return if write.db is True then return input IDs with physical paths; if write.db is False then return just physical paths of extracted ERA5 clim files. +#' @return if write.db is True then return input IDs with physical paths; +#' if write.db is False then return just physical paths of extracted ERA5 clim files. #' @export -#' +#' #' @author Dongchen Zhang #' @importFrom dplyr %>% #' -ERA5_met_process <- function(settings, in.path, out.path, write.db=FALSE, write = TRUE){ - #Initialize the multicore computation. +ERA5_met_process <- function(settings, in.path, out.path, write.db = FALSE, write = TRUE) { + # Initialize the multicore computation. if (future::supportsMulticore()) { future::plan(future::multicore) } else { future::plan(future::multisession) } - - #getting site info - #grab the site info from Bety DB if we can't get the site info directly from the settings object. - if ("try-error" %in% class(try(site_info <- settings$run %>% - purrr::map('site')%>% - purrr::map(function(site.list){ - #conversion from string to number - site.list$lat <- as.numeric(site.list$lat) - site.list$lon <- as.numeric(site.list$lon) - list(site.id=site.list$id, lat=site.list$lat, lon=site.list$lon, site_name=site.list$name) - }) %>% - dplyr::bind_rows() %>% - as.list()))) { - #getting site ID + + # getting site info + # grab the site info from Bety DB if we can't get the site info directly from the settings object. + if ("try-error" %in% class(try(site_info <- settings$run %>% + purrr::map("site") %>% + purrr::map(function(site.list) { + # conversion from string to number + site.list$lat <- as.numeric(site.list$lat) + site.list$lon <- as.numeric(site.list$lon) + list(site.id = site.list$id, lat = site.list$lat, lon = site.list$lon, site_name = site.list$name) + }) %>% + dplyr::bind_rows() %>% + as.list())) + ) { + # getting site ID observations <- c() for (i in 1:length(settings)) { obs <- settings[[i]]$run$site$id - observations <- c(observations,obs) + observations <- c(observations, obs) } - - #query site info - bety <- dplyr::src_postgres(dbname = settings$database$bety$dbname, - host = settings$database$bety$host, - user = settings$database$bety$user, - password = settings$database$bety$password) + + # query site info + bety <- dplyr::src_postgres( + dbname = settings$database$bety$dbname, + host = settings$database$bety$host, + user = settings$database$bety$user, + password = settings$database$bety$password + ) con <- bety$con site_ID <- observations suppressWarnings(site_qry <- glue::glue_sql("SELECT *, ST_X(ST_CENTROID(geometry)) AS lon, ST_Y(ST_CENTROID(geometry)) AS lat FROM sites WHERE id IN ({ids*})", - ids = site_ID, .con = con)) - suppressWarnings(qry_results <- PEcAn.DB::db.query(con = con, query = site_qry))#use PEcAn.DB instead - site_info <- list(site_id=qry_results$id, site_name=qry_results$sitename, lat=qry_results$lat, - lon=qry_results$lon, time_zone=qry_results$time_zone) + ids = site_ID, .con = con + )) + suppressWarnings(qry_results <- PEcAn.DB::db.query(con = con, query = site_qry)) # use PEcAn.DB instead + site_info <- list( + site_id = qry_results$id, site_name = qry_results$sitename, lat = qry_results$lat, + lon = qry_results$lon, time_zone = qry_results$time_zone + ) } - - #initialize db query elements - if(write.db){ + + # initialize db query elements + if (write.db) { mimetype <- "application/x-netcdf" formatname <- "CF Meteorology" hostname <- PEcAn.remote::fqdn() # find mimetype, if it does not exist, it will create one mimetypeid <- PEcAn.DB::get.id("mimetypes", "type_string", mimetype, con, create = TRUE) - + # find appropriate format, create if it does not exist formatid <- PEcAn.DB::get.id( table = "formats", @@ -71,131 +78,135 @@ ERA5_met_process <- function(settings, in.path, out.path, write.db=FALSE, write create = TRUE, dates = TRUE ) - + # setup parent part of query if specified parent <- "" - - #initialize Input_IDs object when looping over each site + + # initialize Input_IDs object when looping over each site Input_IDs <- list() } - - #restructure the site_info into list. + + # restructure the site_info into list. site_info$start_date <- start_date <- rep(settings$state.data.assimilation$start.date, length(settings)) site_info$end_date <- end_date <- rep(settings$state.data.assimilation$end.date, length(settings)) site_info$out.path <- rep(out.path, length(settings)) site_info$in.path <- rep(in.path, length(settings)) site_info$model.type <- rep(settings$model$type, length(settings)) new.site.info <- split(as.data.frame(site_info), seq(nrow(as.data.frame(site_info)))) - - #Extract ERA5 for each site. + + # Extract ERA5 for each site. PEcAn.logger::logger.info("Started extracting ERA5 data!\n") - Clim_paths <- furrr::future_map(new.site.info, function(site){ - #check if sub-folder exists, if doesn't then create a new folder specific for each site - site_outFolder <- paste0(site$out.path,'/', site$site.id) - #check if folder already exists, if it does, then jump to the next loop - if(!file.exists(site_outFolder)){ + Clim_paths <- furrr::future_map(new.site.info, function(site) { + # check if sub-folder exists, if doesn't then create a new folder specific for each site + site_outFolder <- paste0(site$out.path, "/", site$site.id) + # check if folder already exists, if it does, then jump to the next loop + if (!file.exists(site_outFolder)) { dir.create(site_outFolder) - }else{ - #grab physical paths of existing ERA5 files - #need to be generalized when more models come in. - clim.paths <- list(in.path=list.files(path=site_outFolder, pattern = '*.clim', full.names = T)) + } else { + # grab physical paths of existing ERA5 files + # need to be generalized when more models come in. + clim.paths <- list(in.path = list.files(path = site_outFolder, pattern = "*.clim", full.names = T)) names(clim.paths) <- site$site.id return(clim.paths) } - - #extract ERA5.nc files - PEcAn.data.atmosphere::extract.nc.ERA5(slat = site$lat, - slon = site$lon, - in.path = site$in.path, - start_date = site$start_date, - end_date = site$end_date, - outfolder = site_outFolder, - in.prefix = 'ERA5_', - newsite = as.character(site$site.id)) - - #starting working on met2model.model function over each ensemble - #setting up met2model function depending on model name from settings + + # extract ERA5.nc files + PEcAn.data.atmosphere::extract.nc.ERA5( + slat = site$lat, + slon = site$lon, + in.path = site$in.path, + start_date = site$start_date, + end_date = site$end_date, + outfolder = site_outFolder, + in.prefix = "ERA5_", + newsite = as.character(site$site.id) + ) + + # starting working on met2model.model function over each ensemble + # setting up met2model function depending on model name from settings met2model_method <- do.call("::", list(paste0("PEcAn.", site$model.type), paste0("met2model.", site$model.type))) - #grab the rbind.xts function + # grab the rbind.xts function rbind.xts <- do.call("::", list("xts", "rbind.xts")) - #find every path associated with each ensemble member + # find every path associated with each ensemble member ens_nc <- list.files(path = site_outFolder, full.names = T) - #loop over each ensemble member + # loop over each ensemble member for (i in 1:length(ens_nc)) { nc_path <- ens_nc[i] - - #find a proper in prefix for each ensemble member - ens_num <- strsplit(basename(nc_path),"_")[[1]][3] + + # find a proper in prefix for each ensemble member + ens_num <- strsplit(basename(nc_path), "_")[[1]][3] in_prefix <- paste0("ERA5.", ens_num) - - #preparing for the met2model.SIPNET function - met2model_method(in.path = nc_path, - in.prefix = in_prefix, - outfolder = site_outFolder, - start_date = site$start_date, - end_date = site$end_date) + + # preparing for the met2model.SIPNET function + met2model_method( + in.path = nc_path, + in.prefix = in_prefix, + outfolder = site_outFolder, + start_date = site$start_date, + end_date = site$end_date + ) } # grab physical paths of ERA5 files - clim.paths <- list(in.path=list.files(path=site_outFolder, pattern = '*.clim', full.names = T)) + clim.paths <- list(in.path = list.files(path = site_outFolder, pattern = "*.clim", full.names = T)) names(clim.paths) <- site$site.id return(clim.paths) }, .progress = TRUE) PEcAn.logger::logger.info("\nFinished!") - - #write the paths into settings. + + # write the paths into settings. if (write) { - #write paths into settings. + # write paths into settings. for (i in seq_along(settings)) { - #fill in dates related to met files. - settings[[i]]$run$site$met.start <- - settings[[i]]$run$start.date <- + # fill in dates related to met files. + settings[[i]]$run$site$met.start <- + settings[[i]]$run$start.date <- settings[[i]]$state.data.assimilation$start.date - settings[[i]]$run$site$met.end <- - settings[[i]]$run$end.date <- + settings[[i]]$run$site$met.end <- + settings[[i]]$run$end.date <- settings[[i]]$state.data.assimilation$end.date settings[[i]]$run$inputs$met$path <- as.list(unlist(Clim_paths[[i]])) %>% purrr::set_names(rep("path", length(Clim_paths[[i]]))) } - - #write settings into xml file. + + # write settings into xml file. PEcAn.logger::logger.info(paste0("Write updated pecan.xml file into: ", file.path(settings$outdir, "pecan.xml"))) PEcAn.settings::write.settings(settings, outputfile = "pecan.xml") } - - #write into bety - if(write.db){ + + # write into bety + if (write.db) { PEcAn.logger::logger.info("Write into database!") - #loop over each site + # loop over each site for (i in 1:length(site_info$site_id)) { - #loop over each ensemble - #initialize arrays to store input and dbfile IDs. + # loop over each ensemble + # initialize arrays to store input and dbfile IDs. dbfile_IDs <- c() input_IDs <- c() - for(j in 1:length(Clim_paths[[i]])){ - #create input record for each ensemble member - #insert into inputs table + for (j in 1:length(Clim_paths[[i]])) { + # create input record for each ensemble member + # insert into inputs table cmd <- paste0( "INSERT INTO inputs ", "(site_id, format_id, start_date, end_date, name) VALUES (", - site_info$site_id[i], ", ", formatid, ", '", start_date, "', '", end_date, "','", paste0('ERA5_',site_info$site_id[i],"_",as.character(j)), + site_info$site_id[i], ", ", formatid, ", '", start_date, "', '", end_date, "','", paste0("ERA5_", site_info$site_id[i], "_", as.character(j)), "') RETURNING id" ) # This is the id that we just registered inputid <- PEcAn.DB::db.query(query = cmd, con = con) input_IDs <- c(input_IDs, inputid) - - #create dbfiles associated with each ensemble ID + + # create dbfiles associated with each ensemble ID dbfileid <- PEcAn.DB::dbfile.insert( in.path = Clim_paths[[i]][j], in.prefix = paste0("ERA5.", as.character(j)), type = "Input", id = inputid, con = con, reuse = TRUE, hostname = hostname ) dbfile_IDs <- c(dbfile_IDs, dbfileid) } - Input_IDs[[i]] <- list(input_ID=inputid$id, dbfile_IDs=dbfile_IDs, Site_ID=site_info$site_id[i], in.path=Clim_paths[[i]]) + Input_IDs[[i]] <- list(input_ID = inputid$id, dbfile_IDs = dbfile_IDs, Site_ID = site_info$site_id[i], in.path = Clim_paths[[i]]) } - save(Input_IDs, file=paste0(out.path, '/', 'Inputs.RData')) + save(Input_IDs, file = paste0(out.path, "/", "Inputs.RData")) return(Input_IDs) - }else{ - save(Clim_paths, file=paste0(out.path, '/', 'Inputs.RData')) + } else { + save(Clim_paths, file = paste0(out.path, "/", "Inputs.RData")) return(Clim_paths) } -} \ No newline at end of file +} diff --git a/modules/data.atmosphere/R/extract_ERA5.R b/modules/data.atmosphere/R/extract_ERA5.R index fe5e5b7b5ab..72472ba631c 100644 --- a/modules/data.atmosphere/R/extract_ERA5.R +++ b/modules/data.atmosphere/R/extract_ERA5.R @@ -21,18 +21,17 @@ #' @export #' @examples #' \dontrun{ -#' point.data <- ERA5_extract(sslat=40, slon=-120, years=c(1990:1995), vars=NULL) -#' -# point.data %>% -#' purrr::map(~xts::apply.daily(.x, mean)) +#' point.data <- ERA5_extract(sslat = 40, slon = -120, years = c(1990:1995), vars = NULL) #' +#' # point.data %>% +#' purrr::map(~ xts::apply.daily(.x, mean)) #' } extract.nc.ERA5 <- - function(slat , - slon , - in.path , + function(slat, + slon, + in.path, start_date, - end_date , + end_date, outfolder, in.prefix, newsite, @@ -40,127 +39,126 @@ extract.nc.ERA5 <- overwrite = FALSE, verbose = FALSE, ...) { - # library(xts) - # Distributing the job between whatever core is available. - - years <- seq(lubridate::year(start_date), - lubridate::year(end_date), - 1 + # Distributing the job between whatever core is available. + + years <- seq( + lubridate::year(start_date), + lubridate::year(end_date), + 1 ) ensemblesN <- seq(1, 10) - - - tryCatch({ - #for each ensemble - one.year.out <- years %>% - purrr::map(function(year) { - - # for each year - point.data <- ensemblesN %>% - purrr::map(function(ens) { - - - ncfile <- file.path(in.path, paste0(in.prefix, year, ".nc")) - - #printing out initial information. - if (verbose) { - PEcAn.logger::logger.info(paste0("Trying to open :", ncfile, " ")) - - if (!file.exists(ncfile)) - PEcAn.logger::logger.severe("The nc file was not found.") - - #msg - PEcAn.logger::logger.info(paste0(year, " is being processed ", "for ensemble #", ens, " ")) - } - - #open the file - nc_data <- ncdf4::nc_open(ncfile) - # time stamp - - t <- ncdf4::ncvar_get(nc_data, "time") - tunits <- ncdf4::ncatt_get(nc_data, 'time') - tustr <- strsplit(tunits$units, " ") - timestamp <- - as.POSIXct(t * 3600, tz = "UTC", origin = tustr[[1]][3]) - try(ncdf4::nc_close(nc_data)) - - - # set the vars - if (is.null(vars)) - vars <- names(nc_data$var) - # for the variables extract the data - all.data.point <- vars %>% - purrr::set_names(vars) %>% - purrr::map_dfc(function(vname) { - if (verbose) { - PEcAn.logger::logger.info(paste0(" \t ",vname, "is being extracted ! ")) + + + tryCatch( + { + # for each ensemble + one.year.out <- years %>% + purrr::map(function(year) { + # for each year + point.data <- ensemblesN %>% + purrr::map(function(ens) { + ncfile <- file.path(in.path, paste0(in.prefix, year, ".nc")) + + # printing out initial information. + if (verbose) { + PEcAn.logger::logger.info(paste0("Trying to open :", ncfile, " ")) + + if (!file.exists(ncfile)) { + PEcAn.logger::logger.severe("The nc file was not found.") } - - brick.tmp <- - raster::brick(ncfile, varname = vname, level = ens) - nn <- - raster::extract(brick.tmp, - sp::SpatialPoints(cbind(slon, slat)), - method = 'simple') - if (verbose) { - if (!is.numeric(nn)) { - PEcAn.logger::logger.severe(paste0( - "Expected raster object to be numeric, but it has type `", - paste0(typeof(nn), collapse = " "), - "`" - )) + + # msg + PEcAn.logger::logger.info(paste0(year, " is being processed ", "for ensemble #", ens, " ")) + } + + # open the file + nc_data <- ncdf4::nc_open(ncfile) + # time stamp + + t <- ncdf4::ncvar_get(nc_data, "time") + tunits <- ncdf4::ncatt_get(nc_data, "time") + tustr <- strsplit(tunits$units, " ") + timestamp <- + as.POSIXct(t * 3600, tz = "UTC", origin = tustr[[1]][3]) + try(ncdf4::nc_close(nc_data)) + + + # set the vars + if (is.null(vars)) { + vars <- names(nc_data$var) + } + # for the variables extract the data + all.data.point <- vars %>% + purrr::set_names(vars) %>% + purrr::map_dfc(function(vname) { + if (verbose) { + PEcAn.logger::logger.info(paste0(" \t ", vname, "is being extracted ! ")) } - } - - # replacing the missing/filled values with NA - nn[nn == nc_data$var[[vname]]$missval] <- NA - # send out the extracted var as a new col - t(nn) - - }) - - #close the connection - - # send out as xts object - xts::xts(all.data.point, order.by = timestamp) - }) %>% - stats::setNames(paste0("ERA_ensemble_", ensemblesN)) - - #Merge mean and the speard - return(point.data) - - }) %>% - stats::setNames(years) - - - # The order of one.year.out is year and then Ens - Mainly because of the spead / I wanted to touch each file just once. - # This now changes the order to ens - year - point.data <- ensemblesN %>% - purrr::map(function(Ensn) { - rbind.xts <- do.call("::", list("xts", "rbind.xts")) - one.year.out %>% - purrr::map( ~ .x [[Ensn]]) %>% - do.call("rbind.xts", .) - }) - - - # Calling the met2CF inside extract bc in met process met2CF comes before extract ! - out <-met2CF.ERA5( - slat, - slon, - start_date, - end_date, - sitename=newsite, - outfolder, - point.data, - overwrite = FALSE, - verbose = verbose - ) - return(out) - - }, error = function(e) { - PEcAn.logger::logger.severe(paste0(conditionMessage(e))) - }) - - } \ No newline at end of file + + brick.tmp <- + raster::brick(ncfile, varname = vname, level = ens) + nn <- + raster::extract(brick.tmp, + sp::SpatialPoints(cbind(slon, slat)), + method = "simple" + ) + if (verbose) { + if (!is.numeric(nn)) { + PEcAn.logger::logger.severe(paste0( + "Expected raster object to be numeric, but it has type `", + paste0(typeof(nn), collapse = " "), + "`" + )) + } + } + + # replacing the missing/filled values with NA + nn[nn == nc_data$var[[vname]]$missval] <- NA + # send out the extracted var as a new col + t(nn) + }) + + # close the connection + + # send out as xts object + xts::xts(all.data.point, order.by = timestamp) + }) %>% + stats::setNames(paste0("ERA_ensemble_", ensemblesN)) + + # Merge mean and the speard + return(point.data) + }) %>% + stats::setNames(years) + + + # The order of one.year.out is year and then Ens - Mainly because of the spead / I wanted to touch each file just once. + # This now changes the order to ens - year + point.data <- ensemblesN %>% + purrr::map(function(Ensn) { + rbind.xts <- do.call("::", list("xts", "rbind.xts")) + one.year.out %>% + purrr::map(~ .x[[Ensn]]) %>% + do.call("rbind.xts", .) + }) + + + # Calling the met2CF inside extract bc in met process met2CF comes before extract ! + out <- met2CF.ERA5( + slat, + slon, + start_date, + end_date, + sitename = newsite, + outfolder, + point.data, + overwrite = FALSE, + verbose = verbose + ) + return(out) + }, + error = function(e) { + PEcAn.logger::logger.severe(paste0(conditionMessage(e))) + } + ) + } diff --git a/modules/data.atmosphere/man/ERA5_met_process.Rd b/modules/data.atmosphere/man/ERA5_met_process.Rd index 29710098dd0..61be241433f 100644 --- a/modules/data.atmosphere/man/ERA5_met_process.Rd +++ b/modules/data.atmosphere/man/ERA5_met_process.Rd @@ -18,7 +18,8 @@ ERA5_met_process(settings, in.path, out.path, write.db = FALSE, write = TRUE) \item{write}{if write the settings into pecan.xml file in the outdir of settings.} } \value{ -if write.db is True then return input IDs with physical paths; if write.db is False then return just physical paths of extracted ERA5 clim files. +if write.db is True then return input IDs with physical paths; + if write.db is False then return just physical paths of extracted ERA5 clim files. } \description{ Met Processes for ERA5 data diff --git a/modules/data.atmosphere/man/extract.nc.ERA5.Rd b/modules/data.atmosphere/man/extract.nc.ERA5.Rd index d377f131bd4..9b5d55217f3 100644 --- a/modules/data.atmosphere/man/extract.nc.ERA5.Rd +++ b/modules/data.atmosphere/man/extract.nc.ERA5.Rd @@ -58,9 +58,9 @@ For the list of variables check out the documentation at \url{ } \examples{ \dontrun{ -point.data <- ERA5_extract(sslat=40, slon=-120, years=c(1990:1995), vars=NULL) - - purrr::map(~xts::apply.daily(.x, mean)) +point.data <- ERA5_extract(sslat = 40, slon = -120, years = c(1990:1995), vars = NULL) +# point.data \%>\% +purrr::map(~ xts::apply.daily(.x, mean)) } } From 04139f8be62c7c6e9dbf3632a100a5a9a6837d56 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Fri, 3 Jan 2025 23:26:44 -0800 Subject: [PATCH 3/4] run styler on all files in models/sipnet --- models/sipnet/R/met2model.SIPNET.R | 277 +++++----- models/sipnet/R/model2netcdf.SIPNET.R | 471 +++++++++++------- models/sipnet/R/read_restart.SIPNET.R | 159 +++--- models/sipnet/R/sample.IC.SIPNET.R | 68 +-- models/sipnet/R/split_inputs.SIPNET.R | 46 +- models/sipnet/R/veg2model.SIPNET.R | 15 +- models/sipnet/R/write.configs.SIPNET.R | 361 +++++++------- models/sipnet/R/write_restart.SIPNET.R | 85 ++-- models/sipnet/inst/SIPNET.lyford.SDA.R | 28 +- models/sipnet/man/mergeNC.Rd | 14 +- models/sipnet/man/model2netcdf.SIPNET.Rd | 3 +- models/sipnet/tests/testthat.R | 2 +- models/sipnet/tests/testthat/test.met2model.R | 3 +- 13 files changed, 861 insertions(+), 671 deletions(-) diff --git a/models/sipnet/R/met2model.SIPNET.R b/models/sipnet/R/met2model.SIPNET.R index d1c57084318..352e9592555 100644 --- a/models/sipnet/R/met2model.SIPNET.R +++ b/models/sipnet/R/met2model.SIPNET.R @@ -20,21 +20,19 @@ ##' @author Luke Dramko, Michael Dietze, Alexey Shiklomanov, Rob Kooper met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date, overwrite = FALSE, verbose = FALSE, year.fragment = FALSE, ...) { - - - PEcAn.logger::logger.info("START met2model.SIPNET") start_date <- as.POSIXlt(start_date, tz = "UTC") end_date <- as.POSIXlt(end_date, tz = "UTC") if (year.fragment) { # in.prefix is not guaranteed to contain the file extension. escaped <- gsub("(\\W)", "\\\\\\1", in.prefix) # The file name may contain special characters that could mess up the regular expression. - matching_files <- grep(escaped, list.files(in.path), value=TRUE) + matching_files <- grep(escaped, list.files(in.path), value = TRUE) if (length(matching_files) == 0) { PEcAn.logger::logger.severe(paste0("No files found matching ", in.prefix, "; cannot process data.")) } - - # This function is supposed to process netcdf files, so we'll search for files with the extension .nc and use those first. - nc_file = grep("\\.nc$", matching_files) + + # This function is supposed to process netcdf files, + # so we'll search for files with the extension .nc and use those first. + nc_file <- grep("\\.nc$", matching_files) if (length(nc_file) > 0) { if (grepl("\\.nc$", in.prefix)) { out.file <- sub("\\.nc$", ".clim", in.prefix) @@ -49,36 +47,39 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date } } else { # Default behavior out.file <- paste(in.prefix, strptime(start_date, "%Y-%m-%d"), - strptime(end_date, "%Y-%m-%d"), - "clim", - sep = ".") + strptime(end_date, "%Y-%m-%d"), + "clim", + sep = "." + ) } - + out.file.full <- file.path(outfolder, out.file) - - results <- data.frame(file = out.file.full, - host = PEcAn.remote::fqdn(), - mimetype = "text/csv", - formatname = "Sipnet.climna", - startdate = start_date, - enddate = end_date, - dbfile.name = out.file, - stringsAsFactors = FALSE) + + results <- data.frame( + file = out.file.full, + host = PEcAn.remote::fqdn(), + mimetype = "text/csv", + formatname = "Sipnet.climna", + startdate = start_date, + enddate = end_date, + dbfile.name = out.file, + stringsAsFactors = FALSE + ) PEcAn.logger::logger.info("internal results") PEcAn.logger::logger.info(results) - + if (file.exists(out.file.full) && !overwrite) { PEcAn.logger::logger.debug("File '", out.file.full, "' already exists, skipping to next file.") return(invisible(results)) } - + ## check to see if the outfolder is defined, if not create directory for output if (!file.exists(outfolder)) { dir.create(outfolder) } - + out <- NULL - + # get start/end year since inputs are specified on year basis # only if year.fragment = FALSE start_year <- lubridate::year(start_date) @@ -88,67 +89,66 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date } else { end_year <- lubridate::year(end_date) } - + ## loop over files for (year in start_year:end_year) { - skip <- FALSE PEcAn.logger::logger.info(year) - + diy <- PEcAn.utils::days_in_year(year) - + if (!year.fragment) { # default behavior old.file <- file.path(in.path, paste(in.prefix, year, "nc", sep = ".")) } else { # Use the supplied file name old.file <- file.path(in.path, in.prefix) } - + if (file.exists(old.file)) { ## open netcdf - nc <- ncdf4::nc_open(old.file) - + nc <- ncdf4::nc_open(old.file) + ## convert time to seconds sec <- nc$dim$time$vals sec <- PEcAn.utils::ud_convert(sec, unlist(strsplit(nc$dim$time$units, " "))[1], "seconds") - - # Calculate the delta time. If using whole-year data, the appropriate length in seconds is - # fetched; otherwise, it is assumed that the length of time provided in the time dimension of - # the input file is correct. + + # Calculate the delta time. + # If using whole-year data, the appropriate length in seconds is + # fetched; otherwise, it is assumed that the length of time provided in + # the time dimension of the input file is correct. if (year.fragment) { - dt <- mean(diff(sec), na.rm=TRUE) - + dt <- mean(diff(sec), na.rm = TRUE) } else { dt <- PEcAn.utils::seconds_in_year(year) / length(sec) } tstep <- round(86400 / dt) dt <- 86400 / tstep - + ## extract variables lat <- ncdf4::ncvar_get(nc, "latitude") lon <- ncdf4::ncvar_get(nc, "longitude") - Tair <-ncdf4::ncvar_get(nc, "air_temperature") ## in Kelvin + Tair <- ncdf4::ncvar_get(nc, "air_temperature") ## in Kelvin Tair_C <- PEcAn.utils::ud_convert(Tair, "K", "degC") - Qair <-ncdf4::ncvar_get(nc, "specific_humidity") #humidity (kg/kg) + Qair <- ncdf4::ncvar_get(nc, "specific_humidity") # humidity (kg/kg) ws <- try(ncdf4::ncvar_get(nc, "wind_speed")) if (!is.numeric(ws)) { U <- ncdf4::ncvar_get(nc, "eastward_wind") V <- ncdf4::ncvar_get(nc, "northward_wind") - ws <- sqrt(U ^ 2 + V ^ 2) + ws <- sqrt(U^2 + V^2) PEcAn.logger::logger.info("wind_speed absent; calculated from eastward_wind and northward_wind") } - + Rain <- ncdf4::ncvar_get(nc, "precipitation_flux") - - press <- ncdf4::ncvar_get(nc,'air_pressure') ## in pascal - SW <- ncdf4::ncvar_get(nc, "surface_downwelling_shortwave_flux_in_air") ## in W/m2 - - PAR <- try(ncdf4::ncvar_get(nc, "surface_downwelling_photosynthetic_photon_flux_in_air")) ## in umol/m2/s + press <- ncdf4::ncvar_get(nc, "air_pressure") ## in pascal + + SW <- ncdf4::ncvar_get(nc, "surface_downwelling_shortwave_flux_in_air") ## in W/m2 + + PAR <- try(ncdf4::ncvar_get(nc, "surface_downwelling_photosynthetic_photon_flux_in_air")) ## in umol/m2/s if (!is.numeric(PAR)) { PAR <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::sw2ppfd(SW), "umol ", "mol") PEcAn.logger::logger.info("surface_downwelling_photosynthetic_photon_flux_in_air absent; PAR set to SW * 0.45") } - + soilT <- try(ncdf4::ncvar_get(nc, "soil_temperature")) if (!is.numeric(soilT)) { # approximation borrowed from SIPNET CRUNCEP preprocessing's tsoil.py @@ -161,25 +161,24 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date } else { soilT <- PEcAn.utils::ud_convert(soilT, "K", "degC") } - - SVP <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::get.es(Tair_C), "millibar", "Pa") ## Saturation vapor pressure - VPD <- try(ncdf4::ncvar_get(nc, "water_vapor_saturation_deficit")) ## in Pa - if (!is.numeric(VPD)) { - VPD <- SVP * (1 - PEcAn.data.atmosphere::qair2rh(Qair, Tair_C, press = press/100)) + SVP <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::get.es(Tair_C), "millibar", "Pa") ## Saturation vapor pressure + VPD <- try(ncdf4::ncvar_get(nc, "water_vapor_saturation_deficit")) ## in Pa + if (!is.numeric(VPD)) { + VPD <- SVP * (1 - PEcAn.data.atmosphere::qair2rh(Qair, Tair_C, press = press / 100)) PEcAn.logger::logger.info("water_vapor_saturation_deficit absent; VPD calculated from Qair, Tair, and SVP (saturation vapor pressure) ") } e_a <- SVP - VPD VPDsoil <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::get.es(soilT), "millibar", "Pa") * - (1 - PEcAn.data.atmosphere::qair2rh(Qair, soilT, press/100)) - + (1 - PEcAn.data.atmosphere::qair2rh(Qair, soilT, press / 100)) + ncdf4::nc_close(nc) } else { PEcAn.logger::logger.info("Skipping to next year") next } - + ## build time variables (year, month, day of year) nyr <- floor(length(sec) / 86400 / 365 * dt) yr <- NULL @@ -216,103 +215,135 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date } yr[rng] <- rep(y + 1, length(rng)) doy[rng] <- rep(1:366, each = 86400 / dt)[1:length(rng)] - hr[rng] <- rep(seq(0, length = 86400 / dt, by = dt/86400 * 24), 366)[1:length(rng)] + hr[rng] <- rep(seq(0, length = 86400 / dt, by = dt / 86400 * 24), 366)[1:length(rng)] } if (skip) { PEcAn.logger::logger.info("Skipping to next year") next } - - ## 0 YEAR DAY HOUR TIMESTEP AirT SoilT PAR PRECIP VPD VPD_Soil AirVP(e_a) WIND SoilM build data - ## matrix + + ## build data matrix + ## 0 YEAR DAY HOUR TIMESTEP AirT SoilT PAR PRECIP VPD VPD_Soil AirVP(e_a) WIND SoilM n <- length(Tair) - tmp <- cbind(rep(0, n), - yr[1:n], - doy[1:n], - hr[1:n], - rep(dt / 86400, n), - Tair_C, - soilT, - PAR * dt, # mol/m2/hr - Rain * dt, # converts from mm/s to mm - VPD, - VPDsoil, - e_a, - ws, # wind - rep(0.6, n)) # put soil water at a constant. Don't use, set SIPNET to MODEL_WATER = 1 - + tmp <- cbind( + rep(0, n), + yr[1:n], + doy[1:n], + hr[1:n], + rep(dt / 86400, n), + Tair_C, + soilT, + PAR * dt, # mol/m2/hr + Rain * dt, # converts from mm/s to mm + VPD, + VPDsoil, + e_a, + ws, # wind + rep(0.6, n) # put soil water at a constant. Don't use, set SIPNET to MODEL_WATER = 1 + ) + ## quick error check, sometimes get a NA in the last hr hr.na <- which(is.na(tmp[, 4])) if (length(hr.na) > 0) { - tmp[hr.na, 4] <- tmp[hr.na - 1, 4] + dt/86400 * 24 + tmp[hr.na, 4] <- tmp[hr.na - 1, 4] + dt / 86400 * 24 } - - ## filter out days not included in start or end date if not a year fragment. (This procedure would be nonsensible for a year - ## fragment, as it would filter out all of the days.) - if(year == start_year && !year.fragment){ - extra.days <- length(as.Date(paste0(start_year, "-01-01")):as.Date(start_date)) #extra days length includes the start date - if (extra.days > 1){ + + ## filter out days not included in start or end date if not a year fragment. + ## (This procedure would be nonsensible for a year fragment, + ## as it would filter out all of the days.) + if (year == start_year && !year.fragment) { + # extra days length includes the start date + extra.days <- length( + as.Date(paste0(start_year, "-01-01")):as.Date(start_date) + ) + if (extra.days > 1) { PEcAn.logger::logger.info("Subsetting SIPNET met to match start date") - start.row <- ((extra.days - 1) * 86400 / dt) + 1 #subtract to include start.date, add to exclude last half hour of day before - tmp <- tmp[start.row:nrow(tmp),] + # subtract to include start.date, add to exclude last half hour of day before + start.row <- ((extra.days - 1) * 86400 / dt) + 1 + tmp <- tmp[start.row:nrow(tmp), ] } } - if (year == end_year && !year.fragment){ - if(year == start_year){ - extra.days <- length(as.Date(start_date):as.Date(end_date)) - if (extra.days > 1){ + if (year == end_year && !year.fragment) { + if (year == start_year) { + extra.days <- length(as.Date(start_date):as.Date(end_date)) + if (extra.days > 1) { PEcAn.logger::logger.info("Subsetting SIPNET met to match end date") - end.row <- extra.days * 86400 / dt #subtract to include end.date - tmp <- tmp[1:end.row,] + end.row <- extra.days * 86400 / dt # subtract to include end.date + tmp <- tmp[1:end.row, ] } - } else{ - extra.days <- length(as.Date(end_date):as.Date(paste0(end_year, "-12-31"))) #extra days length includes the end date - if (extra.days > 1){ + } else { + # extra days length includes the end date + extra.days <- length( + as.Date(end_date):as.Date(paste0(end_year, "-12-31")) + ) + if (extra.days > 1) { PEcAn.logger::logger.info("Subsetting SIPNET met to match end date") - end.row <- nrow(tmp) - ((extra.days - 1) * 86400 / dt) #subtract to include end.date - tmp <- tmp[1:end.row,] + end.row <- nrow(tmp) - ((extra.days - 1) * 86400 / dt) # subtract to include end.date + tmp <- tmp[1:end.row, ] } } } - - if(year.fragment){ #gets correct DOY for fragmented years - - doy.start.index <- which(doy == lubridate::yday(start_date)) #which part of full doy set matches the start and end date - doy.end.index <- which(doy == lubridate::yday(end_date)) - #need to use the start and end time to figure out how many time steps to include in the doy subset - doy.start <- doy.start.index[ifelse(lubridate::hour(start_date) == 0, 1, lubridate::hour(start_date) / (24 / (86400 / dt)))] - doy.end <- doy.end.index[ifelse(lubridate::hour(end_date) == 0, 1, lubridate::hour(end_date) / (24 / (86400 / dt)))] - #check to see if doy matches with downloaded data dims, if not last time is removed - if(length(doy) != n){d<-doy[doy.start:(doy.end - 1)] }else{d<-(doy[doy.start:(doy.end)]) } - - if(year.fragment){ #gets correct DOY for fragmented years using start date, time since start date and end date - doy.seq <- as.Date(seq(from = start_date + sec[1], to = end_date, length.out = length(sec))) - doy <- as.numeric(strftime(doy.seq, format = "%j")) #starts with 1 on 1-01 - #doy.start <- length(as.Date(paste0(start_year, "-01-01")):as.Date(start_date)) * (86400 / dt) + 1 #subtract to include start.date, add to exclude last half hour of day before - #doy.end <- length(as.Date(paste0(start_year, "-01-01")):as.Date(end_date)) * (86400 / dt) - #doy <- doy[doy.start:doy.end] - year <- as.numeric(strftime(doy.seq, format = "%Y")) - tmp[,3] <- doy - tmp[,2] <- year + + if (year.fragment) { # gets correct DOY for fragmented years + + doy.start.index <- which(doy == lubridate::yday(start_date)) # which part of full doy set matches the start and end date + doy.end.index <- which(doy == lubridate::yday(end_date)) + # need to use the start and end time to figure out how many time steps to include in the doy subset + doy.start <- doy.start.index[ifelse( + lubridate::hour(start_date) == 0, + 1, + lubridate::hour(start_date) / (24 / (86400 / dt)) + )] + doy.end <- doy.end.index[ifelse( + lubridate::hour(end_date) == 0, + 1, + lubridate::hour(end_date) / (24 / (86400 / dt)) + )] + # check to see if doy matches with downloaded data dims, + # if not last time is removed + if (length(doy) != n) { + d <- doy[doy.start:(doy.end - 1)] + } else { + d <- (doy[doy.start:(doy.end)]) + } + + if (year.fragment) { # gets correct DOY for fragmented years using start date, time since start date and end date + doy.seq <- as.Date(seq( + from = start_date + sec[1], + to = end_date, + length.out = length(sec) + )) + doy <- as.numeric(strftime(doy.seq, format = "%j")) # starts with 1 on 1-01 + # doy.start <- length(as.Date(paste0(start_year, "-01-01")):as.Date(start_date)) * (86400 / dt) + 1 #subtract to include start.date, add to exclude last half hour of day before + # doy.end <- length(as.Date(paste0(start_year, "-01-01")):as.Date(end_date)) * (86400 / dt) + # doy <- doy[doy.start:doy.end] + year <- as.numeric(strftime(doy.seq, format = "%Y")) + tmp[, 3] <- doy + tmp[, 2] <- year + } } - } - + if (is.null(out)) { out <- tmp } else { out <- rbind(out, tmp) } - - } ## end loop over years - + } ## end loop over years + if (!is.null(out)) { - ## write output - utils::write.table(out, out.file.full, quote = FALSE, sep = "\t", row.names = FALSE, col.names = FALSE) + utils::write.table( + out, + out.file.full, + quote = FALSE, + sep = "\t", + row.names = FALSE, + col.names = FALSE + ) return(invisible(results)) } else { PEcAn.logger::logger.info("NO MET TO OUTPUT") return(invisible(NULL)) } -} # met2model.SIPNET \ No newline at end of file +} # met2model.SIPNET diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R index 438a5c43de3..ecf6c8b489f 100644 --- a/models/sipnet/R/model2netcdf.SIPNET.R +++ b/models/sipnet/R/model2netcdf.SIPNET.R @@ -1,60 +1,64 @@ #' Merge multiple NetCDF files into one -#' -#' @param files \code{character}. List of filepaths, which should lead to NetCDF files. +#' +#' @param files \code{character}. List of filepaths, which should lead to +#' NetCDF files. #' @param outfile \code{character}. Output filename of the merged data. #' @return A NetCDF file containing all of the merged data. #' @examples #' \dontrun{ -#' files <- list.files(paste0(system.file(package="processNC"), "/extdata"), -#' pattern="tas.*\\.nc", full.names=TRUE) -#' temp <- tempfile(fileext=".nc") -#' mergeNC(files=files, outfile=temp) -#' terra::rast(temp) +#' files <- list.files(paste0(system.file(package = "processNC"), "/extdata"), +#' pattern = "tas.*\\.nc", full.names = TRUE +#' ) +#' temp <- tempfile(fileext = ".nc") +#' mergeNC(files = files, outfile = temp) +#' terra::rast(temp) #' } #' @export mergeNC #' @name mergeNC #' @source https://github.com/RS-eco/processNC/blob/main/R/mergeNC.R -mergeNC <- function( - ##title<< Aggregate data in netCDF files - files ##<< character vector: names of the files to merge - , outfile ##<< character: path to save the results files to. -) - ##description<< - ## This function aggregates time periods in netCDF files. Basically it is just a - ## wrapper around the respective cdo function. -{ - ##test input - #if (system("cdo -V")==0) +## description<< +## This function aggregates time periods in netCDF files. Basically it is just a +## wrapper around the respective cdo function. +mergeNC <- function(files, outfile) { + ## test input + # if (system("cdo -V")==0) # stop('cdo not found. Please install it.') - + ## supply cdo command - cdoCmd <- paste('cdo -cat', paste(files, collapse=" "), outfile, sep=' ') - - ##run command + cdoCmd <- paste("cdo -cat", paste(files, collapse = " "), outfile, sep = " ") + + ## run command system(cdoCmd) - cat(paste('Created file ', outfile, '.\n', sep = '')) - - ## character string: name of the file created. + cat(paste("Created file ", outfile, ".\n", sep = "")) + + ## character string: name of the file created. invisible(outfile) } -#--------------------------------------------------------------------------------------------------# + + + + + +#------------------------------------------------------------------------------# ##' ##' Convert SIPNET DOY to datetime -##' +##' ##' @param sipnet_tval vector of SIPNET DOY values ##' @param base_year base year to calculate datetime from DOY -##' @param base_month reference month for converting from DOY to datetime +##' @param base_month reference month for converting from DOY to datetime ##' @param force_cf force output to follow CF convention. Default FALSE ##' ##' @export ##' ##' @author Alexey Shiklomanov, Shawn Serbin -##' +##' sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1, force_cf = FALSE) { - base_date <- ISOdatetime(base_year, base_month, 1, - 0, 0, 0, "UTC") + base_date <- ISOdatetime( + base_year, base_month, 1, + 0, 0, 0, "UTC" + ) base_date_str <- strftime(base_date, "%F %T %z", tz = "UTC") if (force_cf) { is_cf <- TRUE @@ -63,17 +67,22 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1, # Is CF if first time step is zero is_cf <- sipnet_tval[[1]] == 0 } - + if (is_cf) { cfval <- sipnet_tval } else { cfval <- sipnet_tval - 1 } - + PEcAn.utils::cf2datetime(cfval, paste("days since", base_date_str)) } -#--------------------------------------------------------------------------------------------------# + + + + + +#------------------------------------------------------------------------------# ##' Convert SIPNET output to netCDF ##' ##' Converts all output contained in a folder to netCDF. @@ -85,157 +94,233 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1, ##' @param end_date End time of the simulation ##' @param revision model revision ##' @param overwrite Flag for overwriting nc files or not -##' @param conflict Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next. +##' @param conflict Flag for dealing with conflicted nc files, +##' if T we then will merge those, if F we will jump to the next. ##' @param prefix prefix to read the output files ##' @param delete.raw logical: remove sipnet.out files after converting? ##' ##' @export ##' @author Shawn Serbin, Michael Dietze -model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out", +model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, + delete.raw = FALSE, revision, + prefix = "sipnet.out", overwrite = FALSE, conflict = FALSE) { - ### Read in model output in SIPNET format sipnet_out_file <- file.path(outdir, prefix) - sipnet_output <- utils::read.table(sipnet_out_file, header = T, skip = 1, sep = "") - #sipnet_output_dims <- dim(sipnet_output) + sipnet_output <- utils::read.table( + sipnet_out_file, + header = TRUE, + skip = 1, + sep = "" + ) + # sipnet_output_dims <- dim(sipnet_output) ### Determine number of years and output timestep - #start.day <- sipnet_output$day[1] + # start.day <- sipnet_output$day[1] num_years <- length(unique(sipnet_output$year)) simulation_years <- unique(sipnet_output$year) - + # get all years that we want data from year_seq <- seq(lubridate::year(start_date), lubridate::year(end_date)) # check that specified years and output years match if (!all(year_seq %in% simulation_years)) { - PEcAn.logger::logger.severe("Years selected for model run and SIPNET output years do not match ") + PEcAn.logger::logger.severe( + "Years selected for model run and SIPNET output years do not match " + ) } # get number of model timesteps per day - # outday is the number of time steps in a day - for example 6 hours would have out_day of 4 + # outday is the number of time steps in a day - for example 6 hours would have + # out_day of 4 + + out_day <- sum( + sipnet_output$year == simulation_years[1] & + sipnet_output$day == unique(sipnet_output$day)[1], + na.rm = TRUE + ) # switched to day 2 in case first day is partial - out_day <- sum( - sipnet_output$year == simulation_years[1] & - sipnet_output$day == unique(sipnet_output$day)[1], - na.rm = TRUE - ) # switched to day 2 in case first day is partial - timestep.s <- 86400 / out_day - - + + ### Loop over years in SIPNET output to create separate netCDF outputs for (y in year_seq) { - #initialize the conflicted as FALSE + # initialize the conflicted as FALSE conflicted <- FALSE - conflict <- TRUE #conflict is set to TRUE to enable the rename of yearly nc file for merging SDA results with sub-annual data - #if we have conflicts on this file. - if (file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & overwrite == FALSE & conflict == FALSE) { + # conflict is set to TRUE to enable the rename of yearly nc file + # for merging SDA results with sub-annual data + conflict <- TRUE + # if we have conflicts on this file. + if (file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & + overwrite == FALSE & + conflict == FALSE) { next - }else if(file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & conflict){ + } else if (file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & + conflict) { conflicted <- TRUE - file.rename(file.path(outdir, paste(y, "nc", sep = ".")), file.path(outdir, "previous.nc")) + file.rename( + file.path(outdir, paste(y, "nc", sep = ".")), + file.path(outdir, "previous.nc") + ) } - print(paste("---- Processing year: ", y)) # turn on for debugging + print(paste("---- Processing year: ", y)) # turn on for debugging ## Subset data for processing sub.sipnet.output <- subset(sipnet_output, year == y) sub.sipnet.output.dims <- dim(sub.sipnet.output) dayfrac <- 1 / out_day - step <- utils::head(seq(0, 1, by = dayfrac), -1) ## probably dont want to use - ## hard-coded "step" because - ## leap years may not contain - ## all "steps", or - ## if model run doesnt start - ## at 00:00:00 - - # try to determine if DOY is CF compliant (i.e. 0 based index) or not (1 base index) - pecan_start_doy <- PEcAn.utils::datetime2cf(start_date, paste0("days since ",lubridate::year(start_date),"-01-01"), - tz = "UTC") + + # probably dont want to use hard-coded "step" because leap years may not + # contain all "steps", or if model run doesnt start at 00:00:00 + step <- utils::head(seq(0, 1, by = dayfrac), -1) + + # try to determine if DOY is CF compliant (i.e. 0 based index) + # or not (1 base index) + pecan_start_doy <- PEcAn.utils::datetime2cf( + start_date, + paste0("days since ", lubridate::year(start_date), "-01-01"), + tz = "UTC" + ) tvals <- sub.sipnet.output[["day"]] + sub.sipnet.output[["time"]] / 24 - if (sub.sipnet.output[["day"]][1]-pecan_start_doy==1) { + if (sub.sipnet.output[["day"]][1] - pecan_start_doy == 1) { sub_dates <- sipnet2datetime(tvals, y, force_cf = FALSE) } else { sub_dates <- sipnet2datetime(tvals, y, force_cf = TRUE) } - sub_dates_cf <- PEcAn.utils::datetime2cf(sub_dates, paste0("days since ",paste0(y,"-01-01"))) + sub_dates_cf <- PEcAn.utils::datetime2cf( + sub_dates, + paste0("days since ", paste0(y, "-01-01")) + ) # create netCDF time.bounds variable - bounds <- array(data=NA, dim=c(length(sub_dates_cf),2)) - bounds[,1] <- sub_dates_cf - bounds[,2] <- bounds[,1]+dayfrac + bounds <- array(data = NA, dim = c(length(sub_dates_cf), 2)) + bounds[, 1] <- sub_dates_cf + bounds[, 2] <- bounds[, 1] + dayfrac # create time bounds for each timestep in t, t+1; t+1, t+2... format - bounds <- round(bounds,4) - + bounds <- round(bounds, 4) + ## Setup outputs for netCDF file in appropriate units - output <- list( - "GPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s, # GPP in kgC/m2/s - "NPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s - ((sub.sipnet.output$rAboveground * - 0.001) / timestep.s + (sub.sipnet.output$rRoot * 0.001) / timestep.s), # NPP in kgC/m2/s. Post SIPNET calculation - "TotalResp" = (sub.sipnet.output$rtot * 0.001) / timestep.s, # Total Respiration in kgC/m2/s - "AutoResp" = (sub.sipnet.output$rAboveground * 0.001) / timestep.s + (sub.sipnet.output$rRoot * - 0.001) / timestep.s, # Autotrophic Respiration in kgC/m2/s - "HeteroResp" = ((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot) * 0.001) / timestep.s, # Heterotrophic Respiration in kgC/m2/s - "SoilResp" = (sub.sipnet.output$rSoil * 0.001) / timestep.s, # Soil Respiration in kgC/m2/s - "NEE" = (sub.sipnet.output$nee * 0.001) / timestep.s, # NEE in kgC/m2/s - "AbvGrndWood" = (sub.sipnet.output$plantWoodC * 0.001), # Above ground wood kgC/m2 - "leaf_carbon_content" = (sub.sipnet.output$plantLeafC * 0.001), # Leaf C kgC/m2 - "TotLivBiom" = (sub.sipnet.output$plantWoodC * 0.001) + (sub.sipnet.output$plantLeafC * 0.001) + - (sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001, # Total living C kgC/m2 - "TotSoilCarb" = (sub.sipnet.output$soil * 0.001) + (sub.sipnet.output$litter * 0.001) # Total soil C kgC/m2 + output <- list( + # GPP in kgC/m2/s + "GPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s, + # NPP in kgC/m2/s. Post SIPNET calculation + "NPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s - ( + (sub.sipnet.output$rAboveground * 0.001) / timestep.s + + (sub.sipnet.output$rRoot * 0.001) / timestep.s + ), + # Total Respiration in kgC/m2/s + "TotalResp" = (sub.sipnet.output$rtot * 0.001) / timestep.s, + # Autotrophic Respiration in kgC/m2/s + "AutoResp" = (sub.sipnet.output$rAboveground * 0.001) / timestep.s + + (sub.sipnet.output$rRoot * 0.001) / timestep.s, + # Heterotrophic Respiration in kgC/m2/s + "HeteroResp" = ((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot) * + 0.001) / timestep.s, + # Soil Respiration in kgC/m2/s + "SoilResp" = (sub.sipnet.output$rSoil * 0.001) / timestep.s, + # NEE in kgC/m2/s + "NEE" = (sub.sipnet.output$nee * 0.001) / timestep.s, + # Above ground wood kgC/m2 + "AbvGrndWood" = (sub.sipnet.output$plantWoodC * 0.001), + # Leaf C kgC/m2 + "leaf_carbon_content" = (sub.sipnet.output$plantLeafC * 0.001), + # Total living C kgC/m2 + "TotLivBiom" = (sub.sipnet.output$plantWoodC * 0.001) + + (sub.sipnet.output$plantLeafC * 0.001) + + (sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001, + # Total soil C kgC/m2 + "TotSoilCarb" = (sub.sipnet.output$soil * 0.001) + + (sub.sipnet.output$litter * 0.001) ) + # Qle W/m2 if (revision == "unk") { - ## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : *** - ## *** it says "npp" in the header (written by L774) but the values being written are trackers.evapotranspiration (L806) *** - ## evapotranspiration in SIPNET is cm^3 water per cm^2 of area, to convert it to latent heat units W/m2 multiply with : - ## 0.01 (cm2m) * 1000 (water density, kg m-3) * latent heat of vaporization (J kg-1) - ## latent heat of vaporization is not constant and it varies slightly with temperature, get.lv() returns 2.5e6 J kg-1 by default - output[["Qle"]] <- (sub.sipnet.output$npp * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s # Qle W/m2 + ## *** NOTE : npp in the sipnet output file is actually evapotranspiration + ## this is due to a bug in sipnet.c : + ## *** it says "npp" in the header (written by L774) but the values being + ## written are trackers.evapotranspiration (L806) + ## evapotranspiration in SIPNET is cm^3 water per cm^2 of area. + ## To convert it to latent heat units W/m2 multiply with : + ## 0.01 (cm2m) + ## * 1000 (water density, kg m-3) + ## * latent heat of vaporization (J kg-1) + ## latent heat of vaporization is not constant and it varies slightly with + ## temperature. get.lv() returns 2.5e6 J kg-1 by default + output[["Qle"]] <- (sub.sipnet.output$npp * 10 * + PEcAn.data.atmosphere::get.lv()) / timestep.s } else { - output[["Qle"]] <- (sub.sipnet.output$evapotranspiration * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s # Qle W/m2 + output[["Qle"]] <- (sub.sipnet.output$evapotranspiration * 10 * + PEcAn.data.atmosphere::get.lv()) / timestep.s } - output[["Transp"]] <- (sub.sipnet.output$fluxestranspiration * 10) / timestep.s # Transpiration kgW/m2/s - output[["SoilMoist"]] <- (sub.sipnet.output$soilWater * 10) # Soil moisture kgW/m2 - output[["SoilMoistFrac"]] <- (sub.sipnet.output$soilWetnessFrac) # Fractional soil wetness - output[["SWE"]] <- (sub.sipnet.output$snow * 10) # SWE - output[["litter_carbon_content"]] <- sub.sipnet.output$litter * 0.001 ## litter kgC/m2 - output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # Litter water kgW/m2 - #calculate LAI for standard output - param <- utils::read.table(file.path(gsub(pattern = "/out/", - replacement = "/run/", x = outdir), - "sipnet.param"), stringsAsFactors = FALSE) + # Transpiration kgW/m2/s + output[["Transp"]] <- (sub.sipnet.output$fluxestranspiration * 10) / + timestep.s + # Soil moisture kgW/m2 + output[["SoilMoist"]] <- (sub.sipnet.output$soilWater * 10) + # Fractional soil wetness + output[["SoilMoistFrac"]] <- (sub.sipnet.output$soilWetnessFrac) + # Snow water equivalent + output[["SWE"]] <- (sub.sipnet.output$snow * 10) + # litter kgC/m2 + output[["litter_carbon_content"]] <- sub.sipnet.output$litter * 0.001 + # Litter water kgW/m2 + output[["litter_mass_content_of_water"]] <- ( + sub.sipnet.output$litterWater * 10 + ) + # calculate LAI for standard output + param <- utils::read.table(file.path( + gsub( + pattern = "/out/", + replacement = "/run/", x = outdir + ), + "sipnet.param" + ), stringsAsFactors = FALSE) + # LAI id <- which(param[, 1] == "leafCSpWt") leafC <- 0.48 - SLA <- 1000 / param[id, 2] #SLA, m2/kgC - output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA # LAI - output[["fine_root_carbon_content"]] <- sub.sipnet.output$fineRootC * 0.001 ## fine_root_carbon_content kgC/m2 - output[["coarse_root_carbon_content"]] <- sub.sipnet.output$coarseRootC * 0.001 ## coarse_root_carbon_content kgC/m2 - output[["GWBI"]] <- (sub.sipnet.output$woodCreation * 0.001) / 86400 ## kgC/m2/s - this is daily in SIPNET - output[["AGB"]] <- (sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC) * 0.001 # Total aboveground biomass kgC/m2 - output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2])) - + SLA <- 1000 / param[id, 2] # SLA, m2/kgC + output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA + # fine_root_carbon_content kgC/m2 + output[["fine_root_carbon_content"]] <- sub.sipnet.output$fineRootC * 0.001 + # coarse_root_carbon_content kgC/m2 + output[["coarse_root_carbon_content"]] <- sub.sipnet.output$coarseRootC * 0.001 + # kgC/m2/s - this is daily in SIPNET + output[["GWBI"]] <- (sub.sipnet.output$woodCreation * 0.001) / 86400 + # Total aboveground biomass kgC/m2 + output[["AGB"]] <- (sub.sipnet.output$plantWoodC + + sub.sipnet.output$plantLeafC) * 0.001 + output[["time_bounds"]] <- c(rbind(bounds[, 1], bounds[, 2])) + # ******************** Declare netCDF variables ********************# - t <- ncdf4::ncdim_def(name = "time", - longname = "time", - units = paste0("days since ", y, "-01-01 00:00:00"), - vals = sub_dates_cf, - calendar = "standard", - unlim = TRUE) - lat <- ncdf4::ncdim_def("lat", "degrees_north", vals = as.numeric(sitelat), - longname = "station_latitude") - lon <- ncdf4::ncdim_def("lon", "degrees_east", vals = as.numeric(sitelon), - longname = "station_longitude") + t <- ncdf4::ncdim_def( + name = "time", + longname = "time", + units = paste0("days since ", y, "-01-01 00:00:00"), + vals = sub_dates_cf, + calendar = "standard", + unlim = TRUE + ) + lat <- ncdf4::ncdim_def("lat", "degrees_north", + vals = as.numeric(sitelat), + longname = "station_latitude" + ) + lon <- ncdf4::ncdim_def("lon", "degrees_east", + vals = as.numeric(sitelon), + longname = "station_longitude" + ) dims <- list(lon = lon, lat = lat, time = t) - time_interval <- ncdf4::ncdim_def(name = "hist_interval", - longname="history time interval endpoint dimensions", - vals = 1:2, units="") - + time_interval <- ncdf4::ncdim_def( + name = "hist_interval", + longname = "history time interval endpoint dimensions", + vals = 1:2, units = "" + ) + ## ***** Need to dynamically update the UTC offset here ***** for (i in seq_along(output)) { - if (length(output[[i]]) == 0) + if (length(output[[i]]) == 0) { output[[i]] <- rep(-999, length(t$vals)) + } } # ******************** Declare netCDF variables ********************# @@ -246,11 +331,17 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, "TotalResp" = PEcAn.utils::to_ncvar("TotalResp", dims), "AutoResp" = PEcAn.utils::to_ncvar("AutoResp", dims), "HeteroResp" = PEcAn.utils::to_ncvar("HeteroResp", dims), - "SoilResp" = ncdf4::ncvar_def("SoilResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, - longname = "Soil Respiration"), #need to figure out standard variable for this output + # need to figure out standard variable for this output + "SoilResp" = ncdf4::ncvar_def("SoilResp", + units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, + longname = "Soil Respiration" + ), "NEE" = PEcAn.utils::to_ncvar("NEE", dims), "AbvGrndWood" = PEcAn.utils::to_ncvar("AbvGrndWood", dims), - "leaf_carbon_content" = PEcAn.utils::to_ncvar("leaf_carbon_content", dims), + "leaf_carbon_content" = PEcAn.utils::to_ncvar( + "leaf_carbon_content", + dims + ), "TotLivBiom" = PEcAn.utils::to_ncvar("TotLivBiom", dims), "TotSoilCarb" = PEcAn.utils::to_ncvar("TotSoilCarb", dims), "Qle" = PEcAn.utils::to_ncvar("Qle", dims), @@ -258,64 +349,106 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, "SoilMoist" = PEcAn.utils::to_ncvar("SoilMoist", dims), "SoilMoistFrac" = PEcAn.utils::to_ncvar("SoilMoistFrac", dims), "SWE" = PEcAn.utils::to_ncvar("SWE", dims), - "litter_carbon_content" = PEcAn.utils::to_ncvar("litter_carbon_content", dims), - "litter_mass_content_of_water" = PEcAn.utils::to_ncvar("litter_mass_content_of_water", dims), + "litter_carbon_content" = PEcAn.utils::to_ncvar( + "litter_carbon_content", + dims + ), + "litter_mass_content_of_water" = PEcAn.utils::to_ncvar( + "litter_mass_content_of_water", + dims + ), "LAI" = PEcAn.utils::to_ncvar("LAI", dims), - "fine_root_carbon_content" = PEcAn.utils::to_ncvar("fine_root_carbon_content", dims), - "coarse_root_carbon_content" = PEcAn.utils::to_ncvar("coarse_root_carbon_content", dims), - "GWBI" = ncdf4::ncvar_def("GWBI", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, - longname = "Gross Woody Biomass Increment"), - "AGB" = ncdf4::ncvar_def("AGB", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, - longname = "Total aboveground biomass"), - "time_bounds" = ncdf4::ncvar_def(name="time_bounds", units='', - longname = "history time interval endpoints", dim=list(time_interval,time = t), - prec = "double") + "fine_root_carbon_content" = PEcAn.utils::to_ncvar( + "fine_root_carbon_content", + dims + ), + "coarse_root_carbon_content" = PEcAn.utils::to_ncvar( + "coarse_root_carbon_content", + dims + ), + "GWBI" = ncdf4::ncvar_def("GWBI", + units = "kg C m-2", dim = list(lon, lat, t), missval = -999, + longname = "Gross Woody Biomass Increment" + ), + "AGB" = ncdf4::ncvar_def("AGB", + units = "kg C m-2", dim = list(lon, lat, t), missval = -999, + longname = "Total aboveground biomass" + ), + "time_bounds" = ncdf4::ncvar_def( + name = "time_bounds", units = "", + longname = "history time interval endpoints", + dim = list(time_interval, time = t), + prec = "double" + ) ) - - # ******************** Create netCDF and output variables ********************# + + # ******************** Create netCDF and output variables *****************# ### Output netCDF data - if(conflicted & conflict){ - nc <- ncdf4::nc_create(file.path(outdir, paste("current", "nc", sep = ".")), nc_var) - ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec=NA) + if (conflicted & conflict) { + nc <- ncdf4::nc_create( + file.path(outdir, paste("current", "nc", sep = ".")), + nc_var + ) + ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec = NA) varfile <- file(file.path(outdir, paste(y, "nc", "var", sep = ".")), "w") for (key in names(nc_var)) { ncdf4::ncvar_put(nc, nc_var[[key]], output[[key]]) - cat(paste(nc_var[[key]]$name, nc_var[[key]]$longname), file = varfile, sep = "\n") + cat( + paste(nc_var[[key]]$name, nc_var[[key]]$longname), + file = varfile, + sep = "\n" + ) } close(varfile) ncdf4::nc_close(nc) - - #merge nc files of the same year together to enable the assimilation of sub-annual data - if(file.exists(file.path(outdir, "previous.nc"))){ - files <- c(file.path(outdir, "previous.nc"), file.path(outdir, "current.nc")) - }else{ + + # merge nc files of the same year together to enable the assimilation of + # sub-annual data + if (file.exists(file.path(outdir, "previous.nc"))) { + files <- c( + file.path(outdir, "previous.nc"), + file.path(outdir, "current.nc") + ) + } else { files <- file.path(outdir, "current.nc") } - mergeNC(files = files, outfile = file.path(outdir, paste(y, "nc", sep = "."))) - #The command "cdo" in mergeNC will automatically rename "time_bounds" to "time_bnds". However, "time_bounds" is used - #in read_restart codes later. So we need to read the new NetCDF file and convert the variable name back. - nc<- ncdf4::nc_open(file.path(outdir, paste(y, "nc", sep = ".")),write=TRUE) - nc<-ncdf4::ncvar_rename(nc,"time_bnds","time_bounds") - ncdf4::ncatt_put(nc, "time", "bounds","time_bounds", prec=NA) + mergeNC( + files = files, + outfile = file.path(outdir, paste(y, "nc", sep = ".")) + ) + # The command "cdo" in mergeNC will automatically rename "time_bounds" to + # "time_bnds". However, "time_bounds" is used in read_restart codes later. + # So we need to read the new NetCDF file and convert the variable name back. + nc <- ncdf4::nc_open( + file.path(outdir, paste(y, "nc", sep = ".")), + write = TRUE + ) + nc <- ncdf4::ncvar_rename(nc, "time_bnds", "time_bounds") + ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec = NA) ncdf4::nc_close(nc) - unlink(files, recursive = T) - }else{ - nc <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var) - ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec=NA) + unlink(files, recursive = TRUE) + } else { + nc <- ncdf4::nc_create( + file.path(outdir, paste(y, "nc", sep = ".")), + nc_var + ) + ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec = NA) varfile <- file(file.path(outdir, paste(y, "nc", "var", sep = ".")), "w") for (i in seq_along(nc_var)) { ncdf4::ncvar_put(nc, nc_var[[i]], output[[i]]) - cat(paste(nc_var[[i]]$name, nc_var[[i]]$longname), file = varfile, sep = "\n") + cat( + paste(nc_var[[i]]$name, nc_var[[i]]$longname), + file = varfile, + sep = "\n" + ) } close(varfile) ncdf4::nc_close(nc) } - } ### End of year loop + } ### End of year loop ## Delete raw output, if requested if (delete.raw) { file.remove(sipnet_out_file) } } # model2netcdf.SIPNET -#--------------------------------------------------------------------------------------------------# -### EOF diff --git a/models/sipnet/R/read_restart.SIPNET.R b/models/sipnet/R/read_restart.SIPNET.R index 2260c71a476..042f2aab1f4 100755 --- a/models/sipnet/R/read_restart.SIPNET.R +++ b/models/sipnet/R/read_restart.SIPNET.R @@ -1,140 +1,143 @@ ##' @title Read restart function for SDA with SIPNET -##' +##' ##' @author Ann Raiho \email{araiho@@nd.edu} -##' +##' ##' @inheritParams PEcAn.ModelName::read_restart.ModelName -##' +##' ##' @description Read Restart for SIPNET -##' +##' ##' @return X.vec vector of forecasts ##' @export read_restart.SIPNET <- function(outdir, runid, stop.time, settings, var.names, params) { - prior.sla <- params[[which(!names(params) %in% c("soil", "soil_SDA", "restart"))[1]]]$SLA - + forecast <- list() - params$restart <-c() #state.vars not in var.names will be added here - #SIPNET inital states refer to models/sipnet/inst/template.param - state.vars <- c("SWE", "SoilMoistFrac", "AbvGrndWood", "TotSoilCarb", "LAI", - "litter_carbon_content", "fine_root_carbon_content", - "coarse_root_carbon_content", "litter_mass_content_of_water") - #when adding new state variables make sure the naming is consistent across read_restart, write_restart and write.configs - #pre-populate parsm$restart with NAs so state names can be added + params$restart <- c() # state.vars not in var.names will be added here + # SIPNET inital states refer to models/sipnet/inst/template.param + state.vars <- c( + "SWE", "SoilMoistFrac", "AbvGrndWood", "TotSoilCarb", "LAI", + "litter_carbon_content", "fine_root_carbon_content", + "coarse_root_carbon_content", "litter_mass_content_of_water" + ) + # when adding new state variables make sure the naming is consistent across read_restart, write_restart and write.configs + # pre-populate parsm$restart with NAs so state names can be added params$restart <- rep(NA, length(setdiff(state.vars, var.names))) - #add states to params$restart NOT in var.names + # add states to params$restart NOT in var.names names(params$restart) <- setdiff(state.vars, var.names) # Read ensemble output - ens <- PEcAn.utils::read.output(runid = runid, - outdir = file.path(outdir, runid), - start.year = lubridate::year(stop.time), - end.year = lubridate::year(stop.time), - variables = c(state.vars,"time_bounds")) - #calculate last - start.time <- as.Date(paste0(lubridate::year(stop.time),"-01-01")) - time_var <- ens$time_bounds[1,] - real_time <- as.POSIXct(time_var*3600*24, origin = start.time) + ens <- PEcAn.utils::read.output( + runid = runid, + outdir = file.path(outdir, runid), + start.year = lubridate::year(stop.time), + end.year = lubridate::year(stop.time), + variables = c(state.vars, "time_bounds") + ) + # calculate last + start.time <- as.Date(paste0(lubridate::year(stop.time), "-01-01")) + time_var <- ens$time_bounds[1, ] + real_time <- as.POSIXct(time_var * 3600 * 24, origin = start.time) # last <- which(as.Date(real_time)==as.Date(stop.time))[1] - last <- which(as.Date(real_time)==as.Date(stop.time))[length(which(as.Date(real_time)==as.Date(stop.time)))] - + last <- which(as.Date(real_time) == as.Date(stop.time))[length(which(as.Date(real_time) == as.Date(stop.time)))] + #### PEcAn Standard Outputs if ("AbvGrndWood" %in% var.names) { - forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(ens$AbvGrndWood[last], "kg/m^2", "Mg/ha") + forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(ens$AbvGrndWood[last], "kg/m^2", "Mg/ha") names(forecast[[length(forecast)]]) <- c("AbvGrndWood") - - wood_total_C <- ens$AbvGrndWood[last] + ens$fine_root_carbon_content[last] + ens$coarse_root_carbon_content[last] - if (wood_total_C<=0) wood_total_C <- 0.0001 # Making sure we are not making Nans in case there is no plant living there. - - params$restart["abvGrndWoodFrac"] <- ens$AbvGrndWood[last] / wood_total_C - params$restart["coarseRootFrac"] <- ens$coarse_root_carbon_content[last] / wood_total_C - params$restart["fineRootFrac"] <- ens$fine_root_carbon_content[last] / wood_total_C - }else{ - params$restart["AbvGrndWood"] <- PEcAn.utils::ud_convert(ens$AbvGrndWood[last], "kg/m^2", "g/m^2") + + wood_total_C <- ens$AbvGrndWood[last] + ens$fine_root_carbon_content[last] + ens$coarse_root_carbon_content[last] + if (wood_total_C <= 0) wood_total_C <- 0.0001 # Making sure we are not making Nans in case there is no plant living there. + + params$restart["abvGrndWoodFrac"] <- ens$AbvGrndWood[last] / wood_total_C + params$restart["coarseRootFrac"] <- ens$coarse_root_carbon_content[last] / wood_total_C + params$restart["fineRootFrac"] <- ens$fine_root_carbon_content[last] / wood_total_C + } else { + params$restart["AbvGrndWood"] <- PEcAn.utils::ud_convert(ens$AbvGrndWood[last], "kg/m^2", "g/m^2") # calculate fractions, store in params, will use in write_restart - wood_total_C <- ens$AbvGrndWood[last] + ens$fine_root_carbon_content[last] + ens$coarse_root_carbon_content[last] - if (wood_total_C<=0) wood_total_C <- 0.0001 # Making sure we are not making Nans in case there is no plant living there. - - params$restart["abvGrndWoodFrac"] <- ens$AbvGrndWood[last] / wood_total_C - params$restart["coarseRootFrac"] <- ens$coarse_root_carbon_content[last] / wood_total_C - params$restart["fineRootFrac"] <- ens$fine_root_carbon_content[last] / wood_total_C + wood_total_C <- ens$AbvGrndWood[last] + ens$fine_root_carbon_content[last] + ens$coarse_root_carbon_content[last] + if (wood_total_C <= 0) wood_total_C <- 0.0001 # Making sure we are not making Nans in case there is no plant living there. + + params$restart["abvGrndWoodFrac"] <- ens$AbvGrndWood[last] / wood_total_C + params$restart["coarseRootFrac"] <- ens$coarse_root_carbon_content[last] / wood_total_C + params$restart["fineRootFrac"] <- ens$fine_root_carbon_content[last] / wood_total_C } - + if ("GWBI" %in% var.names) { - forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(mean(ens$GWBI), "kg/m^2/s", "Mg/ha/yr") + forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(mean(ens$GWBI), "kg/m^2/s", "Mg/ha/yr") names(forecast[[length(forecast)]]) <- c("GWBI") } - + # Reading in NET Ecosystem Exchange for SDA - unit is kg C m-2 s-1 and the average is estimated if ("NEE" %in% var.names) { - forecast[[length(forecast) + 1]] <- mean(ens$NEE) ## + forecast[[length(forecast) + 1]] <- mean(ens$NEE) ## names(forecast[[length(forecast)]]) <- c("NEE") } - - + + # Reading in Latent heat flux for SDA - unit is MW m-2 if ("Qle" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$Qle[last]*1e-6 ## + forecast[[length(forecast) + 1]] <- ens$Qle[last] * 1e-6 ## names(forecast[[length(forecast)]]) <- c("Qle") } - + if ("leaf_carbon_content" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$leaf_carbon_content[last] ## kgC/m2*m2/kg*2kg/kgC + forecast[[length(forecast) + 1]] <- ens$leaf_carbon_content[last] ## kgC/m2*m2/kg*2kg/kgC names(forecast[[length(forecast)]]) <- c("LeafC") } - + if ("LAI" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$LAI[last] ## m2/m2 + forecast[[length(forecast) + 1]] <- ens$LAI[last] ## m2/m2 names(forecast[[length(forecast)]]) <- c("LAI") - }else{ + } else { params$restart["LAI"] <- ens$LAI[last] } - + if ("litter_carbon_content" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$litter_carbon_content[last] ##kgC/m2 + forecast[[length(forecast) + 1]] <- ens$litter_carbon_content[last] ## kgC/m2 names(forecast[[length(forecast)]]) <- c("litter_carbon_content") - }else{ - params$restart["litter_carbon_content"] <- PEcAn.utils::ud_convert(ens$litter_carbon_content[last], 'kg m-2', 'g m-2') # kgC/m2 -> gC/m2 + } else { + params$restart["litter_carbon_content"] <- PEcAn.utils::ud_convert(ens$litter_carbon_content[last], "kg m-2", "g m-2") # kgC/m2 -> gC/m2 } - + if ("litter_mass_content_of_water" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$litter_mass_content_of_water[last] ##kgC/m2 + forecast[[length(forecast) + 1]] <- ens$litter_mass_content_of_water[last] ## kgC/m2 names(forecast[[length(forecast)]]) <- c("litter_mass_content_of_water") - }else{ + } else { params$restart["litter_mass_content_of_water"] <- ens$litter_mass_content_of_water[last] } - + if ("SoilMoistFrac" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$SoilMoistFrac[last]*100 ## here we multiply it by 100 to convert from proportion to percentage. + forecast[[length(forecast) + 1]] <- ens$SoilMoistFrac[last] * 100 ## here we multiply it by 100 to convert from proportion to percentage. names(forecast[[length(forecast)]]) <- c("SoilMoistFrac") - }else{ + } else { params$restart["SoilMoistFrac"] <- ens$SoilMoistFrac[last] } - + # This is snow if ("SWE" %in% var.names) { - forecast[[length(forecast) + 1]] <- ens$SWE[last] ## kgC/m2 + forecast[[length(forecast) + 1]] <- ens$SWE[last] ## kgC/m2 names(forecast[[length(forecast)]]) <- c("SWE") - }else{ - params$restart["SWE"] <- ens$SWE[last]/10 + } else { + params$restart["SWE"] <- ens$SWE[last] / 10 } - + if ("TotLivBiom" %in% var.names) { - forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(ens$TotLivBiom[last], "kg/m^2", "Mg/ha") + forecast[[length(forecast) + 1]] <- PEcAn.utils::ud_convert(ens$TotLivBiom[last], "kg/m^2", "Mg/ha") names(forecast[[length(forecast)]]) <- c("TotLivBiom") } - + if ("TotSoilCarb" %in% var.names) { forecast[[length(forecast) + 1]] <- ens$TotSoilCarb[last] names(forecast[[length(forecast)]]) <- c("TotSoilCarb") - }else{ - params$restart["TotSoilCarb"] <- PEcAn.utils::ud_convert(ens$TotSoilCarb[last], 'kg m-2', 'g m-2') # kgC/m2 -> gC/m2 + } else { + params$restart["TotSoilCarb"] <- PEcAn.utils::ud_convert(ens$TotSoilCarb[last], "kg m-2", "g m-2") # kgC/m2 -> gC/m2 } - - #remove any remaining NAs from params$restart + + # remove any remaining NAs from params$restart params$restart <- stats::na.omit(params$restart) - + print(runid) - + X_tmp <- list(X = unlist(forecast), params = params) - + return(X_tmp) -} # read_restart.SIPNET \ No newline at end of file +} # read_restart.SIPNET diff --git a/models/sipnet/R/sample.IC.SIPNET.R b/models/sipnet/R/sample.IC.SIPNET.R index 540a7bf5006..db8dce6eeca 100644 --- a/models/sipnet/R/sample.IC.SIPNET.R +++ b/models/sipnet/R/sample.IC.SIPNET.R @@ -10,67 +10,77 @@ #' @export #' sample.IC.SIPNET <- function(ne, state, year = 1) { - ## Mg C / ha / yr GWBI ## no conversion needed because SIPNET doesn't take GWBI as IC anyway - GWBI <- ifelse(rep("GWBI" %in% names(state), ne), - state$GWBI[sample.int(length(state$GWBI), ne)], ## unit MgC ha-1 yr-1 - stats::runif(ne, 0, 10)) ## prior - + GWBI <- ifelse(rep("GWBI" %in% names(state), ne), + state$GWBI[sample.int(length(state$GWBI), ne)], ## unit MgC ha-1 yr-1 + stats::runif(ne, 0, 10) + ) ## prior + # g C * m-2 ground area in wood (above-ground + roots) Mgha2gm <- (1000000) / (10000) # these unit conversions are for testing # reminder : when working with kgC m-2 s-1 as NPP units singularity issues pop up in sda.enkf # using MgC ha-1 yr-1 for NPP in SDA and also brought back AbvGrndWood to MgC ha-1 for sanity reasons - AbvGrndWood <- ifelse(rep("AbvGrndWood" %in% names(state), ne), - PEcAn.utils::ud_convert(state$AbvGrndWood[sample.int(length(state$AbvGrndWood), ne)], "Mg/ha", "g/m^2"), - stats::runif(ne, 700, 15000)) ## prior - + AbvGrndWood <- ifelse(rep("AbvGrndWood" %in% names(state), ne), + PEcAn.utils::ud_convert(state$AbvGrndWood[sample.int(length(state$AbvGrndWood), ne)], "Mg/ha", "g/m^2"), + stats::runif(ne, 700, 15000) + ) ## prior + # sipnet accepts a plantWoodC pool that is above-ground + roots # instead of roots having their own state, we'll pass around fractions to update them deterministically fine_root_carbon_content <- stats::runif(ne, 100, 1000) coarse_root_carbon_content <- stats::runif(ne, 200, 2000) wood_total_C <- AbvGrndWood + fine_root_carbon_content + coarse_root_carbon_content - + abvGrndWoodFrac <- AbvGrndWood / wood_total_C - coarseRootFrac <- coarse_root_carbon_content / wood_total_C - fineRootFrac <- fine_root_carbon_content / wood_total_C - + coarseRootFrac <- coarse_root_carbon_content / wood_total_C + fineRootFrac <- fine_root_carbon_content / wood_total_C + # initial leaf area, m2 leaves * m-2 ground area (multiply by leafCSpWt to ## get initial plant leaf C) lai <- ifelse(rep("LAI" %in% names(state), ne), - state$LAI[1, sample.int(ncol(state$LAI), ne), year], - stats::runif(ne, 0, 7)) ## prior + state$LAI[1, sample.int(ncol(state$LAI), ne), year], + stats::runif(ne, 0, 7) + ) ## prior ## g C * m-2 ground area litter <- ifelse(rep("litter" %in% names(state), ne), - state$litter[1, sample.int(ncol(state$litter), ne), year], - stats::runif(ne, 130, 1200)) ## prior + state$litter[1, sample.int(ncol(state$litter), ne), year], + stats::runif(ne, 130, 1200) + ) ## prior ## g C * m-2 ground area soil <- ifelse(rep("soil" %in% names(state), ne), - state$soil[1, sample.int(ncol(state$soil), ne), year], - stats::runif(ne, 1200, 2000)) ## prior + state$soil[1, sample.int(ncol(state$soil), ne), year], + stats::runif(ne, 1200, 2000) + ) ## prior ## unitless: fraction of litterWHC litterWFrac <- ifelse(rep("litterW" %in% names(state), ne), - state$litterW[1, sample.int(ncol(state$litterW), ne), year], - stats::runif(ne)) ## prior + state$litterW[1, sample.int(ncol(state$litterW), ne), year], + stats::runif(ne) + ) ## prior ## unitless: fraction of soilWHC soilWFrac <- ifelse(rep("soilW" %in% names(state), ne), - state$soilW[1, sample.int(ncol(state$soilW), ne), year], - stats::runif(ne)) ## prior + state$soilW[1, sample.int(ncol(state$soilW), ne), year], + stats::runif(ne) + ) ## prior ## cm water equiv snow <- ifelse(rep("snow" %in% names(state), ne), - state$snow[1, sample.int(ncol(state$snow), ne), year], - stats::runif(ne, 0, 2000)) ## prior + state$snow[1, sample.int(ncol(state$snow), ne), year], + stats::runif(ne, 0, 2000) + ) ## prior microbe <- ifelse(rep("microbe" %in% names(state), ne), - state$microbe[1, sample.int(ncol(state$microbe), ne), year], - stats::runif(ne, 0.02, 1)) ## prior + state$microbe[1, sample.int(ncol(state$microbe), ne), year], + stats::runif(ne, 0.02, 1) + ) ## prior - return(data.frame(GWBI, AbvGrndWood, abvGrndWoodFrac, coarseRootFrac, fineRootFrac, lai, litter, - soil, litterWFrac, soilWFrac, snow, microbe)) + return(data.frame( + GWBI, AbvGrndWood, abvGrndWoodFrac, coarseRootFrac, fineRootFrac, lai, litter, + soil, litterWFrac, soilWFrac, snow, microbe + )) } # sample.IC.SIPNET diff --git a/models/sipnet/R/split_inputs.SIPNET.R b/models/sipnet/R/split_inputs.SIPNET.R index 4bc6de7dfa6..e82b46727d7 100644 --- a/models/sipnet/R/split_inputs.SIPNET.R +++ b/models/sipnet/R/split_inputs.SIPNET.R @@ -2,7 +2,7 @@ ##' @title split_inputs.SIPNET ##' @name split_inputs.SIPNET ##' @author Mike Dietze and Ann Raiho -##' +##' ##' @param settings PEcAn settings object ##' @param start.time start date and time for each SDA ensemble ##' @param stop.time stop date and time for each SDA ensemble @@ -10,7 +10,7 @@ ##' @param overwrite Default FALSE ##' @param outpath if specified, write output to a new directory. Default NULL writes back to the directory being read ##' @description Splits climate met for SIPNET -##' +##' ##' @return file split up climate file ##' ##' @importFrom dplyr %>% @@ -20,40 +20,44 @@ split_inputs.SIPNET <- function(settings, start.time, stop.time, inputs, overwri met <- inputs path <- dirname(met) prefix <- sub(".clim", "", basename(met), fixed = TRUE) - if(is.null(outpath)){ + if (is.null(outpath)) { outpath <- path } - if(!dir.exists(outpath)) dir.create(outpath) - + if (!dir.exists(outpath)) dir.create(outpath) + file <- NA names(file) <- paste(start.time, "-", stop.time) - - #Changing the name of the files, so it would contain the name of the hour as well. - file <- paste0(outpath, "/", prefix, ".", - paste0(start.time%>% as.character() %>% gsub(' ',"_",.), - "-", - stop.time%>% as.character() %>% gsub(' ',"_",.)), ".clim") - - if(file.exists(file) & !overwrite){ + + # Changing the name of the files, so it would contain the name of the hour as well. + file <- paste0( + outpath, "/", prefix, ".", + paste0( + start.time %>% as.character() %>% gsub(" ", "_", .), + "-", + stop.time %>% as.character() %>% gsub(" ", "_", .) + ), ".clim" + ) + + if (file.exists(file) & !overwrite) { return(file) } input.dat <- utils::read.table(met, header = FALSE) - #@Hamze, I added the Date variable by using year, doy, and hour and filtered the clim based that and then removed it afterwards. - dat<-input.dat %>% - dplyr::mutate(Date = strptime(paste(V2, V3), format = "%Y %j", tz = "UTC")%>% as.POSIXct()) %>% - dplyr::mutate(Date = as.POSIXct(paste0(Date, ceiling(V4), ":00"), format = "%Y-%m-%d %H:%M", tz = "UTC")) %>% - dplyr::filter(Date >= start.time, Date < stop.time) %>% + # @Hamze, I added the Date variable by using year, doy, and hour and filtered the clim based that and then removed it afterwards. + dat <- input.dat %>% + dplyr::mutate(Date = strptime(paste(V2, V3), format = "%Y %j", tz = "UTC") %>% as.POSIXct()) %>% + dplyr::mutate(Date = as.POSIXct(paste0(Date, ceiling(V4), ":00"), format = "%Y-%m-%d %H:%M", tz = "UTC")) %>% + dplyr::filter(Date >= start.time, Date < stop.time) %>% dplyr::select(-Date) - - + + ###### Write Met to file utils::write.table(dat, file, row.names = FALSE, col.names = FALSE) ###### Output input path to inputs - #settings$run$inputs$met$path <- file + # settings$run$inputs$met$path <- file return(file) } # split_inputs.SIPNET diff --git a/models/sipnet/R/veg2model.SIPNET.R b/models/sipnet/R/veg2model.SIPNET.R index 741e37a6748..508e3463064 100644 --- a/models/sipnet/R/veg2model.SIPNET.R +++ b/models/sipnet/R/veg2model.SIPNET.R @@ -1,18 +1,17 @@ #' veg2model.SIPNET -#' +#' #' @param outfolder location to store ncdf files #' @param poolinfo object passed from write_ic contains output from cohort2pool function #' @param siteid object passed from write_ic contains site id #' @param ens number of ensemble members -#' +#' #' @return result object with filepaths to ncdf files #' @export #' @author Alexis Helgeson -#' -veg2model.SIPNET <- function(outfolder, poolinfo, siteid, ens){ - +#' +veg2model.SIPNET <- function(outfolder, poolinfo, siteid, ens) { # Execute pool_ic function -result <- PEcAn.data.land::pool_ic_list2netcdf(input = poolinfo, outdir = outfolder, siteid = siteid, ens = ens) + result <- PEcAn.data.land::pool_ic_list2netcdf(input = poolinfo, outdir = outfolder, siteid = siteid, ens = ens) -return(result) -} \ No newline at end of file + return(result) +} diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R index c294808b7b5..fee1c735a50 100755 --- a/models/sipnet/R/write.configs.SIPNET.R +++ b/models/sipnet/R/write.configs.SIPNET.R @@ -17,10 +17,10 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET") config.text <- readLines(con = template.in, n = -1) writeLines(config.text, con = file.path(settings$rundir, run.id, "sipnet.in")) - + ### WRITE *.clim - template.clim <- settings$run$inputs$met$path ## read from settings - + template.clim <- settings$run$inputs$met$path ## read from settings + if (!is.null(inputs)) { ## override if specified in inputs if ("met" %in% names(inputs)) { @@ -28,7 +28,7 @@ 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)) @@ -36,14 +36,14 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs 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) } else { jobsh <- readLines(con = system.file("template.job", package = "PEcAn.SIPNET"), n = -1) } - + # create host specific setttings hostsetup <- "" if (!is.null(settings$model$prerun)) { @@ -52,13 +52,13 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if (!is.null(settings$host$prerun)) { hostsetup <- paste(hostsetup, sep = "\n", paste(settings$host$prerun, collapse = "\n")) } - + # create cdo specific settings cdosetup <- "" if (!is.null(settings$host$cdosetup)) { cdosetup <- paste(cdosetup, sep = "\n", paste(settings$host$cdosetup, collapse = "\n")) } - + hostteardown <- "" if (!is.null(settings$model$postrun)) { hostteardown <- paste(hostteardown, sep = "\n", paste(settings$model$postrun, collapse = "\n")) @@ -66,92 +66,92 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if (!is.null(settings$host$postrun)) { hostteardown <- paste(hostteardown, sep = "\n", paste(settings$host$postrun, collapse = "\n")) } - + # create rabbitmq specific setup. cpruncmd <- cpoutcmd <- rmoutdircmd <- rmrundircmd <- "" if (!is.null(settings$host$rabbitmq)) { - #rsync cmd from remote to local host. + # rsync cmd from remote to local host. settings$host$rabbitmq$cpfcmd <- ifelse(is.null(settings$host$rabbitmq$cpfcmd), "", settings$host$rabbitmq$cpfcmd) cpruncmd <- gsub("@OUTDIR@", settings$host$rundir, settings$host$rabbitmq$cpfcmd) cpruncmd <- gsub("@OUTFOLDER@", rundir, cpruncmd) - + cpoutcmd <- gsub("@OUTDIR@", settings$host$outdir, settings$host$rabbitmq$cpfcmd) cpoutcmd <- gsub("@OUTFOLDER@", outdir, cpoutcmd) - - #delete files within rundir and outdir. + + # delete files within rundir and outdir. rmoutdircmd <- paste("rm", file.path(outdir, "*")) rmrundircmd <- paste("rm", file.path(rundir, "*")) } - + # create job.sh jobsh <- gsub("@HOST_SETUP@", hostsetup, jobsh) jobsh <- gsub("@CDO_SETUP@", cdosetup, jobsh) jobsh <- gsub("@HOST_TEARDOWN@", hostteardown, jobsh) - + jobsh <- gsub("@SITE_LAT@", settings$run$site$lat, jobsh) jobsh <- gsub("@SITE_LON@", settings$run$site$lon, jobsh) jobsh <- gsub("@SITE_MET@", template.clim, jobsh) - + jobsh <- gsub("@OUTDIR@", outdir, jobsh) jobsh <- gsub("@RUNDIR@", rundir, jobsh) - + jobsh <- gsub("@START_DATE@", settings$run$start.date, jobsh) - jobsh <- gsub("@END_DATE@",settings$run$end.date , jobsh) - + jobsh <- gsub("@END_DATE@", settings$run$end.date, jobsh) + jobsh <- gsub("@BINARY@", settings$model$binary, jobsh) jobsh <- gsub("@REVISION@", settings$model$revision, jobsh) - + jobsh <- gsub("@CPRUNCMD@", cpruncmd, jobsh) jobsh <- gsub("@CPOUTCMD@", cpoutcmd, jobsh) jobsh <- gsub("@RMOUTDIRCMD@", rmoutdircmd, jobsh) jobsh <- gsub("@RMRUNDIRCMD@", rmrundircmd, jobsh) - - if(is.null(settings$state.data.assimilation$NC.Prefix)){ + + if (is.null(settings$state.data.assimilation$NC.Prefix)) { settings$state.data.assimilation$NC.Prefix <- "sipnet.out" } jobsh <- gsub("@PREFIX@", settings$state.data.assimilation$NC.Prefix, jobsh) - - #overwrite argument - if(is.null(settings$state.data.assimilation$NC.Overwrite)){ + + # overwrite argument + if (is.null(settings$state.data.assimilation$NC.Overwrite)) { settings$state.data.assimilation$NC.Overwrite <- FALSE } jobsh <- gsub("@OVERWRITE@", settings$state.data.assimilation$NC.Overwrite, jobsh) - - #allow conflict? meaning allow full year nc export. - if(is.null(settings$state.data.assimilation$FullYearNC)){ + + # allow conflict? meaning allow full year nc export. + if (is.null(settings$state.data.assimilation$FullYearNC)) { settings$state.data.assimilation$FullYearNC <- FALSE } jobsh <- gsub("@CONFLICT@", settings$state.data.assimilation$FullYearNC, jobsh) - + if (is.null(settings$model$delete.raw)) { settings$model$delete.raw <- FALSE } 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")) - + ### WRITE *.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")) - + ### WRITE *.param template.param <- system.file("template.param", package = "PEcAn.SIPNET") if ("default.param" %in% names(settings$model)) { template.param <- settings$model$default.param } - + param <- utils::read.table(template.param) - + #### write run-specific PFT parameters here #### Get parameters being handled by PEcAn for (pft in seq_along(trait.values)) { pft.traits <- unlist(trait.values[[pft]]) pft.names <- names(pft.traits) - + ## Append/replace params specified as constants constant.traits <- unlist(defaults[[1]]$constants) constant.names <- names(constant.traits) - + # Replace matches for (i in seq_along(constant.traits)) { ind <- match(constant.names[i], pft.names) @@ -164,21 +164,21 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs pft.traits[ind] <- constant.traits[i] } } - + # Remove NAs. Constants may be specified as NA to request template defaults. Note that it is 'NA' # (character) not actual NA due to being read in as XML pft.names <- pft.names[pft.traits != "NA" & !is.na(pft.traits)] pft.traits <- pft.traits[pft.traits != "NA" & !is.na(pft.traits)] pft.traits <- as.numeric(pft.traits) - + # Leaf carbon concentration - leafC <- 0.48 #0.5 + leafC <- 0.48 # 0.5 if ("leafC" %in% pft.names) { leafC <- pft.traits[which(pft.names == "leafC")] id <- which(param[, 1] == "cFracLeaf") - param[id, 2] <- leafC * 0.01 # convert to percentage from 0 to 1 + param[id, 2] <- leafC * 0.01 # convert to percentage from 0 to 1 } - + # Specific leaf area converted to SLW SLA <- NA id <- which(param[, 1] == "leafCSpWt") @@ -188,7 +188,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs } else { SLA <- 1000 * leafC / param[id, 2] } - + # Maximum photosynthesis Amax <- NA id <- which(param[, 1] == "aMax") @@ -202,82 +202,82 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("AmaxFrac" %in% pft.names) { param[which(param[, 1] == "aMaxFrac"), 2] <- pft.traits[which(pft.names == "AmaxFrac")] } - + ### Canopy extinction coefficiet (k) if ("extinction_coefficient" %in% pft.names) { param[which(param[, 1] == "attenuation"), 2] <- pft.traits[which(pft.names == "extinction_coefficient")] } - + # Leaf respiration rate converted to baseFolRespFrac if ("leaf_respiration_rate_m2" %in% pft.names) { Rd <- pft.traits[which(pft.names == "leaf_respiration_rate_m2")] id <- which(param[, 1] == "baseFolRespFrac") - param[id, 2] <- max(min(Rd/Amax, 1), 0) + param[id, 2] <- max(min(Rd / Amax, 1), 0) } - + # Low temp threshold for photosynethsis if ("Vm_low_temp" %in% pft.names) { param[which(param[, 1] == "psnTMin"), 2] <- pft.traits[which(pft.names == "Vm_low_temp")] } - + # Opt. temp for photosynthesis if ("psnTOpt" %in% pft.names) { param[which(param[, 1] == "psnTOpt"), 2] <- pft.traits[which(pft.names == "psnTOpt")] } - + # Growth respiration factor (fraction of GPP) if ("growth_resp_factor" %in% pft.names) { param[which(param[, 1] == "growthRespFrac"), 2] <- pft.traits[which(pft.names == "growth_resp_factor")] } ### !!! NOT YET USED - #Jmax = NA - #if("Jmax" %in% pft.names){ + # Jmax = NA + # if("Jmax" %in% pft.names){ # Jmax = pft.traits[which(pft.names == 'Jmax')] ### Using Jmax scaled to 25 degC. Maybe not be the best approach - #} - - #alpha = NA - #if("quantum_efficiency" %in% pft.names){ + # } + + # alpha = NA + # if("quantum_efficiency" %in% pft.names){ # alpha = pft.traits[which(pft.names == 'quantum_efficiency')] - #} - + # } + # Half saturation of PAR. PAR at which photosynthesis occurs at 1/2 theoretical maximum (Einsteins * m^-2 ground area * day^-1). - #if(!is.na(Jmax) & !is.na(alpha)){ + # if(!is.na(Jmax) & !is.na(alpha)){ # param[which(param[,1] == "halfSatPar"),2] = Jmax/(2*alpha) ### WARNING: this is a very coarse linear approximation and needs improvement ***** ### Yes, we also need to work on doing a paired query where we have both data together. ### Once halfSatPar is calculated, need to remove Jmax and quantum_efficiency from param list so they are not included in SA - #} + # } ### !!! - + # Half saturation of PAR. PAR at which photosynthesis occurs at 1/2 theoretical maximum (Einsteins * m^-2 ground area * day^-1). # Temporary implementation until above is working. if ("half_saturation_PAR" %in% pft.names) { param[which(param[, 1] == "halfSatPar"), 2] <- pft.traits[which(pft.names == "half_saturation_PAR")] } - + # Ball-berry slomatal slope parameter m if ("stomatal_slope.BB" %in% pft.names) { id <- which(param[, 1] == "m_ballBerry") param[id, 2] <- pft.traits[which(pft.names == "stomatal_slope.BB")] } - + # Slope of VPD–photosynthesis relationship. dVpd = 1 - dVpdSlope * vpd^dVpdExp if ("dVPDSlope" %in% pft.names) { param[which(param[, 1] == "dVpdSlope"), 2] <- pft.traits[which(pft.names == "dVPDSlope")] } - + # VPD–water use efficiency relationship. dVpd = 1 - dVpdSlope * vpd^dVpdExp if ("dVpdExp" %in% pft.names) { param[which(param[, 1] == "dVpdExp"), 2] <- pft.traits[which(pft.names == "dVpdExp")] } - + # Leaf turnover rate average turnover rate of leaves, in fraction per day NOTE: read in as # per-year rate! if ("leaf_turnover_rate" %in% pft.names) { param[which(param[, 1] == "leafTurnoverRate"), 2] <- pft.traits[which(pft.names == "leaf_turnover_rate")] } - + if ("wueConst" %in% pft.names) { param[which(param[, 1] == "wueConst"), 2] <- pft.traits[which(pft.names == "wueConst")] } @@ -285,7 +285,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("veg_respiration_Q10" %in% pft.names) { param[which(param[, 1] == "vegRespQ10"), 2] <- pft.traits[which(pft.names == "veg_respiration_Q10")] } - + # Base vegetation respiration. vegetation maintenance respiration at 0 degrees C (g C respired * g^-1 plant C * day^-1) # NOTE: only counts plant wood C - leaves handled elsewhere (both above and below-ground: assumed for now to have same resp. rate) # NOTE: read in as per-year rate! @@ -294,35 +294,35 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs id <- which(param[, 1] == "baseVegResp") ## Convert from umols CO2 kg s-1 to gC g day-1 stem_resp_g <- (((pft.traits[which(pft.names == "stem_respiration_rate")]) * - (44.0096 / 1e+06) * (12.01 / 44.0096)) / 1000) * 86400 + (44.0096 / 1e+06) * (12.01 / 44.0096)) / 1000) * 86400 ## use Q10 to convert stem resp from reference of 25C to 0C param[id,2] = ## pft.traits[which(pft.names=='stem_respiration_rate')]*vegRespQ10^(-25/10) - param[id, 2] <- stem_resp_g * vegRespQ10^(-25/10) + param[id, 2] <- stem_resp_g * vegRespQ10^(-25 / 10) } - + # turnover of fine roots (per year rate) if ("root_turnover_rate" %in% pft.names) { id <- which(param[, 1] == "fineRootTurnoverRate") param[id, 2] <- pft.traits[which(pft.names == "root_turnover_rate")] } - + # fine root respiration Q10 if ("fine_root_respiration_Q10" %in% pft.names) { param[which(param[, 1] == "fineRootQ10"), 2] <- pft.traits[which(pft.names == "fine_root_respiration_Q10")] } - + # base respiration rate of fine roots (per year rate) if ("root_respiration_rate" %in% pft.names) { fineRootQ10 <- param[which(param[, 1] == "fineRootQ10"), 2] id <- which(param[, 1] == "baseFineRootResp") ## Convert from umols CO2 kg s-1 to gC g day-1 root_resp_rate_g <- (((pft.traits[which(pft.names == "root_respiration_rate")]) * - (44.0096/1e+06) * (12.01 / 44.0096)) / 1000) * 86400 + (44.0096 / 1e+06) * (12.01 / 44.0096)) / 1000) * 86400 ## use Q10 to convert stem resp from reference of 25C to 0C param[id,2] = ## pft.traits[which(pft.names=='root_respiration_rate')]*fineRootQ10^(-25/10) - param[id, 2] <- root_resp_rate_g * fineRootQ10 ^ (-25 / 10) + param[id, 2] <- root_resp_rate_g * fineRootQ10^(-25 / 10) } - + # coarse root respiration Q10 if ("coarse_root_respiration_Q10" %in% pft.names) { param[which(param[, 1] == "coarseRootQ10"), 2] <- pft.traits[which(pft.names == "coarse_root_respiration_Q10")] @@ -336,34 +336,36 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs sum_alloc <- pft.traits[which(pft.names == "root_allocation_fraction")] + pft.traits[which(pft.names == "wood_allocation_fraction")] + pft.traits[which(pft.names == "leaf_allocation_fraction")] - if(sum_alloc > 1){ + if (sum_alloc > 1) { # I want this to be a severe for now, lateer can be changed back to warning - PEcAn.logger::logger.warn("Sum of allocation parameters exceeds 1 for runid = ", run.id, - "- This won't break anything since SIPNET has internal check, but notice that such combinations might not take effect in the outputs.") + PEcAn.logger::logger.warn( + "Sum of allocation parameters exceeds 1 for runid = ", run.id, + "- This won't break anything since SIPNET has internal check, but notice that such combinations might not take effect in the outputs." + ) } } - - + + # fineRootAllocation if ("root_allocation_fraction" %in% pft.names) { param[which(param[, 1] == "fineRootAllocation"), 2] <- pft.traits[which(pft.names == "root_allocation_fraction")] } - + # woodAllocation if ("wood_allocation_fraction" %in% pft.names) { param[which(param[, 1] == "woodAllocation"), 2] <- pft.traits[which(pft.names == "wood_allocation_fraction")] } - + # leafAllocation if ("leaf_allocation_fraction" %in% pft.names) { param[which(param[, 1] == "leafAllocation"), 2] <- pft.traits[which(pft.names == "leaf_allocation_fraction")] } - + # wood_turnover_rate if ("wood_turnover_rate" %in% pft.names) { param[which(param[, 1] == "woodTurnoverRate"), 2] <- pft.traits[which(pft.names == "wood_turnover_rate")] } - + ### ----- Soil parameters soil respiration Q10. if ("soil_respiration_Q10" %in% pft.names) { param[which(param[, 1] == "soilRespQ10"), 2] <- pft.traits[which(pft.names == "soil_respiration_Q10")] @@ -372,7 +374,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("som_respiration_rate" %in% pft.names) { param[which(param[, 1] == "baseSoilResp"), 2] <- pft.traits[which(pft.names == "som_respiration_rate")] } - + # litterBreakdownRate if ("turn_over_time" %in% pft.names) { id <- which(param[, 1] == "litterBreakdownRate") @@ -382,12 +384,12 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("frozenSoilEff" %in% pft.names) { param[which(param[, 1] == "frozenSoilEff"), 2] <- pft.traits[which(pft.names == "frozenSoilEff")] } - + # frozenSoilFolREff if ("frozenSoilFolREff" %in% pft.names) { param[which(param[, 1] == "frozenSoilFolREff"), 2] <- pft.traits[which(pft.names == "frozenSoilFolREff")] } - + # soilWHC if ("soilWHC" %in% pft.names) { param[which(param[, 1] == "soilWHC"), 2] <- pft.traits[which(pft.names == "soilWHC")] @@ -395,30 +397,30 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs # 10/31/2017 IF: these were the two assumptions used in the emulator paper in order to reduce dimensionality # These results in improved winter soil respiration values # they don't affect anything when the seasonal soil respiration functionality in SIPNET is turned-off - if(TRUE){ + if (TRUE) { # assume soil resp Q10 cold == soil resp Q10 param[which(param[, 1] == "soilRespQ10Cold"), 2] <- param[which(param[, 1] == "soilRespQ10"), 2] # default SIPNET prior of baseSoilRespCold was 1/4th of baseSoilResp # assuming they will scale accordingly param[which(param[, 1] == "baseSoilRespCold"), 2] <- param[which(param[, 1] == "baseSoilResp"), 2] * 0.25 } - + if ("immedEvapFrac" %in% pft.names) { param[which(param[, 1] == "immedEvapFrac"), 2] <- pft.traits[which(pft.names == "immedEvapFrac")] } - + if ("leafWHC" %in% pft.names) { param[which(param[, 1] == "leafPoolDepth"), 2] <- pft.traits[which(pft.names == "leafWHC")] } - + if ("waterRemoveFrac" %in% pft.names) { param[which(param[, 1] == "waterRemoveFrac"), 2] <- pft.traits[which(pft.names == "waterRemoveFrac")] } - + if ("fastFlowFrac" %in% pft.names) { param[which(param[, 1] == "fastFlowFrac"), 2] <- pft.traits[which(pft.names == "fastFlowFrac")] } - + if ("rdConst" %in% pft.names) { param[which(param[, 1] == "rdConst"), 2] <- pft.traits[which(pft.names == "rdConst")] } @@ -426,56 +428,56 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("GDD" %in% pft.names) { param[which(param[, 1] == "gddLeafOn"), 2] <- pft.traits[which(pft.names == "GDD")] } - + # Fraction of leaf fall per year (should be 1 for decid) if ("fracLeafFall" %in% pft.names) { param[which(param[, 1] == "fracLeafFall"), 2] <- pft.traits[which(pft.names == "fracLeafFall")] } - + # Leaf growth. Amount of C added to the leaf during the greenup period if ("leafGrowth" %in% pft.names) { param[which(param[, 1] == "leafGrowth"), 2] <- pft.traits[which(pft.names == "leafGrowth")] } - #update LeafOnday and LeafOffDay - if (!is.null(settings$run$inputs$leaf_phenology)){ - obs_year_start <- lubridate::year(settings$run$start.date) - obs_year_end <- lubridate::year(settings$run$end.date) - if (obs_year_start != obs_year_end) { - PEcAn.logger::logger.info("Start.date and end.date are not in the same year. Currently start.date is used for refering phenological data") - } - leaf_pheno_path <- settings$run$inputs$leaf_phenology$path ## read from settings - if (!is.null(leaf_pheno_path)){ - ##read data - leafphdata <- utils::read.csv(leaf_pheno_path) - leafOnDay <- leafphdata$leafonday[leafphdata$year == obs_year_start & leafphdata$site_id==settings$run$site$id] - leafOffDay<- leafphdata$leafoffday[leafphdata$year== obs_year_start & leafphdata$site_id==settings$run$site$id] - if (!is.na(leafOnDay)){ - param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay - } - if (!is.na(leafOffDay)){ - param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay - } + # update LeafOnday and LeafOffDay + if (!is.null(settings$run$inputs$leaf_phenology)) { + obs_year_start <- lubridate::year(settings$run$start.date) + obs_year_end <- lubridate::year(settings$run$end.date) + if (obs_year_start != obs_year_end) { + PEcAn.logger::logger.info("Start.date and end.date are not in the same year. Currently start.date is used for refering phenological data") + } + leaf_pheno_path <- settings$run$inputs$leaf_phenology$path ## read from settings + if (!is.null(leaf_pheno_path)) { + ## read data + leafphdata <- utils::read.csv(leaf_pheno_path) + leafOnDay <- leafphdata$leafonday[leafphdata$year == obs_year_start & leafphdata$site_id == settings$run$site$id] + leafOffDay <- leafphdata$leafoffday[leafphdata$year == obs_year_start & leafphdata$site_id == settings$run$site$id] + if (!is.na(leafOnDay)) { + param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay + } + if (!is.na(leafOffDay)) { + param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay + } } else { - PEcAn.logger::logger.info("No phenology data were found. Please consider running `PEcAn.data.remote::extract_phenology_MODIS` to get the parameter file.") + PEcAn.logger::logger.info("No phenology data were found. Please consider running `PEcAn.data.remote::extract_phenology_MODIS` to get the parameter file.") } } } ## end loop over PFTS ####### end parameter update - #working on reading soil file (only working for 1 soil file) - if(length(settings$run$inputs$soilinitcond$path)==1){ + # working on reading soil file (only working for 1 soil file) + if (length(settings$run$inputs$soilinitcond$path) == 1) { soil_IC_list <- PEcAn.data.land::pool_ic_netcdf2list(settings$run$inputs$soilinitcond$path) - #SoilWHC and LitterWHC - if("volume_fraction_of_water_in_soil_at_saturation"%in%names(soil_IC_list$vals)){ - #SoilWHC - param[which(param[, 1] == "soilWHC"), 2] <- mean(unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"]))*100 - - #LitterWHC - #param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*100 - } - if("soil_hydraulic_conductivity_at_saturation"%in%names(soil_IC_list$vals)){ - #litwaterDrainrate - param[which(param[, 1] == "litWaterDrainRate"), 2] <- unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1]*100/(3600*24) + # SoilWHC and LitterWHC + if ("volume_fraction_of_water_in_soil_at_saturation" %in% names(soil_IC_list$vals)) { + # SoilWHC + param[which(param[, 1] == "soilWHC"), 2] <- mean(unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])) * 100 + + # LitterWHC + # param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*100 + } + if ("soil_hydraulic_conductivity_at_saturation" %in% names(soil_IC_list$vals)) { + # litwaterDrainrate + param[which(param[, 1] == "litWaterDrainRate"), 2] <- unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1] * 100 / (3600 * 24) } } if (!is.null(IC)) { @@ -484,16 +486,16 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs plant_wood_vars <- c("AbvGrndWood", "abvGrndWoodFrac", "coarseRootFrac", "fineRootFrac") if (all(plant_wood_vars %in% ic.names)) { # reconstruct total wood C - if(IC$abvGrndWoodFrac < 0.05){ + if (IC$abvGrndWoodFrac < 0.05) { wood_total_C <- IC$AbvGrndWood - }else{ + } else { wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac } - #Sanity check + # Sanity check if (is.infinite(wood_total_C) | is.nan(wood_total_C) | wood_total_C < 0) { wood_total_C <- 0 - if (round(IC$AbvGrndWood) > 0 & round(IC$abvGrndWoodFrac, 3) == 0) + if (round(IC$AbvGrndWood) > 0 & round(IC$abvGrndWoodFrac, 3) == 0) { PEcAn.logger::logger.warn( paste0( "There is a major problem with ", @@ -506,9 +508,10 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs ) ) } - param[which(param[, 1] == "plantWoodInit"), 2] <- wood_total_C + } + param[which(param[, 1] == "plantWoodInit"), 2] <- wood_total_C param[which(param[, 1] == "coarseRootFrac"), 2] <- IC$coarseRootFrac - param[which(param[, 1] == "fineRootFrac"), 2] <- IC$fineRootFrac + param[which(param[, 1] == "fineRootFrac"), 2] <- IC$fineRootFrac } ## laiInit m2/m2 if ("lai" %in% ic.names) { @@ -524,8 +527,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs } ## litterWFracInit fraction if ("litter_mass_content_of_water" %in% ic.names) { - #here we use litterWaterContent/litterWHC to calculate the litterWFracInit - param[which(param[, 1] == "litterWFracInit"), 2] <- IC$litter_mass_content_of_water/(param[which(param[, 1] == "litterWHC"), 2]*10) + # here we use litterWaterContent/litterWHC to calculate the litterWFracInit + param[which(param[, 1] == "litterWFracInit"), 2] <- IC$litter_mass_content_of_water / (param[which(param[, 1] == "litterWHC"), 2] * 10) } ## soilWFracInit fraction if ("soilWFrac" %in% ic.names) { @@ -539,16 +542,14 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs if ("microbe" %in% ic.names) { param[which(param[, 1] == "microbeInit"), 2] <- IC$microbe } - } - - else if (length(settings$run$inputs$poolinitcond$path)>0) { + } else if (length(settings$run$inputs$poolinitcond$path) > 0) { ICs_num <- length(settings$run$inputs$poolinitcond$path) IC.path <- settings$run$inputs$poolinitcond$path[[sample(1:ICs_num, 1)]] IC.pools <- PEcAn.data.land::prepare_pools(IC.path, constants = list(sla = SLA)) - - if(!is.null(IC.pools)){ - IC.nc <- ncdf4::nc_open(IC.path) #for additional variables specific to SIPNET + + if (!is.null(IC.pools)) { + IC.nc <- ncdf4::nc_open(IC.path) # for additional variables specific to SIPNET ## plantWoodInit gC/m2 if ("wood" %in% names(IC.pools)) { param[which(param[, 1] == "plantWoodInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$wood, "kg m-2", "g m-2") @@ -559,81 +560,82 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs param[which(param[, 1] == "laiInit"), 2] <- lai } - #Initial LAI is set as 0 for deciduous forests and grasslands for non-growing seasons - if (!(lubridate::month(settings$run$start.date) %in% seq(5,9))){ #Growing seasons are coarsely defined as months from May to September for non-conifers in the US - site_pft <- utils::read.csv(settings$run$inputs$pft.site$path) - site.pft.name <- site_pft$pft[site_pft$site == settings$run$site$id] - if (site.pft.name!="boreal.coniferous") { #Currently only excluding boreal conifers. Other evergreen PFTs could be added here later. - param[which(param[, 1] == "laiInit"), 2] <- 0 - } + # Initial LAI is set as 0 for deciduous forests and grasslands for non-growing seasons + if (!(lubridate::month(settings$run$start.date) %in% seq(5, 9))) { # Growing seasons are coarsely defined as months from May to September for non-conifers in the US + site_pft <- utils::read.csv(settings$run$inputs$pft.site$path) + site.pft.name <- site_pft$pft[site_pft$site == settings$run$site$id] + if (site.pft.name != "boreal.coniferous") { # Currently only excluding boreal conifers. Other evergreen PFTs could be added here later. + param[which(param[, 1] == "laiInit"), 2] <- 0 + } } ## neeInit gC/m2 - nee <- try(ncdf4::ncvar_get(IC.nc,"nee"),silent = TRUE) + nee <- try(ncdf4::ncvar_get(IC.nc, "nee"), silent = TRUE) if (!is.na(nee) && is.numeric(nee)) { param[which(param[, 1] == "neeInit"), 2] <- nee } ## litterInit gC/m2 if ("litter" %in% names(IC.pools)) { - param[which(param[, 1] == "litterInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$litter, 'g m-2', 'g m-2') # BETY: kgC m-2 + param[which(param[, 1] == "litterInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$litter, "g m-2", "g m-2") # BETY: kgC m-2 } ## soilInit gC/m2 if ("soil" %in% names(IC.pools)) { - param[which(param[, 1] == "soilInit"), 2] <- PEcAn.utils::ud_convert(sum(IC.pools$soil), 'kg m-2', 'g m-2') # BETY: kgC m-2 + param[which(param[, 1] == "soilInit"), 2] <- PEcAn.utils::ud_convert(sum(IC.pools$soil), "kg m-2", "g m-2") # BETY: kgC m-2 } ## soilWFracInit fraction - soilWFrac <- try(ncdf4::ncvar_get(IC.nc,"SoilMoistFrac"),silent = TRUE) + soilWFrac <- try(ncdf4::ncvar_get(IC.nc, "SoilMoistFrac"), silent = TRUE) if (!"try-error" %in% class(soilWFrac)) { if (!is.na(soilWFrac) && is.numeric(soilWFrac)) { - param[which(param[, 1] == "soilWFracInit"), 2] <- sum(soilWFrac)/100 + param[which(param[, 1] == "soilWFracInit"), 2] <- sum(soilWFrac) / 100 } } ## litterWFracInit fraction litterWFrac <- soilWFrac - + ## snowInit cm water equivalent (cm = g / cm2 because 1 g water = 1 cm3 water) - snow = try(ncdf4::ncvar_get(IC.nc,"SWE"),silent = TRUE) + snow <- try(ncdf4::ncvar_get(IC.nc, "SWE"), silent = TRUE) if (!is.na(snow) && is.numeric(snow)) { - param[which(param[, 1] == "snowInit"), 2] <- PEcAn.utils::ud_convert(snow, "kg m-2", "g cm-2") # BETY: kg m-2 + param[which(param[, 1] == "snowInit"), 2] <- PEcAn.utils::ud_convert(snow, "kg m-2", "g cm-2") # BETY: kg m-2 } ## leafOnDay - leafOnDay <- try(ncdf4::ncvar_get(IC.nc,"date_of_budburst"),silent = TRUE) + leafOnDay <- try(ncdf4::ncvar_get(IC.nc, "date_of_budburst"), silent = TRUE) if (!is.na(leafOnDay) && is.numeric(leafOnDay)) { param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay } ## leafOffDay - leafOffDay <- try(ncdf4::ncvar_get(IC.nc,"date_of_senescence"),silent = TRUE) + leafOffDay <- try(ncdf4::ncvar_get(IC.nc, "date_of_senescence"), silent = TRUE) if (!is.na(leafOffDay) && is.numeric(leafOffDay)) { param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay } - microbe <- try(ncdf4::ncvar_get(IC.nc,"Microbial Biomass C"),silent = TRUE) + microbe <- try(ncdf4::ncvar_get(IC.nc, "Microbial Biomass C"), silent = TRUE) if (!is.na(microbe) && is.numeric(microbe)) { - param[which(param[, 1] == "microbeInit"), 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") #BETY: mg microbial C kg-1 soil + param[which(param[, 1] == "microbeInit"), 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") # BETY: mg microbial C kg-1 soil } - + ncdf4::nc_close(IC.nc) - }else{ + } else { PEcAn.logger::logger.error("Bad initial conditions filepath; keeping defaults") } - }else{ - #some stuff about IC file that we can give in lieu of actual ICs + } else { + # some stuff about IC file that we can give in lieu of actual ICs } - - + + if (!is.null(settings$run$inputs$soilmoisture)) { - #read soil moisture netcdf file, grab closet date to start_date, set equal to soilWFrac - if(!is.null(settings$run$inputs$soilmoisture$path)){ + # read soil moisture netcdf file, grab closet date to start_date, set equal to soilWFrac + if (!is.null(settings$run$inputs$soilmoisture$path)) { soil.path <- settings$run$inputs$soilmoisture$path soilWFrac <- ncdf4::ncvar_get(ncdf4::nc_open(soil.path), varid = "mass_fraction_of_unfrozen_water_in_soil_moisture") - + param[which(param[, 1] == "soilWFracInit"), 2] <- soilWFrac } - } - if(file.exists(file.path(settings$rundir, run.id, "sipnet.param"))) file.rename(file.path(settings$rundir, run.id, "sipnet.param"),file.path(settings$rundir, run.id, paste0("sipnet_",lubridate::year(settings$run$start.date),"_",lubridate::year(settings$run$end.date),".param"))) - + if (file.exists(file.path(settings$rundir, run.id, "sipnet.param"))) file.rename(file.path(settings$rundir, run.id, "sipnet.param"), file.path(settings$rundir, run.id, paste0("sipnet_", lubridate::year(settings$run$start.date), "_", lubridate::year(settings$run$end.date), ".param"))) - utils::write.table(param, file.path(settings$rundir, run.id, "sipnet.param"), row.names = FALSE, col.names = FALSE, - quote = FALSE) + + utils::write.table(param, file.path(settings$rundir, run.id, "sipnet.param"), + row.names = FALSE, col.names = FALSE, + quote = FALSE + ) } # write.config.SIPNET #--------------------------------------------------------------------------------------------------# ##' @@ -648,20 +650,19 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs ##' ##' @author Shawn Serbin, David LeBauer remove.config.SIPNET <- function(main.outdir, settings) { - ### Remove files on localhost if (settings$host$name == "localhost") { - files <- paste0(settings$outdir, list.files(path = settings$outdir, recursive = FALSE)) # Need to change this to the run folder when implemented - files <- files[-grep("*.xml", files)] # Keep pecan.xml file + files <- paste0(settings$outdir, list.files(path = settings$outdir, recursive = FALSE)) # Need to change this to the run folder when implemented + files <- files[-grep("*.xml", files)] # Keep pecan.xml file pft.dir <- strsplit(settings$pfts$pft$outdir, "/")[[1]] ln <- length(pft.dir) pft.dir <- pft.dir[ln] - files <- files[-grep(pft.dir, files)] # Keep pft folder + files <- files[-grep(pft.dir, files)] # Keep pft folder # file.remove(files,recursive=TRUE) - system(paste("rm -r ", files, sep = "", collapse = " "), ignore.stderr = TRUE) # remove files/dirs - + system(paste("rm -r ", files, sep = "", collapse = " "), ignore.stderr = TRUE) # remove files/dirs + ### On remote host } else { print("*** WARNING: Removal of files on remote host not yet implemented ***") } -} # remove.config.SIPNET +} # remove.config.SIPNET diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R index 32d2736312e..ab7193b626b 100755 --- a/models/sipnet/R/write_restart.SIPNET.R +++ b/models/sipnet/R/write_restart.SIPNET.R @@ -12,45 +12,46 @@ ##' @param settings PEcAn settings object ##' @param new.state analysis state vector ##' @param RENAME flag to either rename output file or not -##' @param new.params list of parameters to convert between different states +##' @param new.params list of parameters to convert between different states ##' @param inputs list of model inputs to use in write.configs.SIPNET ##' @param verbose decide if we want to print the outputs. -##' +##' ##' @return NONE ##' ##' @importFrom dplyr %>% ##' @export write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state, RENAME = TRUE, new.params = FALSE, inputs, verbose = FALSE) { - rundir <- settings$host$rundir variables <- colnames(new.state) # values that will be used for updating other states deterministically depending on the SDA states if (length(new.params$restart) > 0) { IC_extra <- data.frame(t(new.params$restart)) - } else{ + } else { IC_extra <- data.frame() - } - + } + if (RENAME) { - file.rename(file.path(outdir, runid, "sipnet.out"), - file.path(outdir, runid, paste0("sipnet.", as.Date(start.time), ".out"))) + file.rename( + file.path(outdir, runid, "sipnet.out"), + file.path(outdir, runid, paste0("sipnet.", as.Date(start.time), ".out")) + ) system(paste("rm", file.path(rundir, runid, "sipnet.clim"))) } else { print(paste("Files not renamed -- Need to rerun timestep", start.time, "before next time step")) } - + settings$run$start.date <- start.time settings$run$end.date <- stop.time - + ## Converting to sipnet units prior.sla <- new.params[[which(!names(new.params) %in% c("soil", "soil_SDA", "restart"))[1]]]$SLA - unit.conv <- 2 * (10000 / 1) * (1 / 1000) * (3.154 * 10^7) # kgC/m2/s -> Mg/ha/yr - + unit.conv <- 2 * (10000 / 1) * (1 / 1000) * (3.154 * 10^7) # kgC/m2/s -> Mg/ha/yr + analysis.save <- list() - + if ("NPP" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr") #*unit.conv -> Mg/ha/yr + analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr") #* unit.conv -> Mg/ha/yr names(analysis.save[[length(analysis.save)]]) <- c("NPP") } @@ -58,60 +59,60 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, analysis.save[[length(analysis.save) + 1]] <- new.state$NEE names(analysis.save[[length(analysis.save)]]) <- c("NEE") } - + if ("AbvGrndWood" %in% variables) { - AbvGrndWood <- PEcAn.utils::ud_convert(new.state$AbvGrndWood, "Mg/ha", "g/m^2") + AbvGrndWood <- PEcAn.utils::ud_convert(new.state$AbvGrndWood, "Mg/ha", "g/m^2") analysis.save[[length(analysis.save) + 1]] <- AbvGrndWood names(analysis.save[[length(analysis.save)]]) <- c("AbvGrndWood") } - + if ("LeafC" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- new.state$LeafC * prior.sla * 2 ## kgC/m2*m2/kg*2kg/kgC -> m2/m2 + analysis.save[[length(analysis.save) + 1]] <- new.state$LeafC * prior.sla * 2 ## kgC/m2*m2/kg*2kg/kgC -> m2/m2 if (new.state$LeafC < 0) analysis.save[[length(analysis.save)]] <- 0 names(analysis.save[[length(analysis.save)]]) <- c("lai") } - + if ("litter_carbon_content" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$litter_carbon_content, 'kg m-2', 'g m-2') # kgC/m2 -> gC/m2 + analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$litter_carbon_content, "kg m-2", "g m-2") # kgC/m2 -> gC/m2 if (new.state$litter_carbon_content < 0) analysis.save[[length(analysis.save)]] <- 0 names(analysis.save[[length(analysis.save)]]) <- c("litter_carbon_content") } - + if ("TotSoilCarb" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$TotSoilCarb, 'kg m-2', 'g m-2') # kgC/m2 -> gC/m2 + analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$TotSoilCarb, "kg m-2", "g m-2") # kgC/m2 -> gC/m2 if (new.state$TotSoilCarb < 0) analysis.save[[length(analysis.save)]] <- 0 names(analysis.save[[length(analysis.save)]]) <- c("soil") } - - if("litter_mass_content_of_water" %in% variables){ - analysis.save[[length(analysis.save) + 1]] <- new.state$litter_mass_content_of_water ## unitless + + if ("litter_mass_content_of_water" %in% variables) { + analysis.save[[length(analysis.save) + 1]] <- new.state$litter_mass_content_of_water ## unitless if (new.state$litter_mass_content_of_water < 0 || new.state$litter_mass_content_of_water > 1) analysis.save[[length(analysis.save)]] <- 0.5 names(analysis.save[[length(analysis.save)]]) <- c("litter_mass_content_of_water") } - + if ("SoilMoistFrac" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- new.state$SoilMoistFrac/100 ## unitless + analysis.save[[length(analysis.save) + 1]] <- new.state$SoilMoistFrac / 100 ## unitless if (analysis.save[[length(analysis.save)]] < 0 || analysis.save[[length(analysis.save)]] > 1) analysis.save[[length(analysis.save)]] <- 0.5 names(analysis.save[[length(analysis.save)]]) <- c("soilWFrac") } - + if ("SWE" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- new.state$SWE/10 + analysis.save[[length(analysis.save) + 1]] <- new.state$SWE / 10 if (analysis.save[[length(analysis.save)]] < 0) analysis.save[[length(analysis.save)]] <- 0 names(analysis.save[[length(analysis.save)]]) <- c("SWE") } if ("LAI" %in% variables) { - analysis.save[[length(analysis.save) + 1]] <- new.state$LAI + analysis.save[[length(analysis.save) + 1]] <- new.state$LAI if (new.state$LAI < 0) analysis.save[[length(analysis.save)]] <- 0 names(analysis.save[[length(analysis.save)]]) <- c("lai") } - - if (!is.null(analysis.save) && length(analysis.save)>0){ + + if (!is.null(analysis.save) && length(analysis.save) > 0) { analysis.save.mat <- data.frame(matrix(unlist(analysis.save, use.names = TRUE), nrow = 1)) colnames(analysis.save.mat) <- names(unlist(analysis.save)) - analysis.save.mat <- cbind(analysis.save.mat, IC_extra) #add in all restart values - }else{ + analysis.save.mat <- cbind(analysis.save.mat, IC_extra) # add in all restart values + } else { analysis.save.mat <- NULL } @@ -119,11 +120,13 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, print(runid %>% as.character()) print(analysis.save.mat) } - do.call(write.config.SIPNET, args = list(defaults = NULL, - trait.values = new.params, - settings = settings, - run.id = runid, - inputs = inputs, - IC = analysis.save.mat)) + do.call(write.config.SIPNET, args = list( + defaults = NULL, + trait.values = new.params, + settings = settings, + run.id = runid, + inputs = inputs, + IC = analysis.save.mat + )) print(runid) -} # write_restart.SIPNET \ No newline at end of file +} # write_restart.SIPNET diff --git a/models/sipnet/inst/SIPNET.lyford.SDA.R b/models/sipnet/inst/SIPNET.lyford.SDA.R index 0bd93acece5..cbff691e95c 100644 --- a/models/sipnet/inst/SIPNET.lyford.SDA.R +++ b/models/sipnet/inst/SIPNET.lyford.SDA.R @@ -6,11 +6,11 @@ options(warn = 1, keep.source = TRUE, error = quote({ })) status.start <- function(name) { - cat(paste(name, format(Sys.time(), "%F %T"), sep="\t"), file=file.path(settings$outdir, "STATUS"), append=TRUE) + cat(paste(name, format(Sys.time(), "%F %T"), sep = "\t"), file = file.path(settings$outdir, "STATUS"), append = TRUE) } -status.end <- function(status="DONE") { - cat(paste("", format(Sys.time(), "%F %T"), status, "\n", sep="\t"), file=file.path(settings$outdir, "STATUS"), append=TRUE) +status.end <- function(status = "DONE") { + cat(paste("", format(Sys.time(), "%F %T"), status, "\n", sep = "\t"), file = file.path(settings$outdir, "STATUS"), append = TRUE) } #---------------- Load libraries. -----------------------------------------------------------------# @@ -23,7 +23,7 @@ library(rjags) library(reshape2) #--------------------------------------------------------------------------------------------------# # -# +# # # 35 # FALSE @@ -65,22 +65,24 @@ library(reshape2) #---------------- Load PEcAn settings file. -------------------------------------------------------# # Open and read in settings file for PEcAn run. -settings <- read.settings("pecan.SDA.xml") +settings <- read.settings("pecan.SDA.xml") #--------------------------------------------------------------------------------------------------# #---------------- Load data. -------------------------------------------------------# -load('~/sipnet_lyford_summary.Rdata') -years<-1962:2015 -names(obs.mean) <- paste0(years,'/12/31') +load("~/sipnet_lyford_summary.Rdata") +years <- 1962:2015 +names(obs.mean) <- paste0(years, "/12/31") #---------------- Build Initial Conditions ----------------------------------------------------------------------# status.start("IC") -ne = as.numeric(settings$state.data.assimilation$n.ensemble) -IC = sample.IC.SIPNET(ne,state,year=1) +ne <- as.numeric(settings$state.data.assimilation$n.ensemble) +IC <- sample.IC.SIPNET(ne, state, year = 1) status.end() #--------------- Assimilation -------------------------------------------------------# status.start("EnKF") -sda.enkf(settings=settings, obs.mean = obs.mean, - obs.cov = obs.cov, IC = IC, Q = NULL) -status.end() \ No newline at end of file +sda.enkf( + settings = settings, obs.mean = obs.mean, + obs.cov = obs.cov, IC = IC, Q = NULL +) +status.end() diff --git a/models/sipnet/man/mergeNC.Rd b/models/sipnet/man/mergeNC.Rd index 011a8c8e46f..1cf69f87ae1 100644 --- a/models/sipnet/man/mergeNC.Rd +++ b/models/sipnet/man/mergeNC.Rd @@ -10,7 +10,8 @@ https://github.com/RS-eco/processNC/blob/main/R/mergeNC.R mergeNC(files, outfile) } \arguments{ -\item{files}{\code{character}. List of filepaths, which should lead to NetCDF files.} +\item{files}{\code{character}. List of filepaths, which should lead to +NetCDF files.} \item{outfile}{\code{character}. Output filename of the merged data.} } @@ -22,10 +23,11 @@ Merge multiple NetCDF files into one } \examples{ \dontrun{ -files <- list.files(paste0(system.file(package="processNC"), "/extdata"), - pattern="tas.*\\\\.nc", full.names=TRUE) -temp <- tempfile(fileext=".nc") -mergeNC(files=files, outfile=temp) -terra::rast(temp) +files <- list.files(paste0(system.file(package = "processNC"), "/extdata"), + pattern = "tas.*\\\\.nc", full.names = TRUE +) +temp <- tempfile(fileext = ".nc") +mergeNC(files = files, outfile = temp) +terra::rast(temp) } } diff --git a/models/sipnet/man/model2netcdf.SIPNET.Rd b/models/sipnet/man/model2netcdf.SIPNET.Rd index b88a96198cf..657ba430854 100644 --- a/models/sipnet/man/model2netcdf.SIPNET.Rd +++ b/models/sipnet/man/model2netcdf.SIPNET.Rd @@ -36,7 +36,8 @@ model2netcdf.SIPNET( \item{overwrite}{Flag for overwriting nc files or not} -\item{conflict}{Flag for dealing with conflicted nc files, if T we then will merge those, if F we will jump to the next.} +\item{conflict}{Flag for dealing with conflicted nc files, +if T we then will merge those, if F we will jump to the next.} } \description{ Converts all output contained in a folder to netCDF. diff --git a/models/sipnet/tests/testthat.R b/models/sipnet/tests/testthat.R index 7a1fc04289f..c750f864487 100644 --- a/models/sipnet/tests/testthat.R +++ b/models/sipnet/tests/testthat.R @@ -2,4 +2,4 @@ library(testthat) library(PEcAn.utils) PEcAn.logger::logger.setQuitOnSevere(FALSE) -#test_check("PEcAn.SIPNET") +# test_check("PEcAn.SIPNET") diff --git a/models/sipnet/tests/testthat/test.met2model.R b/models/sipnet/tests/testthat/test.met2model.R index 9619f05e3ac..7dc5d4f0687 100644 --- a/models/sipnet/tests/testthat/test.met2model.R +++ b/models/sipnet/tests/testthat/test.met2model.R @@ -6,7 +6,8 @@ teardown(unlink(outfolder, recursive = TRUE)) test_that("Met conversion runs without error", { nc_path <- system.file("test-data", "CRUNCEP.2000.nc", - package = "PEcAn.utils") + package = "PEcAn.utils" + ) in.path <- dirname(nc_path) in.prefix <- "CRUNCEP" start_date <- "2000-01-01" From 6b95332440b919daaf8f1b837da05141eed41329 Mon Sep 17 00:00:00 2001 From: Chris Black Date: Sat, 4 Jan 2025 13:27:52 -0800 Subject: [PATCH 4/4] run styler on workflow pkg --- base/workflow/R/create_execute_test_xml.R | 106 ++++++----- base/workflow/R/do_conversions.R | 76 ++++---- base/workflow/R/run.write.configs.R | 149 ++++++++------- base/workflow/R/runModule.get.trait.data.R | 12 +- base/workflow/R/runModule.run.write.configs.R | 4 +- base/workflow/R/start_model_runs.R | 174 +++++++++--------- base/workflow/man/start_model_runs.Rd | 2 +- .../testthat/test-runModule.get.trait.data.R | 11 +- .../tests/testthat/test.do_conversions.R | 10 +- .../tests/testthat/test.start_model_runs.R | 1 - 10 files changed, 294 insertions(+), 251 deletions(-) diff --git a/base/workflow/R/create_execute_test_xml.R b/base/workflow/R/create_execute_test_xml.R index aacb05d8b08..9317ac9f208 100644 --- a/base/workflow/R/create_execute_test_xml.R +++ b/base/workflow/R/create_execute_test_xml.R @@ -46,15 +46,14 @@ create_execute_test_xml <- function(model_id, db_bety_hostname = NULL, db_bety_port = NULL, db_bety_driver = "Postgres") { - php_file <- file.path(pecan_path, "web", "config.php") config.list <- PEcAn.utils::read_web_config(php_file) if (is.null(db_bety_username)) db_bety_username <- config.list$db_bety_username if (is.null(db_bety_password)) db_bety_password <- config.list$db_bety_password if (is.null(db_bety_hostname)) db_bety_hostname <- config.list$db_bety_hostname if (is.null(db_bety_port)) db_bety_port <- config.list$db_bety_port - - #opening a connection to bety + + # opening a connection to bety con <- PEcAn.DB::db.open(list( user = db_bety_username, password = db_bety_password, @@ -65,17 +64,19 @@ create_execute_test_xml <- function(model_id, on.exit(PEcAn.DB::db.close(con), add = TRUE) settings <- list( - info = list(notes = "Test_Run", - userid = user_id, - username = "None", - dates = Sys.Date()) + info = list( + notes = "Test_Run", + userid = user_id, + username = "None", + dates = Sys.Date() + ) ) - #Outdir + # Outdir model.new <- dplyr::tbl(con, "models") %>% dplyr::filter(.data$id == !!model_id) %>% dplyr::collect() - + outdir_pre <- paste( model.new[["model_name"]], format(as.Date(start_date), "%Y-%m"), @@ -90,75 +91,84 @@ create_execute_test_xml <- function(model_id, outdir <- normalizePath(outdir) settings$outdir <- outdir - #Database BETY + # Database BETY settings$database <- list( - bety = list(user = db_bety_username, - password = db_bety_password, - host = db_bety_hostname, - dbname = "bety", - driver = db_bety_driver, - write = FALSE), + bety = list( + user = db_bety_username, + password = db_bety_password, + host = db_bety_hostname, + dbname = "bety", + driver = db_bety_driver, + write = FALSE + ), dbfiles = dbfiles_folder ) - #PFT - if (is.null(pft)){ + # PFT + if (is.null(pft)) { # Select the first PFT in the model list. pft <- dplyr::tbl(con, "pfts") %>% dplyr::filter(.data$modeltype_id == !!model.new$modeltype_id) %>% dplyr::collect() - + pft <- pft$name[[1]] - message("PFT is `NULL`. Defaulting to the following PFT: ", - pft) + message( + "PFT is `NULL`. Defaulting to the following PFT: ", + pft + ) } ## Putting multiple PFTs separated by semicolon settings$pfts <- strsplit(pft, ";")[[1]] %>% - purrr::map( ~ list(name = .x, - constants = list(num = 1) - ) - ) %>% - stats::setNames(rep("pft", length(.data))) + purrr::map(~ list( + name = .x, + constants = list(num = 1) + )) %>% + stats::setNames(rep("pft", length(.data))) - #Meta Analysis + # Meta Analysis settings$meta.analysis <- list(iter = 3000, random.effects = FALSE) - #Ensemble + # Ensemble settings$ensemble <- list( size = ensemble_size, variable = sensitivity_variable, - samplingspace = list(met = list(method = "sampling"), - parameters = list(method = "uniform")) + samplingspace = list( + met = list(method = "sampling"), + parameters = list(method = "uniform") + ) ) - #Sensitivity + # Sensitivity if (sensitivity) { settings$sensitivity.analysis <- list( quantiles = list(sigma1 = -2, sigma2 = -1, sigma3 = 1, sigma4 = 2) ) } - #Model + # Model settings$model$id <- model.new[["id"]] - #Workflow + # Workflow settings$workflow$id - settings$workflow$id <- paste0("Test_run_","_",model.new$model_name) + settings$workflow$id <- paste0("Test_run_", "_", model.new$model_name) settings$run <- list( site = list(id = site_id, met.start = start_date, met.end = end_date), - inputs = list(met = list(source = met, output = model.new[["model_name"]], - username = "pecan")), + inputs = list(met = list( + source = met, output = model.new[["model_name"]], + username = "pecan" + )), start.date = start_date, end.date = end_date ) settings$host$name <- "localhost" - + # Add model specific options - settings<-model_specific_tags(settings, model.new) - #create file and Run + settings <- model_specific_tags(settings, model.new) + # create file and Run XML::saveXML(PEcAn.settings::listToXml(settings, "pecan"), - file = file.path(outdir, "pecan.xml")) + file = file.path(outdir, "pecan.xml") + ) file.copy(file.path(pecan_path, "web", "workflow.R"), outdir) cwd <- getwd() setwd(outdir) @@ -180,14 +190,14 @@ create_execute_test_xml <- function(model_id, #' @return updated settings list #' @export #' -model_specific_tags <- function(settings, model.info){ - - #some extra settings for LPJ-GUESS - if(model.info$model_name=="LPJ-GUESS"){ - settings$run$inputs <- c(settings$run$inputs , - list(soil=list(id=1000000903)) - ) +model_specific_tags <- function(settings, model.info) { + # some extra settings for LPJ-GUESS + if (model.info$model_name == "LPJ-GUESS") { + settings$run$inputs <- c( + settings$run$inputs, + list(soil = list(id = 1000000903)) + ) } - + return(settings) } diff --git a/base/workflow/R/do_conversions.R b/base/workflow/R/do_conversions.R index 50fd71e812e..027ded6234f 100644 --- a/base/workflow/R/do_conversions.R +++ b/base/workflow/R/do_conversions.R @@ -12,61 +12,61 @@ do_conversions <- function(settings, overwrite.met = FALSE, overwrite.fia = FALS if (PEcAn.settings::is.MultiSettings(settings)) { return(PEcAn.settings::papply(settings, do_conversions)) } - + needsave <- FALSE if (is.character(settings$run$inputs)) { - settings$run$inputs <- NULL ## check for empty set + settings$run$inputs <- NULL ## check for empty set } - + dbfiles.local <- settings$database$dbfiles dbfiles <- ifelse(!PEcAn.remote::is.localhost(settings$host) & !is.null(settings$host$folder), settings$host$folder, dbfiles.local) - PEcAn.logger::logger.debug("do.conversion outdir",dbfiles) + PEcAn.logger::logger.debug("do.conversion outdir", dbfiles) # For each input for (i in seq_along(settings$run$inputs)) { input <- settings$run$inputs[[i]] if (is.null(input)) { next } - + input.tag <- names(settings$run$input)[i] - PEcAn.logger::logger.info("PROCESSING: ",input.tag) - - + PEcAn.logger::logger.info("PROCESSING: ", input.tag) + + ic.flag <- fia.flag <- FALSE - - if ((input.tag %in% c("css", "pss", "site")) && - is.null(input$path) && !is.null(input$source)) { - if(!is.null(input$useic)){ # set TRUE if IC Workflow, leave empty if not - ic.flag <- TRUE - }else if(input$source == "FIA"){ + + if ((input.tag %in% c("css", "pss", "site")) && + is.null(input$path) && !is.null(input$source)) { + if (!is.null(input$useic)) { # set TRUE if IC Workflow, leave empty if not + ic.flag <- TRUE + } else if (input$source == "FIA") { fia.flag <- TRUE # possibly a warning for deprecation in the future } } - + # BADM IC - if(input.tag == "poolinitcond" && is.null(input$path)){ - ic.flag <- TRUE + if (input.tag == "poolinitcond" && is.null(input$path)) { + ic.flag <- TRUE } - + # IC conversion : for now for ED only, hence the css/pss/site check # TRUE if (ic.flag) { - settings <- PEcAn.data.land::ic_process(settings, input, dir = dbfiles, overwrite = overwrite.ic) + settings <- PEcAn.data.land::ic_process(settings, input, dir = dbfiles, overwrite = overwrite.ic) needsave <- TRUE } - + # keep fia.to.psscss if (fia.flag) { settings <- PEcAn.data.land::fia.to.psscss(settings, overwrite = overwrite.fia) needsave <- TRUE } - - + + # soil extraction - if(input.tag == "soil" && is.null(input$path)){ - settings$run$inputs[[i]]$path <- PEcAn.data.land::soil_process(settings, input, dbfiles.local, overwrite=FALSE) + if (input.tag == "soil" && is.null(input$path)) { + settings$run$inputs[[i]]$path <- PEcAn.data.land::soil_process(settings, input, dbfiles.local, overwrite = FALSE) needsave <- TRUE ## NOTES: at the moment only processing soil locally. Need to think about how to generalize this ## because many models will read PEcAn standard in write.configs and write out into settings @@ -75,30 +75,38 @@ do_conversions <- function(settings, overwrite.met = FALSE, overwrite.fia = FALS } # Phenology data extraction - if(input.tag == "leaf_phenology" && is.null(input$path)){ - #settings$run$inputs[[i]]$path <- PEcAn.data.remote::extract_phenology_MODIS(site_info,start_date,end_date,outdir,run_parallel = TRUE,ncores = NULL) + if (input.tag == "leaf_phenology" && is.null(input$path)) { + # settings$run$inputs[[i]]$path <- PEcAn.data.remote::extract_phenology_MODIS( + # site_info, + # start_date, + # end_date, + # outdir, + # run_parallel = TRUE, + # ncores = NULL + # ) needsave <- TRUE } # met conversion - + if (input.tag == "met") { name <- "MET Process" - if ( (PEcAn.utils::status.check(name) == 0)) { ## previously is.null(input$path) && - PEcAn.logger::logger.info("calling met.process: ",settings$run$inputs[[i]][['path']]) - settings$run$inputs[[i]] <- + if ((PEcAn.utils::status.check(name) == 0)) { ## previously is.null(input$path) && + PEcAn.logger::logger.info("calling met.process: ", settings$run$inputs[[i]][["path"]]) + settings$run$inputs[[i]] <- PEcAn.data.atmosphere::met.process( - site = settings$run$site, + site = settings$run$site, input_met = settings$run$inputs$met, start_date = settings$run$start.date, end_date = settings$run$end.date, model = settings$model$type, host = settings$host, - dbparms = settings$database$bety, + dbparms = settings$database$bety, dir = dbfiles, spin = settings$spin, - overwrite = overwrite.met) - PEcAn.logger::logger.debug("updated met path: ",settings$run$inputs[[i]][['path']]) + overwrite = overwrite.met + ) + PEcAn.logger::logger.debug("updated met path: ", settings$run$inputs[[i]][["path"]]) needsave <- TRUE } } diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R index 22ed6f3b729..f8036400989 100644 --- a/base/workflow/R/run.write.configs.R +++ b/base/workflow/R/run.write.configs.R @@ -22,51 +22,58 @@ #' @export #' #' @author David LeBauer, Shawn Serbin, Ryan Kelly, Mike Dietze -run.write.configs <- function(settings, write = TRUE, ens.sample.method = "uniform", - posterior.files = rep(NA, length(settings$pfts)), +run.write.configs <- function(settings, write = TRUE, ens.sample.method = "uniform", + posterior.files = rep(NA, length(settings$pfts)), overwrite = TRUE) { - tryCatch({ - con <- PEcAn.DB::db.open(settings$database$bety) - on.exit(PEcAn.DB::db.close(con), add = TRUE) - }, error = function(e) { - PEcAn.logger::logger.severe( - "Connection requested, but failed to open with the following error: ", - conditionMessage(e)) - }) - + tryCatch( + { + con <- PEcAn.DB::db.open(settings$database$bety) + on.exit(PEcAn.DB::db.close(con), add = TRUE) + }, + error = function(e) { + PEcAn.logger::logger.severe( + "Connection requested, but failed to open with the following error: ", + conditionMessage(e) + ) + } + ) + ## Which posterior to use? for (i in seq_along(settings$pfts)) { ## if posterior.files is specified us that if (is.na(posterior.files[i])) { ## otherwise, check to see if posteriorid exists if (!is.null(settings$pfts[[i]]$posteriorid)) { - #TODO: sometimes `files` is a 0x0 tibble and other operations with it fail. + # TODO: sometimes `files` is a 0x0 tibble and other operations with it fail. files <- PEcAn.DB::dbfile.check("Posterior", - settings$pfts[[i]]$posteriorid, - con, settings$host$name, return.all = TRUE) - pid <- grep("post.distns.*Rdata", files$file_name) ## is there a posterior file? + settings$pfts[[i]]$posteriorid, + con, settings$host$name, + return.all = TRUE + ) + pid <- grep("post.distns.*Rdata", files$file_name) ## is there a posterior file? if (length(pid) == 0) { - pid <- grep("prior.distns.Rdata", files$file_name) ## is there a prior file? + pid <- grep("prior.distns.Rdata", files$file_name) ## is there a prior file? } if (length(pid) > 0) { posterior.files[i] <- file.path(files$file_path[pid], files$file_name[pid]) - } ## otherwise leave posteriors as NA + } ## otherwise leave posteriors as NA } ## otherwise leave NA and get.parameter.samples will look for local } else { ## does posterior.files point to a directory instead of a file? - if(utils::file_test("-d",posterior.files[i])){ - pfiles = dir(posterior.files[i],pattern = "post.distns.*Rdata",full.names = TRUE) - if(length(pfiles)>1){ - pid = grep("post.distns.Rdata",pfiles) - if(length(pid > 0)){ - pfiles = pfiles[grep("post.distns.Rdata",pfiles)] + if (utils::file_test("-d", posterior.files[i])) { + pfiles <- dir(posterior.files[i], pattern = "post.distns.*Rdata", full.names = TRUE) + if (length(pfiles) > 1) { + pid <- grep("post.distns.Rdata", pfiles) + if (length(pid > 0)) { + pfiles <- pfiles[grep("post.distns.Rdata", pfiles)] } else { PEcAn.logger::logger.error( "run.write.configs: could uniquely identify posterior files within", - posterior.files[i]) + posterior.files[i] + ) } - posterior.files[i] = pfiles + posterior.files[i] <- pfiles } } ## also, double check PFT outdir exists @@ -74,9 +81,9 @@ run.write.configs <- function(settings, write = TRUE, ens.sample.method = "unifo ## no outdir settings$pfts[[i]]$outdir <- file.path(settings$outdir, "pfts", settings$pfts[[i]]$name) } - } ## end else + } ## end else } ## end for loop over pfts - + ## Sample parameters model <- settings$model$type scipen <- getOption("scipen") @@ -95,90 +102,96 @@ run.write.configs <- function(settings, write = TRUE, ens.sample.method = "unifo } else { PEcAn.logger::logger.error(samples.file, "not found, this file is required by the run.write.configs function") } - + ## remove previous runs.txt if (overwrite && file.exists(file.path(settings$rundir, "runs.txt"))) { PEcAn.logger::logger.warn("Existing runs.txt file will be removed.") unlink(file.path(settings$rundir, "runs.txt")) } - + PEcAn.utils::load.modelpkg(model) - + ## Check for model-specific write configs - - my.write.config <- paste0("write.config.",model) + + my.write.config <- paste0("write.config.", model) if (!exists(my.write.config)) { - PEcAn.logger::logger.error(my.write.config, - "does not exist, please make sure that the model package contains a function called", - my.write.config) + PEcAn.logger::logger.error( + my.write.config, + "does not exist, please make sure that the model package contains a function called", + my.write.config + ) } - + ## Prepare for model output. Clean up any old config files (if exists) - #TODO: shouldn't this check if the files exist before removing them? + # TODO: shouldn't this check if the files exist before removing them? my.remove.config <- paste0("remove.config.", model) if (exists(my.remove.config)) { do.call(my.remove.config, args = list(settings$rundir, settings)) } - + # TODO RK : need to write to runs_inputs table - + # Save names pft.names <- names(trait.samples) trait.names <- lapply(trait.samples, names) - + ### NEED TO IMPLEMENT: Load Environmental Priors and Posteriors - + ### Sensitivity Analysis if ("sensitivity.analysis" %in% names(settings)) { - ### Write out SA config files PEcAn.logger::logger.info("\n ----- Writing model run config files ----") - sa.runs <- PEcAn.uncertainty::write.sa.configs(defaults = settings$pfts, - quantile.samples = sa.samples, - settings = settings, - model = model, - write.to.db = write) - + sa.runs <- PEcAn.uncertainty::write.sa.configs( + defaults = settings$pfts, + quantile.samples = sa.samples, + settings = settings, + model = model, + write.to.db = write + ) + # Store output in settings and output variables runs.samples$sa <- sa.run.ids <- sa.runs$runs settings$sensitivity.analysis$ensemble.id <- sa.ensemble.id <- sa.runs$ensemble.id - + # Save sensitivity analysis info fname <- PEcAn.uncertainty::sensitivity.filename(settings, "sensitivity.samples", "Rdata", - all.var.yr = TRUE, pft = NULL) + all.var.yr = TRUE, pft = NULL + ) save(sa.run.ids, sa.ensemble.id, sa.samples, pft.names, trait.names, file = fname) - - } ### End of SA - + } ### End of SA + ### Write ENSEMBLE if ("ensemble" %in% names(settings)) { - ens.runs <- PEcAn.uncertainty::write.ensemble.configs(defaults = settings$pfts, - ensemble.samples = ensemble.samples, - settings = settings, - model = model, - write.to.db = write) - + ens.runs <- PEcAn.uncertainty::write.ensemble.configs( + defaults = settings$pfts, + ensemble.samples = ensemble.samples, + settings = settings, + model = model, + write.to.db = write + ) + # Store output in settings and output variables runs.samples$ensemble <- ens.run.ids <- ens.runs$runs settings$ensemble$ensemble.id <- ens.ensemble.id <- ens.runs$ensemble.id - ens.samples <- ensemble.samples # rename just for consistency - + ens.samples <- ensemble.samples # rename just for consistency + # Save ensemble analysis info fname <- PEcAn.uncertainty::ensemble.filename(settings, "ensemble.samples", "Rdata", all.var.yr = TRUE) save(ens.run.ids, ens.ensemble.id, ens.samples, pft.names, trait.names, file = fname) } else { PEcAn.logger::logger.info("not writing config files for ensemble, settings are NULL") - } ### End of Ensemble - + } ### End of Ensemble + PEcAn.logger::logger.info("###### Finished writing model run config files #####") PEcAn.logger::logger.info("config files samples in ", file.path(settings$outdir, "run")) - + ### Save output from SA/Ensemble runs - # A lot of this is duplicate with the ensemble/sa specific output above, but kept for backwards compatibility. - save(ensemble.samples, trait.samples, sa.samples, runs.samples, pft.names, trait.names, - file = file.path(settings$outdir, "samples.Rdata")) + # A lot of this is duplicate with the ensemble/sa specific output above, but kept for backwards compatibility. + save(ensemble.samples, trait.samples, sa.samples, runs.samples, pft.names, trait.names, + file = file.path(settings$outdir, "samples.Rdata") + ) PEcAn.logger::logger.info("parameter values for runs in ", file.path(settings$outdir, "samples.RData")) options(scipen = scipen) - + return(invisible(settings)) } diff --git a/base/workflow/R/runModule.get.trait.data.R b/base/workflow/R/runModule.get.trait.data.R index ce084567714..873163c0d0e 100644 --- a/base/workflow/R/runModule.get.trait.data.R +++ b/base/workflow/R/runModule.get.trait.data.R @@ -4,7 +4,9 @@ ##' `MultiSettings` ##' @export runModule.get.trait.data <- function(settings) { - if (is.null(settings$meta.analysis)) return(settings) ## if there's no MA, there's no need for traits + if (is.null(settings$meta.analysis)) { + return(settings) # if there's no MA, there's no need for traits + } if (PEcAn.settings::is.MultiSettings(settings)) { pfts <- list() pft.names <- character(0) @@ -36,9 +38,11 @@ runModule.get.trait.data <- function(settings) { modeltype = settings$model$type, dbfiles = settings$database$dbfiles, database = settings$database$bety, - forceupdate = ifelse(is.null(settings$meta.analysis$update), - FALSE, - settings$meta.analysis$update), + forceupdate = ifelse( + is.null(settings$meta.analysis$update), + FALSE, + settings$meta.analysis$update + ), write = settings$database$bety$write ) diff --git a/base/workflow/R/runModule.run.write.configs.R b/base/workflow/R/runModule.run.write.configs.R index 61b7d0bed28..ebb992130e6 100644 --- a/base/workflow/R/runModule.run.write.configs.R +++ b/base/workflow/R/runModule.run.write.configs.R @@ -5,7 +5,6 @@ #' @return A modified settings object, invisibly #' @export runModule.run.write.configs <- function(settings, overwrite = TRUE) { - if (PEcAn.settings::is.MultiSettings(settings)) { if (overwrite && file.exists(file.path(settings$rundir, "runs.txt"))) { PEcAn.logger::logger.warn("Existing runs.txt file will be removed.") @@ -32,7 +31,8 @@ runModule.run.write.configs <- function(settings, overwrite = TRUE) { write = settings$database$bety$write, ens.sample.method = settings$ensemble$samplingspace$parameters$method, posterior.files = posterior.files, - overwrite = overwrite)) + overwrite = overwrite + )) } else { stop("runModule.run.write.configs only works with Settings or MultiSettings") } diff --git a/base/workflow/R/start_model_runs.R b/base/workflow/R/start_model_runs.R index 12ef228704e..e4ec46fbbbc 100644 --- a/base/workflow/R/start_model_runs.R +++ b/base/workflow/R/start_model_runs.R @@ -6,42 +6,45 @@ #' @export #' @examples #' \dontrun{ -#' start_model_runs(settings) +#' start_model_runs(settings) #' } #' @author Shawn Serbin, Rob Kooper, David LeBauer, Alexey Shiklomanov #' start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { - run_file <- file.path(settings$rundir, "runs.txt") # check if runs need to be done if (!file.exists(run_file)) { PEcAn.logger::logger.warn( - "runs.txt not found, assuming no runs need to be done") + "runs.txt not found, assuming no runs need to be done" + ) return() } run_list <- readLines(con = run_file) nruns <- length(run_list) if (nruns == 0) { PEcAn.logger::logger.warn( - "runs.txt found, but is empty. Assuming no runs need to be done") + "runs.txt found, but is empty. Assuming no runs need to be done" + ) return() } - + model <- settings$model$type PEcAn.logger::logger.info( - "-------------------------------------------------------------------") + "-------------------------------------------------------------------" + ) PEcAn.logger::logger.info(paste(" Starting model runs", model)) PEcAn.logger::logger.info( - "-------------------------------------------------------------------") - + "-------------------------------------------------------------------" + ) + is_local <- PEcAn.remote::is.localhost(settings$host) is_qsub <- !is.null(settings$host$qsub) is_rabbitmq <- !is.null(settings$host$rabbitmq) is_modellauncher <- !is.null(settings$host$modellauncher) - + # Check if Njobmax tag exists in seetings - if (is_modellauncher){ - if (!is.null(settings$host$modellauncher$Njobmax)){ + if (is_modellauncher) { + if (!is.null(settings$host$modellauncher$Njobmax)) { Njobmax <- settings$host$modellauncher$Njobmax } else { Njobmax <- nruns @@ -50,14 +53,14 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { compt_run_modellauncher <- 1 job_modellauncher <- list() } - + # loop through runs and either call start run, or launch job on remote machine jobids <- list() - + ## setup progressbar pb <- utils::txtProgressBar(min = 0, max = nruns, style = 3) pbi <- 0 - + # create database connection if (write) { dbcon <- PEcAn.DB::db.open(settings$database$bety) @@ -65,53 +68,53 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { } else { dbcon <- NULL } - + # launcher folder jobfile <- NULL firstrun <- NULL - - #Copy all run directories over if not local + + # Copy all run directories over if not local if (!is_local) { # copy over run directories PEcAn.utils::retry.func( PEcAn.remote::remote.copy.to( host = settings$host, - src = settings$rundir, - dst = dirname(settings$host$rundir), + src = settings$rundir, + dst = dirname(settings$host$rundir), delete = TRUE - ), + ), sleep = 2 ) - + # copy over out directories PEcAn.utils::retry.func( PEcAn.remote::remote.copy.to( host = settings$host, src = settings$modeloutdir, dst = dirname(settings$host$outdir), - #include all directories, exclude all files + # include all directories, exclude all files options = c("--include=*/", "--exclude=*"), delete = TRUE ), sleep = 2 ) } - + # launch each of the jobs for (run in run_list) { run_id_string <- format(run, scientific = FALSE) # write start time to database PEcAn.DB::stamp_started(con = dbcon, run = run) - + # check to see if we use the model launcher if (is_rabbitmq) { run_id_string <- format(run, scientific = FALSE) folder <- file.path(settings$rundir, run_id_string) out <- PEcAn.remote::start_rabbitmq( - folder, settings$host$rabbitmq$uri, settings$host$rabbitmq$queue) + folder, settings$host$rabbitmq$uri, settings$host$rabbitmq$queue + ) PEcAn.logger::logger.debug("JOB.SH submit status:", out) jobids[run] <- folder - } else if (is_modellauncher) { # set up launcher script if we use modellauncher if (is.null(firstrun)) { @@ -121,15 +124,16 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { rundir = settings$rundir, host_rundir = settings$host$rundir, mpirun = settings$host$modellauncher$mpirun, - binary = settings$host$modellauncher$binary) + binary = settings$host$modellauncher$binary + ) job_modellauncher[compt_run_modellauncher] <- run - compt_run_modellauncher <- compt_run_modellauncher+1 + compt_run_modellauncher <- compt_run_modellauncher + 1 } writeLines( c(file.path(settings$host$rundir, run_id_string)), - con = jobfile) + con = jobfile + ) pbi <- pbi + 1 - } else if (is_qsub) { out <- PEcAn.remote::start_qsub( run = run, @@ -140,13 +144,14 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { host_outdir = settings$host$outdir, stdout_log = "stdout.log", stderr_log = "stderr.log", - job_script = "job.sh") + job_script = "job.sh" + ) PEcAn.logger::logger.debug("JOB.SH submit status:", out) jobids[run] <- PEcAn.remote::qsub_get_jobid( out = out[length(out)], qsub.jobid = settings$host$qsub.jobid, - stop.on.error = stop.on.error) - + stop.on.error = stop.on.error + ) } else { # if qsub option is not invoked. just start model runs in serial. out <- PEcAn.remote::start_serial( @@ -154,53 +159,53 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { host = settings$host, rundir = settings$rundir, host_rundir = settings$host$rundir, - job_script = "job.sh") - + job_script = "job.sh" + ) + # check output to see if an error occurred during the model run PEcAn.remote::check_model_run(out = out, stop.on.error = stop.on.error) - + if (!is_local) { # copy data back to local PEcAn.utils::retry.func( PEcAn.remote::remote.copy.from( host = settings$host, src = file.path(settings$host$outdir, run_id_string), - dst = settings$modeloutdir), + dst = settings$modeloutdir + ), sleep = 2 ) } - + # write finished time to database PEcAn.DB::stamp_finished(con = dbcon, run = run) - + pbi <- pbi + 1 utils::setTxtProgressBar(pb, pbi) } - + # Check if compt_run has reached Njobmax - if (is_modellauncher){ + if (is_modellauncher) { compt_run <- compt_run + 1 - if (compt_run == Njobmax){ + if (compt_run == Njobmax) { close(jobfile) firstrun <- NULL compt_run <- 0 jobfile <- NULL - } + } } - } # end loop over runs close(pb) - + # need to actually launch the model launcher if (is_modellauncher) { - # Only close if not already closed - if (compt_run != 0){ + if (compt_run != 0) { close(jobfile) } - + if (!is_local) { - for (run in run_list){ #only re-copy run dirs that have launcher and job list + for (run in run_list) { # only re-copy run dirs that have launcher and job list if (run %in% job_modellauncher) { # copy launcher and joblist PEcAn.utils::retry.func( @@ -208,15 +213,15 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { host = settings$host, src = file.path(settings$rundir, format(run, scientific = FALSE)), dst = settings$host$rundir, - delete = TRUE), + delete = TRUE + ), sleep = 2 ) - } } } if (is_qsub) { - for (run in run_list){ + for (run in run_list) { if (run %in% job_modellauncher) { out <- PEcAn.remote::start_qsub( run = run, @@ -228,71 +233,73 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { stdout_log = "launcher.out.log", stderr_log = "launcher.err.log", job_script = "launcher.sh", - qsub_extra = settings$host$modellauncher$qsub) + qsub_extra = settings$host$modellauncher$qsub + ) } # HACK: Code below gets 'run' from names(jobids) so need an entry for # each run. But when using modellauncher all runs have the same jobid jobids[run] <- sub(settings$host$qsub.jobid, "\\1", out[length(out)]) } - } else { out <- PEcAn.remote::start_serial( run = run, host = settings$host, rundir = settings$rundir, host_rundir = settings$host$rundir, - job_script = "launcher.sh") - + job_script = "launcher.sh" + ) + # check output to see if an error occurred during the model run PEcAn.remote::check_model_run(out = out, stop.on.error = TRUE) - + # write finished time to database for (run in run_list) { PEcAn.DB::stamp_finished(con = dbcon, run = run) } - + utils::setTxtProgressBar(pb, pbi) } } - + # wait for all jobs to finish if (length(jobids) > 0) { PEcAn.logger::logger.debug( "Waiting for the following jobs:", - unlist(unique(jobids))) + unlist(unique(jobids)) + ) } - - #TODO figure out a way to do this while for unique(jobids) instead of jobids + + # TODO figure out a way to do this while for unique(jobids) instead of jobids while (length(jobids) > 0) { Sys.sleep(10) - + if (!is_local) { - #Copy over log files to check progress + # Copy over log files to check progress try(PEcAn.remote::remote.copy.from( host = settings$host, src = settings$host$outdir, dst = dirname(settings$modeloutdir), - options = c('--exclude=*.h5') + options = c("--exclude=*.h5") )) } - + for (run in names(jobids)) { run_id_string <- format(run, scientific = FALSE) - + # check to see if job is done job_finished <- FALSE if (is_rabbitmq) { - job_finished <- + job_finished <- file.exists(file.path(jobids[run], "rabbitmq.out")) } else if (is_qsub) { job_finished <- PEcAn.remote::qsub_run_finished( run = jobids[run], host = settings$host, - qstat = settings$host$qstat) + qstat = settings$host$qstat + ) } - + if (job_finished) { - # TODO check output log if (is_rabbitmq) { data <- readLines(file.path(jobids[run], "rabbitmq.out")) @@ -305,9 +312,9 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { } } } - + # Write finish time to database - #TODO this repeats for every run in `jobids` writing every run's time stamp every time. This actually takes quite a long time with a lot of ensembles and should either 1) not be a for loop (no `for(x in run_list)`) or 2) if `is_modellauncher`, be done outside of the jobids for loop after all jobs are finished. + # TODO this repeats for every run in `jobids` writing every run's time stamp every time. This actually takes quite a long time with a lot of ensembles and should either 1) not be a for loop (no `for(x in run_list)`) or 2) if `is_modellauncher`, be done outside of the jobids for loop after all jobs are finished. if (is_modellauncher) { for (x in run_list) { PEcAn.DB::stamp_finished(con = dbcon, run = x) @@ -315,25 +322,25 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { } else { PEcAn.DB::stamp_finished(con = dbcon, run = run) } - + # move progress bar if (!is_modellauncher) { pbi <- pbi + 1 } utils::setTxtProgressBar(pb, pbi) - + # remove job if (is_modellauncher) { for (x in run_list) { jobids[x] <- NULL - } + } } else { jobids[run] <- NULL } } # End job finished - } # end loop over runs - } # end while loop checking runs - + } # end loop over runs + } # end while loop checking runs + # Copy data back to local if (!is_local) { PEcAn.utils::retry.func( @@ -354,13 +361,14 @@ start_model_runs <- function(settings, write = TRUE, stop.on.error = TRUE) { #' `settings` instead of as a separate argument. #' #' @export -runModule_start_model_runs <- function(settings, stop.on.error=TRUE) { +runModule_start_model_runs <- function(settings, stop.on.error = TRUE) { if (PEcAn.settings::is.MultiSettings(settings) || - PEcAn.settings::is.Settings(settings)) { + PEcAn.settings::is.Settings(settings)) { write <- settings$database$bety$write return(start_model_runs(settings, write, stop.on.error)) } else { PEcAn.logger::logger.severe( - "runModule_start_model_runs only works with Settings or MultiSettings") + "runModule_start_model_runs only works with Settings or MultiSettings" + ) } -} # runModule_start_model_runs \ No newline at end of file +} # runModule_start_model_runs diff --git a/base/workflow/man/start_model_runs.Rd b/base/workflow/man/start_model_runs.Rd index a1c171073bb..f6e34703150 100644 --- a/base/workflow/man/start_model_runs.Rd +++ b/base/workflow/man/start_model_runs.Rd @@ -27,7 +27,7 @@ Start selected ecosystem model runs within PEcAn workflow }} \examples{ \dontrun{ - start_model_runs(settings) +start_model_runs(settings) } } \author{ diff --git a/base/workflow/tests/testthat/test-runModule.get.trait.data.R b/base/workflow/tests/testthat/test-runModule.get.trait.data.R index e1ba65a12a1..7293b64b035 100644 --- a/base/workflow/tests/testthat/test-runModule.get.trait.data.R +++ b/base/workflow/tests/testthat/test-runModule.get.trait.data.R @@ -1,11 +1,11 @@ - context("testing workflow functions") settings <- PEcAn.settings::Settings(list( outdir = "/dev/null", database = list(bety = list(), dbfiles = "fake_path"), pfts = list(pft = list(name = "fake", outdir = "fake_pft_path")), - meta.analysis = list(threshold = 1, iter = 1))) + meta.analysis = list(threshold = 1, iter = 1) +)) test_that("settings not changed if no MA", { settings_noMA <- settings @@ -14,11 +14,12 @@ test_that("settings not changed if no MA", { }) test_that("get.trait.data is called exactly once", { - mock_vals <- mockery::mock(settings$pfts, cycle=TRUE) + mock_vals <- mockery::mock(settings$pfts, cycle = TRUE) mockery::stub( where = runModule.get.trait.data, what = "PEcAn.DB::get.trait.data", - how = mock_vals) + how = mock_vals + ) res <- runModule.get.trait.data(settings) expect_equal(res, settings) @@ -30,4 +31,4 @@ test_that("get.trait.data is called exactly once", { expect_equal(length(res_multi), 3) expect_identical(res_multi[[1]], settings) mockery::expect_called(mock_vals, 2) -}) \ No newline at end of file +}) diff --git a/base/workflow/tests/testthat/test.do_conversions.R b/base/workflow/tests/testthat/test.do_conversions.R index a28a09667ff..8ba03f66ad1 100644 --- a/base/workflow/tests/testthat/test.do_conversions.R +++ b/base/workflow/tests/testthat/test.do_conversions.R @@ -1,6 +1,6 @@ test_that("`do_conversions` able to return settings from pecan.METProcess.xml if it already exists", { withr::with_tempdir({ - settings <- list(host = list(name = 'test', folder = 'test'), outdir = getwd()) + settings <- list(host = list(name = "test", folder = "test"), outdir = getwd()) file_path <- file.path(getwd(), "pecan.METProcess.xml") file.create(file_path) writeLines( @@ -16,16 +16,16 @@ test_that("`do_conversions` able to return settings from pecan.METProcess.xml if test_that("`do_conversions` able to call met.process if the input tag has met, update the met path and save settings to pecan.METProcess.xml", { withr::with_tempdir({ - mocked_res <- mockery::mock(list(path = 'test')) - mockery::stub(do_conversions, 'PEcAn.data.atmosphere::met.process', mocked_res) + mocked_res <- mockery::mock(list(path = "test")) + mockery::stub(do_conversions, "PEcAn.data.atmosphere::met.process", mocked_res) settings <- list( - host = list(name = 'test', folder = 'test'), + host = list(name = "test", folder = "test"), outdir = getwd(), run = list(site = list(id = 0), inputs = list(met = list(id = 1))) ) res <- do_conversions(settings) mockery::expect_called(mocked_res, 1) - expect_equal(res$run$inputs$met$path, 'test') + expect_equal(res$run$inputs$met$path, "test") expect_true(file.exists(file.path(getwd(), "pecan.METProcess.xml"))) }) }) diff --git a/base/workflow/tests/testthat/test.start_model_runs.R b/base/workflow/tests/testthat/test.start_model_runs.R index 052950e2c93..7f876a0431c 100644 --- a/base/workflow/tests/testthat/test.start_model_runs.R +++ b/base/workflow/tests/testthat/test.start_model_runs.R @@ -17,4 +17,3 @@ test_that("`start_model_runs` throws a warning if runs.txt is empty", { expect_output(start_model_runs(settings), "runs.txt found, but is empty") }) }) -