Add functions that can extract and update phenological status based on MODIS data and enable monthly SDA #3249
Conversation
Merge Dongchen's branch
…velop Merge from the upstream.
Qianyuxuan
force-pushed
the
sda_monthly_modex
branch
from
a0ada30
to
9136c43
Compare
…sda_monthly_local Merge from the upstream
|
@Qianyuxuan there's a bunch of errors coming up in the build that suggest that you're missing namespaces on functions |
|
@mdietze |
|
|
| prior.distns <- post.distns | ||
| } | ||
| if (exists("trait.mcmc") & exists("post.distns")) { | ||
| post.distns <- post.distns[!(row.names(post.distns) %in% names(trait.mcmc)),] |
There was a problem hiding this comment.
This appears to be undoing recent changes from #3220 -- is that intentional?
There was a problem hiding this comment.
This function was adapted from Alexis Helgeson's codes to read "basal soil respiration rate" from "trait.mcmc.Rdata" as we found basal soil respiration rate seems too high from "post.distns.Rdata" for NEON sites. I am not sure whether she put forward a PR for this change. But maybe we can revert it for the develop branch as the changes are specific for our study?
There was a problem hiding this comment.
@Qianyuxuan so @infotroph isn't objecting to using trait.mcmc or post.distns, but that when you resolved you merge conflict you undid the effort to load the files into a specific envir and then use that envir name (distns). I think you can easily adapt your code to include the previous changes.
There was a problem hiding this comment.
Yes, correct -- loading files into a specific envir is the important part.
Also, to make sure I'm reading right: Am I correct that all the changes in this file are specific to basal soil respiration, they're simply making sure trait.mcmc and post.distns are available to be processed elsewhere?
There was a problem hiding this comment.
@Qianyuxuan you still need to fix this change before this PR can go in.
| #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. |
There was a problem hiding this comment.
This seems like a reasonable fix for this PR. Long term it might be worth being consistent about calling this time_bnds throughout, which I think is better aligned with the CF standard: ORNL DAAC lists time_bnds as required, and time_bnds appears* in the conventions themseves whereas time_bounds does not. Can we consider a followup patch to make read_restart use time_bnds?
*Though admittedly time_bnds only appears in an non-normative example. The CF conventions are messy!
Co-authored-by: Chris Black <chris@ckblack.org>
Co-authored-by: Chris Black <chris@ckblack.org>
Co-authored-by: Chris Black <chris@ckblack.org>
Co-authored-by: Chris Black <chris@ckblack.org>
Co-authored-by: Chris Black <chris@ckblack.org>
Co-authored-by: Chris Black <chris@ckblack.org>
| @@ -240,8 +240,9 @@ template <- PEcAn.settings::Settings(list( | |||
| type = "SIPNET", | |||
| revision = "ssr", | |||
| delete.raw = FALSE, | |||
| binary = model_binary, | |||
| jobtemplate = "~/sipnet_geo.job" | |||
| binary = "/usr2/postdoc/istfer/SIPNET/trunk//sipnet_if", | |||
There was a problem hiding this comment.
is there a reason the variable model_binary got replaced with a specific path rather than just setting the model_binary variable to the path you wanted?
| @@ -77,11 +78,18 @@ prepare_pools <- function(nc.path, constants = NULL){ | |||
| PEcAn.logger::logger.error("TotLivBiom is less than sum of AbvGrndWood and roots; will use default for leaf biomass") | |||
| } | |||
| } | |||
| #Initial LAI at the beginning of year is set as 0 for deciduous forests and grasslands | |||
There was a problem hiding this comment.
This solution is model specific and the approach you're using to implement it will only work in a sitegroup run (pft.site doesn't exist for single-site runs). However, prepare_pools is a general, across-model function. Furthermore, filtering for one specific PFT name doesn't work, even within SIPNET, as there are other evergreen PFTs in the world other than boreal.coniferous. I'd recommend moving this hack to write.configs.SIPNET where you can more reliably detect the PFT being used, the start date of the run, and the phenological state of the system
| ##' @param constants list of constants; must include SLA in m2 / kg C if providing LAI for leaf carbon | ||
| ##' @return list of pool values in kg C / m2 with generic names | ||
| ##' @author Anne Thomas | ||
| prepare_pools <- function(nc.path, constants = NULL){ | ||
| prepare_pools <- function(nc.path, settings, constants = NULL){ |
There was a problem hiding this comment.
per comment below about relocating your LAI hack, please remove settings as an argument and remake the documentation
| @@ -198,6 +198,8 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples, | |||
| ##' @param write.to.db logical: Record this run in BETY? | |||
| ##' @param restart In case this is a continuation of an old simulation. restart needs to be a list with name tags of runid, inputs, new.params (parameters), new.state (initial condition), ensemble.id (ensemble id), start.time and stop.time.See Details. | |||
| ##' @param rename Decide if we want to rename previous output files, for example convert from sipnet.out to sipnet.2020-07-16.out. | |||
| ##' @param time obervation timepoints | |||
There was a problem hiding this comment.
write.ensemble.configs is not a SDA-specific function in PEcAn, but one that is used throughout the code to write configs. The two variables you've added aren't general, and aren't documented in a way that would make sense to a general user. I'm also surprised that they can't be passed in as part of the settings. please update
| RENAME = rename)#for restart from previous model runs, not sharing the same outdir | ||
| RENAME = rename, | ||
| obs_time = time, | ||
| update_phenology = update_phenology)#for restart from previous model runs, not sharing the same outdir |
There was a problem hiding this comment.
per previous comment, please remove obs_time and update_phenology, otherwise you will have also update every single model's write.restart function too. Also, please don't relocate a comment about one variable (RENAME) to now point to a different variable (update_phenology), as this is going to confuse the heck out of users.
General recommendation -- before you add/remove function arguments, please think about what functions are model-specific and analysis-specific (e.g. SDA, PDA, uncertainty analysis) versus those that are used generally throughout the system.
| ##' @param prefix prefix to read the output files | ||
| ##' @param delete.raw Flag to remove sipnet.out files, FALSE = do not remove files TRUE = remove files | ||
| ##' | ||
| ##' @export | ||
| ##' @author Shawn Serbin, Michael Dietze | ||
| model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out", | ||
| overwrite = FALSE, conflict = FALSE) { | ||
| overwrite = FALSE, conflict = TRUE) { |
There was a problem hiding this comment.
Does it make sense to change the default here rather than just changing the value being passed in when your running the sub-annual SDA specifically?
| } ## end loop over PFTS | ||
|
|
||
| #update LeafOnday and LeafOffDay | ||
| if (update_phenology){ |
There was a problem hiding this comment.
Is there a reason that update_phenology needs to be a global variable rather than being a setting in settings$model$leaf_phenology or something similar?
|
|
||
| #update LeafOnday and LeafOffDay | ||
| if (update_phenology){ | ||
| leaf_pheno_outdir <- settings$model$leaf_phenology$outdir ## read from settings |
There was a problem hiding this comment.
- the leaf phenology data is an INPUT and thus should be referenced like all other inputs, e.g. settings$run$inputs$leaf_phenology.
- As such leaf_phenology should be registered with BETY as a new optional input to SIPNET.
- The SIPNET run is not writing anything new to leaf_phenology$outdir and thus it shouldn't be called
outdirbut should useid,path, orsourcelike all other inputs, depending on whether you're specifying a database ID, a file path, or the name of the input source that is then processed within do.conversions. - It's worth noting that prescribed phenology is not a SIPNET-specific feature, but something that could be valuable to any model, which is an additional argument for moving it to inputs instead of model
| ##read data | ||
| leafphdata <- utils::read.csv(paste0(leaf_pheno_outdir,"/leaf_phenology.csv")) | ||
| leafOnDay <- leafphdata$leafonday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id] | ||
| leafOffDay<- leafphdata$leafoffday[leafphdata$year==obs_time & leafphdata$site_id==settings$run$site$id] |
There was a problem hiding this comment.
- why does this have to be obs_time instead of start.date?
- I would expect obs_time to be in a datetime format, and yet you're matching it against year? Seems like they're a missing extraction of year
- how does SDA update the phenological status if the phenology is prescribed?
- If you're adding a whole new class of inputs to the PEcAn system (prescribed phenology) there needs to be updates to the general documentation explaining any new tags, new sources, new file formats, etc. Indeed, there needs to be updates explaining how to run the sub-annual SDA more generally as well.
| ##' @author Qianyu Li | ||
| ##' | ||
|
|
||
| phenology_MODIS_extract<- function(settings,run_parallel = TRUE,ncores = NULL){ |
There was a problem hiding this comment.
- since this is an input processing function I'd recommend tweaking the function name and arguments to be more inline with existing input processing functions (e.g. starting with extract or download)
- function doesn't seem to be called anywhere from within the existing workflow. I suggest adding to
do.conversions - in keeping with the current philosophy of making PEcAn modules more independent, I'd recommend shifting the site_info bit into do.conversions itself and making the arguments to this function be (outdir, site_info, start_date, end_date, run_parallel, ncores). In other words, don't pass in the overall
settingsobject, just pass the info you need. - Since this functionality is generally useful, not just in SDA, I'd recommend that the code you move to do.conversion grab start_date and end_date from the overall settings, not settings$state.data.assimilation
Description
Motivation and Context
Review Time Estimate
Types of changes
Checklist: