Pecan disturbance da - #4123
Pecan disturbance da#4123DongchenZ wants to merge 3 commits into
Conversation
mdietze
left a comment
mdietze
left a comment
There was a problem hiding this comment.
modules/assim.sequential/inst/anchor/disturbanceDA_scripts would really benifit from a README explaining what the other scripts are and the order they need to be run in. For example, 5 different scripts (everything but downscale) seem to be associated with data prep. It's ambiguous which script actually runs the SDA. I'd also suggest a git mv to move the entire directory out of /inst/anchor and instead just have in /inst as it's a completely distinct analysis from the continental SDA.
| library(lgarch) | ||
| library(parallel) | ||
| #read settings. | ||
| settings <- PEcAn.settings::read.settings("/projectnb/dietzelab/dongchen/SDA_disturbance/example_polygon/pecan.xml") |
There was a problem hiding this comment.
please add this pecan.xml to the PR
| ```{r} | ||
| library(randomForest) | ||
| basedir <- "/projectnb/dietzelab/dongchen/SDA_disturbance/example_polygon/downscale_dist_DA/" | ||
| rf.train.downscale <- function(dat.train, predictor) { |
There was a problem hiding this comment.
probably a future direction, but this seems like a place where we really need to converge on a set of tools for doing the downscaling, rather than reinventing (or copy-pasting) long functions from script-to-script. In the coming months I'd really like to see some harmonization across what you're doing, and what Julian, Yang, and David are doing.
| CA <- array(NA, dim = c(n.comp, n.state, n.state)) | ||
| wA <- array(NA, dim = c(n.comp)) | ||
|
|
||
| # select undisturbed / disturbed forecast |
There was a problem hiding this comment.
Future direction: you don't need to assign each ensemble member to disturbed or undisturbed, that just throws away 1/2 of your SIPNET runs. Instead just make a second ensemble that applies disturbance to a copy of every member of the first ensemble
|
|
||
| # kalman update to compute analysis means, covariances, weights | ||
| n.obs <- length(Y) | ||
| pD <- c(0.9, 0.1) # prior probability of disturbance 0.1, default value which can be tweaked |
There was a problem hiding this comment.
These are really bad defaults. Need to be able to pass in data-driven priors by pixel and year.
| for (i in 1:n.obs) { | ||
| H.mat[i,H[i]] <- 1 | ||
| } | ||
| for (d in 1:2) { |
There was a problem hiding this comment.
Improve documentation. This looks like you're just applying a standard Kalman update to each of the two disturbance classes (disturbed and undisturbed). Also, it's not clear exactly what lik and wA are (they look like the likelihood and some sort of weight, but it's not documented)
| if (as.logical(block$disturbance)) { | ||
| load(block$dist.prior.file) # load disturbance prior file. | ||
| # prescribe mu0 for the disturbance-related reductions in AGB and LAI. | ||
| mu0 <- c(att$agb.bio.reduce[as.numeric(block$site.ids)], att$lai.bio.reduce[as.numeric(block$site.ids)]) |
There was a problem hiding this comment.
I can't find where in the code att is being defined
| ##' @importFrom dplyr %>% | ||
| ##' @export | ||
|
|
||
| disturbance_DA <- function (X, Y, R, H, aqq, bqq, disturbance, mu0=rep(80, 80), V0=diag(rep(0.25, 2))) { |
There was a problem hiding this comment.
Default for mu0 don't make sense as these parameters have to be between 0 and 1.
| new.biomass <- tobit(mvtnorm::rmvnorm(1,mu0*old.biomass,tcrossprod(old.biomass,old.biomass)*V0)) ## draw disturbed leaf and stem | ||
| check.max <- which(new.biomass > old.biomass) | ||
| if(length(check.max)>0) new.biomass[check.max] = old.biomass[check.max] # constrain new.biomass[i] <= old.biomass[i] | ||
| residual = sum(old.biomass-new.biomass) |
There was a problem hiding this comment.
not an intuitive name for biomass lost to disturbance (e.g., for fire this is the combustion). Also, you need to make sure you're explicitly tracking removals as they are important model outputs.
| if(length(check.max)>0) new.biomass[check.max] = old.biomass[check.max] # constrain new.biomass[i] <= old.biomass[i] | ||
| residual = sum(old.biomass-new.biomass) | ||
| x[bio.idxs] <- new.biomass # update state with disturbed biomass | ||
| x[4] <- x[4] + residual*alloc.soil # some of the biomass gets transported to the soil |
There was a problem hiding this comment.
Interesting that you're passing in the vegetative pools via a variable but you're hard coding that x[4] has to be soils. This seems a bit dangerous. More broadly, it's not safe to assume that the same fraction of leaf and stem go to the soil. Better solution and more general solution(not this PR) would be a full transition matrix.
| residual = sum(old.biomass-new.biomass) | ||
| x[bio.idxs] <- new.biomass # update state with disturbed biomass | ||
| x[4] <- x[4] + residual*alloc.soil # some of the biomass gets transported to the soil | ||
| removal <- residual*(1-alloc.soil) |
There was a problem hiding this comment.
similar to earlier comment, losses are being calculated but then never saved/reported
2 participants
Add the disturbance DA algorithm to the existing SDA analysis pipeline
Description
In this PR, I added three things:
/pecan/modules/assim.sequential/inst/anchor/disturbanceDA_scripts.disturbance_DAfunction for incorporating disturbance dynamics into the DA workflow --/pecan/modules/assim.sequential/R/disturbance_DA.R.disturbance_DAfunction inside the SDA analysis function/pecan/modules/assim.sequential/R/Analysis_sda_block.R.Motivation and Context
Review Time Estimate
Types of changes
Checklist: