Add codes to extract SoilGrids soil texture data and derive ensemble … #3406
Conversation
…soil parameter files
| #working on reading soil file (only working for soilWHC) | ||
| if (length(settings$run$inputs$soilinitcond$path) > 0) { | ||
| soil_ICs_num <- length(settings$run$inputs$soilinitcond$path) | ||
| soil_ICs_path <- settings$run$inputs$soilinitcond$path[[sample(1:soil_ICs_num, 1)]] |
There was a problem hiding this comment.
I'm guessing the logic of this bit is that if an ensemble input is passed in then you choose one at random? Personally, I'd prefer to throw an error at this point, as by this point in the workflow ensemble inputs should have already been converted to specific inputs for each run. If they're still ensemble at this point it's an error and we should stop. Proceeding on with a random choice means (a) we don't catch an error and (b) we don't know/record which ensemble member mapped to which run.
There was a problem hiding this comment.
Yeah I followed the logic as choosing one from ensemble initial conditions here:https://github.com/PecanProject/pecan/blob/develop/models/sipnet/R/write.configs.SIPNET.R#L546. So you suggested we should just randomly choose and write one ensemble member in the script "soil_params_ensemble.R" and read that specific path in "write.configs.SIPNET.R"?
There was a problem hiding this comment.
So your soil_params_ensemble should generate a bunch of ensembles and load them ALL into the overall settings object. What I'm saying is that write.configs.SIPNET should only end up seeing ONE of those ensembles because run.write.configs should do the joint sampling of all the different types of ensembles (met, soils, IC, etc), record which ensemble member was assigned which inputs, and then ensure that the settings object reaching write.configs.SIPNET is only for a specific ensemble member. If you're finding that write.configs.SIPNET is being handed all of the soil ensemble members then there's something wrong upstream in the pecan.xml (e.g. not including soils in the section) or in run.write.configs
There was a problem hiding this comment.
@mdietze It seems that "run.write.configs" does not play a role in our current SDA workflow but I did find ways to do the sampling of soil ensemble members through the codes in "ensemble.R" by providing "soilinitcond" tag under "ensemble" in the XML file: https://github.com/PecanProject/pecan/blob/develop/modules/uncertainty/R/ensemble.R#L401-L406. However, it is only for a new fresh run and seems not applicable when it is a restart. The current restart components only include met inputs and parameters as https://github.com/PecanProject/pecan/blob/develop/modules/uncertainty/R/ensemble.R#L422C5-L424. What is your suggestion on that?
There was a problem hiding this comment.
- write.ensemble.configs is what's called by run.write.configs, so if that's what SDA is calling too we should be fine. But as noted earlier, it's important for the code to confirm that it's only receiving on ensemble member for each input, and to throw an error if it's not rather than resampling.
- Restart is grabbing
inputsand soils should be part of inputs (as should phenology, initial conditions, etc), not just met. You should verify this is working correctly, as it is is critical that ensemble sampling is preserved (not repeated) both iteratively within a SDA run and across SDA runs (e.g., when running a forecast today that starts from yesterday's forecast). Resampling things (params, inputs, IC, etc) that have already been samples is going to seriously mess up the covariance structures in our products
There was a problem hiding this comment.
@mdietze Based on the current argument passing of restart in "sda.enkf_MultiSite.R": https://github.com/PecanProject/pecan/blob/develop/modules/assim.sequential/R/sda.enkf_MultiSite.R#L358-L368, it seems we didn't grab soil, phenology and IC other than met now for inputs (note only "met" is specified in "input.ens.gen"). There is also a comment: #TODO: incorporate Phyllis's restart work. #sample all inputs specified in the settings$ensemble not just met". So I guess it is what we should add here?
There was a problem hiding this comment.
I added soil parameter files to inputs for the restart list. Phenology parameters are sampled within parameter sampling and initial conditions don't apply for t>1 situation, so I guess the current solution seems fine.
There was a problem hiding this comment.
Phenology parameters are sampled within parameter sampling
Not sure what you mean here, phenology files are separate from the parameter posterior files.
There was a problem hiding this comment.
So the general script is well established followed by a similar style as the previous SOC extraction function. The questions that I raised here are: 1) for the parallel design when over the site locations and data products, and 2) for the usage of downloaded local raster files (if they exist).
| name_tag <- expand.grid(depths, data_tag, stringsAsFactors = F) #find the combinations between data and depth tags. | ||
| L <- split(as.data.frame(name_tag), seq(nrow(as.data.frame(name_tag))))#convert tags into lists. | ||
|
|
||
| if ("try-error" %in% class(try(soilt_real <- L %>% furrr::future_map(function(l){ |
There was a problem hiding this comment.
When the site number goes up, this parallel will be more likely to produce errors due to the limitation of the SoilGrids server. The extraction speed is similar with or without the parallel in my previous 6400 site runs. Thus, it's recommended to just delete the parallel extraction part from the script.
| #LitterWHC | ||
| #param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*100 | ||
| #working on reading soil file | ||
| if (length(settings$run$inputs$soilinitcond$path) > 0) { |
There was a problem hiding this comment.
Still no clear how you're handling the case if >1 file is passed in
There was a problem hiding this comment.
So the current design is when t=1, only one file path is sampled and passed into the setting based on codes in "write.ensemble.configs". When t>1, file paths of total ensemble size will be passed in but then "template.soilinit" will be assigned with the specific path in the inputs.
There was a problem hiding this comment.
@Qianyuxuan I was asking about n=1 versus n>1, not t=1 vs t>1, and I'm still confused -- if you're passed in a whole vector of paths instead of a single one, you should throw an error not continue on, but I'm not seeing that check anywhere.
There was a problem hiding this comment.
I followed the steps for the pass of met files: https://github.com/PecanProject/pecan/blob/develop/models/sipnet/R/write.configs.SIPNET.R#L21-L29. I did find that a whole vector of met paths (10 ERA5 ensemble met paths defined in the settings) was passed in when t>1. But I don't think it is an error as the specific path will be overwritten by the one listed in the restart met input later.
There was a problem hiding this comment.
per previous conversation, I recommend adding a check to make sure that length(settings$run$inputs$soilinitcond$path) == 1 and throwing an error otherwise
| soil_IC_list <- PEcAn.data.land::pool_ic_netcdf2list(template.soilinit) | ||
| # Calculate the thickness of soil layers based on the assumption that the depth values are at bottoms and the first layer top is at 0 | ||
| if ("depth" %in% names(soil_IC_list$dims)) { | ||
| thickness<-c(soil_IC_list$dims$depth[1],diff(soil_IC_list$dims$depth)) |
There was a problem hiding this comment.
Does this line still work if the product only has one layer?
There was a problem hiding this comment.
Also, should/could there be an "else" here? A bunch of code below continues to check for "depth" being defined, but doesn't ever use depth again -- instead you're just using this to ensure thickness is defined.
There was a problem hiding this comment.
@Qianyuxuan sorry for not clarifying, the request isn't to have the code fail if thickness is missing, but to to only have depth/thickness within the if, to take all the other variable out of the if/else, and to use a default or user-specified soil depth if one isn't present in the file.
| soilWHC_total <- sum(unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])*thickness) | ||
| param[which(param[, 1] == "soilWHC"), 2] <- soilWHC_total | ||
| #LitterWHC | ||
| param[which(param[, 1] == "litterWHC"), 2] <- soilWHC_total |
There was a problem hiding this comment.
This is incorrect. If there's a litter layer it would have a much much smaller WHC than the entire soil pool.
| #litwaterDrainrate | ||
| param[which(param[, 1] == "litWaterDrainRate"), 2] <- unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1]*100/(3600*24) | ||
| param[which(param[, 1] == "litWaterDrainRate"), 2] <- soilHC_total |
There was a problem hiding this comment.
First, I'm not sure why you're using the average hydraulic conductivity of the entire soil column to control the rate at which moisture leaves the litter and enters the first soil layer. Second, we should check that the way SIPNET uses the conductivity (as a proportional decay rate, see https://github.com/PecanProject/sipnet/blob/becec4d2d6d857fa25dc47f974c48621a0b9044b/sipnet.c#L1516) is compatible with the way the parameter is defined in the soils file (Darcy's law) or whether additional assumptions or unit conversions are needed
There was a problem hiding this comment.
Fixed. Only the conductivity for first soil layer is used for litWaterDrainRate and the equation complies to Darcy's law.
| #LitterWHC | ||
| param[which(param[, 1] == "litterWHC"), 2] <- soilWHC_total | ||
| } | ||
| if ("depth" %in% names(soil_IC_list$dims) && "soil_hydraulic_conductivity_at_saturation" %in% names(soil_IC_list$vals)) { |
There was a problem hiding this comment.
per comment below about this being for surface litter, I'm not sure you need "depth" here
| #working on reading soil file (only working for soilWHC) | ||
| if (length(settings$run$inputs$soilinitcond$path) > 0) { | ||
| soil_ICs_num <- length(settings$run$inputs$soilinitcond$path) | ||
| soil_ICs_path <- settings$run$inputs$soilinitcond$path[[sample(1:soil_ICs_num, 1)]] |
There was a problem hiding this comment.
Phenology parameters are sampled within parameter sampling
Not sure what you mean here, phenology files are separate from the parameter posterior files.
| @@ -32,17 +32,18 @@ metSplit <- function(conf.settings, inputs, settings, model, no_split = FALSE, o | |||
| # Loading the model package - this is required bc of the furrr | |||
| library(paste0("PEcAn.",model), character.> | |||
|
|
|||
| inputs.split <- list() | |||
| # Initialize inputs.split with inputs to keep the sub-list of soil parameters "soilinitcond" | |||
There was a problem hiding this comment.
I'm not following the logic here. The metSplit code should only ever be touching the met. The comment here about soil parameters is confusing, as metSplit shouldn't be touching ANY of the other inputs, not just soils. I'm having a hard time telling if the code below is incorrect, or if it's fixing a long standing bug, or if the larger issue is that metSplit shouldn't ever have been passed the entirety of inputs to begin with.
There was a problem hiding this comment.
I didn't edit "metSplit" this time, instead I added similar modified codes in "sda.enkf_MultiSite.R". In my own opinion, the function "metSplit" lacks some flexibility as the formats "inputs = inputs$samples[[I]]", but we need to add another sublist such as "soilintcond" under inputs soil we need to specify "met".
| @@ -355,17 +355,28 @@ sda.enkf.multisite <- function(settings, | |||
| #reformatting params | |||
| new.params <- sda_matchparam(settings, ensemble.samples, site.ids, nens) | |||
| } | |||
| #sample met ensemble members | |||
| #sample met and soil parameter ensemble members | |||
There was a problem hiding this comment.
should be sampling ALL input ensemble members, not just met, not just (met + soils)
There was a problem hiding this comment.
#TODO: incorporate Phyllis's restart work # sample all inputs specified in the settings$ensemble not just met
There was a problem hiding this comment.
FYI Phyllis's work lives on the BU server at /projectnb/dietzelab/yixuanli and still needs going through to make sure (a) write.ensemble.configs is saving ALL the info it needs about ensembles, not just params and (b) that you can have the option to pass in a specific already-sampled ensemble of params, met, soil, IC, etc. without having the code resample it. This is important in the general case to allow for sensitivity analysis designs. It's also important for the SDA case as the sample from t==1 needs to be the same sample that's used for t>1.
There was a problem hiding this comment.
I did find the previous design had two separate sampling in the SDA workflow, leading to met file for t==1 differs from that for t>1. I kept the sampling in SDA codes instead, saved the inputs as "samples", and passed it into "write.ensemble.configs" as an additional argument. If it is a fresh new run, those samples will be read and recorded directly. And those samples will be passed into the restart list when t>1 as well. Although it might not be the perfect solution, but I did test it and found it work well to make sure the ensemble member keeps the same through time.
| parent_ids = NULL | ||
| ) | ||
| }) | ||
| input_name <- c("met","soilinitcond") |
There was a problem hiding this comment.
should not be hard coded. The set of inputs that are ensemble-based needs to be detected from the settings
| input_name <- c("met","soilinitcond") | ||
| inputs <- list() | ||
| new_inputs <- list() | ||
| for (i_input in seq_along(input_name)){ |
There was a problem hiding this comment.
Inputs can't be sampled in an arbitrary order. See the logic in write.ensemble.configs. Indeed, if the SDA can't call write.ensemble.configs directly, then there's a lot of code there than needs to be refactored into a function that can be called by both forward and SDA code, since maintaining two separate versions of the sampling code is clearly resulting in code divergence.
There was a problem hiding this comment.
In my new edits, the ids for samples that were generated in SDA codes and passed into the "write.ensemble.configs" will be recorded in "README.txt".
| @@ -225,6 +225,7 @@ template <- PEcAn.settings::Settings(list( | |||
| ensemble = structure(list(size = 25, variable = "NPP", | |||
| samplingspace = structure(list( | |||
| parameters = structure(list(method = "lhc")), | |||
| soilinitcond = structure(list(method = "sampling")), | |||
There was a problem hiding this comment.
Check to make sure there aren't other inputs that are being left out
There was a problem hiding this comment.
Also, is the operational SDA using a LHC sampling of the posterior parameters? That's not random.
There was a problem hiding this comment.
I think so. The current SDA workflow used "get.parameter.samples(settings, ens.sample.method = settings$ensemble$samplingspace$parameters$method)" to sample parameters and the method is "lhc".
There was a problem hiding this comment.
Code still retains lhc instead of sampling
| #LitterWHC | ||
| #param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*100 | ||
| #working on reading soil file | ||
| if (length(settings$run$inputs$soilinitcond$path) > 0) { |
There was a problem hiding this comment.
@Qianyuxuan I was asking about n=1 versus n>1, not t=1 vs t>1, and I'm still confused -- if you're passed in a whole vector of paths instead of a single one, you should throw an error not continue on, but I'm not seeing that check anywhere.
| if (!is.null(inputs)) { | ||
| ## override if specified in inputs | ||
| if ("soilinitcond" %in% names(inputs)) { | ||
| template.soilinit <- inputs$soilinit$path |
There was a problem hiding this comment.
I'm now realizing that I'm also confused by this bit. Why do settings$run$inputs$soilinitcond and settings$run$inputs$soilinit both exist?
There was a problem hiding this comment.
"soilinit" is only a typo and will be corrected.
…an into sda_soil_params_soilgrids
| ##' @return The individual alphas that work best to fit the observed quantiles | ||
| ##' @export | ||
| ##' @author Qianyu Li | ||
| estimate_dirichlet_parameters <- function(means, quantiles) { |
There was a problem hiding this comment.
I don't think it should be combined with the existing function prior_fn (that is similar in approach but for univariate distributions though it has its own limitations). But perhaps it should be included in the same file in the priors package?
pecan/modules/priors/R/priors.R
Line 90 in edc14b8
Co-authored-by: DongchenZ <82834956+DongchenZ@users.noreply.github.com>
| #LitterWHC | ||
| #param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*100 | ||
| #working on reading soil file | ||
| if (length(settings$run$inputs$soilinitcond$path) > 0) { |
There was a problem hiding this comment.
per previous conversation, I recommend adding a check to make sure that length(settings$run$inputs$soilinitcond$path) == 1 and throwing an error otherwise
| @@ -355,17 +355,33 @@ sda.enkf.multisite <- function(settings, | |||
| #reformatting params | |||
| new.params <- sda_matchparam(settings, ensemble.samples, site.ids, nens) | |||
| } | |||
| #sample met ensemble members | |||
| #TODO: incorporate Phyllis's restart work | |||
There was a problem hiding this comment.
please retain this comment, which will be helpful to the folks working on uncertainty partitioning.
| } | ||
| for (i in seq_along(samp.ordered)){ | ||
| #call the function responsible for generating the ensemble | ||
| inputs[[s]][[names(samp.ordered)[i]]] <- input.ens.gen(settings=conf.settings[[s]], |
There was a problem hiding this comment.
is this sampling occurring for each site independently or across all sites simultanously?
There was a problem hiding this comment.
For each site independently.
There was a problem hiding this comment.
Is independent sampling (which is NOT what we want) something new here or are you just continuing the status quo. If it's the status quo I can see punting the fix to a new PR, but if you moved from joint to independent sampling you'll need to revert it back to joint.
There was a problem hiding this comment.
I didn't do anything new here and just continued the status quo. The ensemble size for met is 10 at each site, does it mean we may need to sample from the joint distribution using bootstrapping?
There was a problem hiding this comment.
No, what I mean is that if model ensemble member 1 at site 1 samples met ensemble k (e.g. 6) then model ensemble member 1 needs to also use met ensemble member k at sites 2 through 6400. Even though the specific met data is different from site to site, the spatiotemporal values within a specific met ensemble member share a covariance structure. Similarly if ensemble member 1 at site 1 samples a particular parameter vector with a specific row number, then all other sites need to use the same posterior row number (though row k for one PFT will obviously point to a different set of parameters than row k of another PFT). By having the same parameters, met, etc across sites we induce the correct spatiotemporal covariance structure across sites.
| @@ -379,7 +395,42 @@ sda.enkf.multisite <- function(settings, | |||
| if (t>1){ | |||
| #for next time step split the met if model requires | |||
| #-Splitting the input for the models that they don't care about the start and end time of simulations and they run as long as their met file. | |||
| inputs.split <- metSplit(conf.settings, inputs, settings, model, no_split = FALSE, obs.times, t, nens, restart_flag = FALSE, my.split_inputs) | |||
| #set start and end date for splitting met | |||
There was a problem hiding this comment.
Conceptually, I'm not seeing why you're replacing a simple call to a function with a big chunk of code, which undo's the concept of modularization and refactoring. If there's something wrong within the function, I'd fix the function.
Also, conceptually note that MOST models don't need their met split, so make sure that any implementation of metSplit
There was a problem hiding this comment.
Previously, I directly modified "metSplit" but it did not seem to be suggested. The challenge is it that we only need to do the split function for "met" sublist but now it lacks the specification: https://github.com/PecanProject/pecan/blob/develop/modules/assim.sequential/R/metSplit.R#L39. Do you suggest we change it to:" inputs.split$met$samples[i] <- do.call(
my.split_inputs,
args = list(
settings = settings,
start.time = (lubridate::ymd_hms(start.time, truncated = 3) + lubridate::second(lubridate::hms("00:00:01"))),
stop.time = lubridate::ymd_hms(stop.time, truncated = 3),
inputs = inputs$met$samples[[i]])
| inputs.split <- | ||
| furrr::future_pmap(list(conf.settings %>% `class<-`(c("list")), inputs, model), function(settings, inputs, model) { | ||
| # Loading the model package - this is required bc of the furrr | ||
| library(paste0("PEcAn.",model), character.> |
There was a problem hiding this comment.
We shouldn't have to load library's like this. But if you do, load it once, before the pmap, not repeatedly.
There was a problem hiding this comment.
Has removed this block of code.
| @@ -225,6 +225,7 @@ template <- PEcAn.settings::Settings(list( | |||
| ensemble = structure(list(size = 25, variable = "NPP", | |||
| samplingspace = structure(list( | |||
| parameters = structure(list(method = "lhc")), | |||
| soilinitcond = structure(list(method = "sampling")), | |||
There was a problem hiding this comment.
Code still retains lhc instead of sampling
modules/uncertainty/R/ensemble.R
Outdated
| @@ -393,6 +394,7 @@ write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, | |||
| "site : ", settings$run$site$name, "\n", | |||
| "site id : ", format(settings$run$site$id, scientific = FALSE), "\n", | |||
| "met data : ", samples$met$samples[[i]], "\n", | |||
| "soil parameter : ", samples$soilinitcond$samples[[i]], "\n", | |||
There was a problem hiding this comment.
This line will fail for any model that doesn't have soilinitcond. You should dynamically detect the set of sampled inputs present
modules/uncertainty/R/ensemble.R
Outdated
| RENAME = rename)#for restart from previous model runs, not sharing the same outdir | ||
| ) | ||
| } | ||
| params<-new.params | ||
| return(invisible(list(runs = data.frame(id=run.id), ensemble.id = ensemble.id, samples=list(met=inputs) | ||
| return(invisible(list(runs = data.frame(id=run.id), ensemble.id = ensemble.id, samples=list(met=inputs$met,soilinitcond=inputs$soilinitcond) |
There was a problem hiding this comment.
same need to generalize to any set of inputs
…an into sda_soil_params_soilgrids
| soil_IC_list <- PEcAn.data.land::pool_ic_netcdf2list(template.soilinit) | ||
| # Calculate the thickness of soil layers based on the assumption that the depth values are at bottoms and the first layer top is at 0 | ||
| if ("depth" %in% names(soil_IC_list$dims)) { | ||
| thickness<-c(soil_IC_list$dims$depth[1],diff(soil_IC_list$dims$depth)) |
There was a problem hiding this comment.
@Qianyuxuan sorry for not clarifying, the request isn't to have the code fail if thickness is missing, but to to only have depth/thickness within the if, to take all the other variable out of the if/else, and to use a default or user-specified soil depth if one isn't present in the file.
| } | ||
| for (i in seq_along(samp.ordered)){ | ||
| #call the function responsible for generating the ensemble | ||
| inputs[[s]][[names(samp.ordered)[i]]] <- input.ens.gen(settings=conf.settings[[s]], |
There was a problem hiding this comment.
Is independent sampling (which is NOT what we want) something new here or are you just continuing the status quo. If it's the status quo I can see punting the fix to a new PR, but if you moved from joint to independent sampling you'll need to revert it back to joint.
| } | ||
| input_path <- conf.settings[[s]]$run$inputs[[tolower(input_tag)]]$path | ||
| inputs[[s]][[input_tag]]$ids<-inputs[[random_site]][[input_tag]]$ids | ||
| inputs[[s]][[input_tag]]$samples<- input_path[inputs[[random_site]][[input_tag]]$ids] |
There was a problem hiding this comment.
This is a good start. Tagging @blesson07asd so he's aware of these changes and can account for them when he further generalizes ensemble sampling. Eventually this will need to be expanded to account for parameter sampling too, but we can do that later.
|
@Qianyuxuan there's a small fix that needs to occur in data.land/DESCRIPTION to add the following new libraries to Suggests: ‘doSNOW’ ‘foreach’ ‘MCMCpack’. After that you'll probably need to run scripts/generate_dependencies.R and add any updated files to the PR (e.g., docker/depends/pecan_package_dependencies.csv) |
Add codes to extract SoilGrids soil texture data and derive ensemble soil parameter files
Description
Motivation and Context
Provide ensemble soil parameter files based on SoilGrids soil texture data for SDA.
Review Time Estimate
Types of changes
Checklist: