Add functions that can extract and update phenological status based on MODIS data and enable monthly SDA

base/workflow/R/do_conversions.R

@@ -73,6 +73,13 @@ do_conversions <- function(settings, overwrite.met = FALSE, overwrite.fia = FALS
       ## which is done locally in rundir and then rsync'ed to remote
       ## rather than having a model-format soils file that is processed remotely
     }
+
+    # 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)
+      needsave <- TRUE
+    }
+
     # met conversion
 
     if (input.tag == "met") {

book_source/03_topical_pages/04_R_workflow.Rmd

@@ -121,6 +121,12 @@ The main script that handles Met Processing, is [`met.process`](https://github.c
 
 ### Converting from PEcAn standard to model-specific format {#workflow-met-model}
 
+## Input phenological Data {#workflow-input-phenology}
+
+To enable the use of MODIS phenology data (MODIS Land Cover Dynamics (MCD12Q2)) to update the phenological parameters (leaf-on date) and (leaf-off date) at each restart timepoint: 
+ 1. Generate the phenological parameter CSV file by running `PEcAn.data.remote::extract_phenology_MODIS`.
+ 2. Provide the generated phenological parameter CSV file to `settings$run$inputs$leaf_phenology$path`.
+
 ## Traits {#workflow-traits}
 
 (TODO: Under construction)
@@ -130,6 +136,7 @@ The main script that handles Met Processing, is [`met.process`](https://github.c
 (TODO: Under construction)
 
 ## Model Configuration {#workflow-modelconfig}
+To enable the state data assimilation with sub-annual data, default `conflict` in `model2netcdf` should be `TRUE`. 
 
 (TODO: Under construction)
 

models/sipnet/R/model2netcdf.SIPNET.R

@@ -140,7 +140,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
   for (y in year_seq) {
     #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) {
       next
@@ -296,13 +296,19 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       close(varfile)
       ncdf4::nc_close(nc)
 
-      #merge nc files
+      #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)
+      ncdf4::nc_close(nc)
       unlink(files, recursive = T)
     }else{
       nc      <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var)
@@ -323,4 +329,4 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
   }
 } # model2netcdf.SIPNET
 #--------------------------------------------------------------------------------------------------#
-### EOF
+### EOF

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

@@ -11,6 +11,14 @@
 ##' Writes a configuration files for your model
 ##' @name write.config.SIPNET
 ##' @title Writes a configuration files for SIPNET model
+##' @param defaults pft
+##' @param trait.values vector of samples for a given trait
+##' @param settings PEcAn settings object
+##' @param run.id run ID
+##' @param inputs list of model inputs
+##' @param IC initial condition
+##' @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 spinup currently unused, included for compatibility with other models
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
@@ -438,7 +446,31 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     if ("leafGrowth" %in% pft.names) {
       param[which(param[, 1] == "leafGrowth"), 2] <- pft.traits[which(pft.names == "leafGrowth")]
     }
-  }  ## end loop over PFTS
+
+    #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 <-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafonday)
+       leafOffDay<-  dplyr::filter(leafphdata, year== obs_year_start & site_id==settings$run$site$id) %>%  dplyr::select (leafoffday)
+       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.")
+      }
+    }
+  } ## 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){
@@ -467,7 +499,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       }else{
         wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
       }
-      
+
       #Sanity check
       if (is.infinite(wood_total_C) | is.nan(wood_total_C) | wood_total_C < 0) {
         wood_total_C <- 0
@@ -536,6 +568,15 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (!is.na(lai) && is.numeric(lai)) {
         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       
+          }
+      }
       ## neeInit gC/m2
       nee <- try(ncdf4::ncvar_get(IC.nc,"nee"),silent = TRUE)
       if (!is.na(nee) && is.numeric(nee)) {

models/sipnet/R/write_restart.SIPNET.R

@@ -28,7 +28,7 @@
 ##' @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
@@ -59,15 +59,15 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
     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")
   }
-  
+
   if ("NEE" %in% variables) {
     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")
-    analysis.save[[length(analysis.save) + 1]] <- AbvGrndWood 	  
+    analysis.save[[length(analysis.save) + 1]] <- AbvGrndWood
     names(analysis.save[[length(analysis.save)]]) <- c("AbvGrndWood")
   }
 
@@ -106,7 +106,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
     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  
     if (new.state$LAI < 0) analysis.save[[length(analysis.save)]] <- 0
@@ -120,7 +120,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
   }else{
     analysis.save.mat <- NULL
   }
-  
+
   if (verbose) {
     print(runid %>% as.character())
     print(analysis.save.mat)

modules/assim.sequential/R/Analysis_sda_block.R

@@ -112,9 +112,9 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
     PEcAn.logger::logger.warn("The zero variances in Pf is being replaced by one fifth of the minimum variance in those matrices respectively.")
   }
   #distance calculations and localization
-  site.locs <- settings$run %>% 
-    purrr::map('site') %>% 
-    purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>% 
+  site.locs <- settings$run %>%
+    purrr::map('site') %>%
+    purrr::map_dfr(~c(.x[['lon']],.x[['lat']]) %>% as.numeric)%>%
     t %>%
     `colnames<-`(c("Lon","Lat")) %>%
     `rownames<-`(site.ids)
@@ -172,9 +172,9 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
     # if there is any site that has zero observation.
     if (any(obs_per_site == 0)) {
       #name matching between observation names and state variable names.
-      f.2.y.ind <- obs.mean[[t]] %>% 
-        purrr::map(\(x)which(var.names %in% names(x))) %>% 
-        unlist %>% 
+      f.2.y.ind <- obs.mean[[t]] %>%
+        purrr::map(\(x)which(var.names %in% names(x))) %>%
+        unlist %>%
         unique
       H <- list(ind = f.2.y.ind %>% purrr::map(function(start){
         seq(start, length(site.ids) * length(var.names), length(var.names))
@@ -283,9 +283,9 @@ build.block.xy <- function(settings, block.list.all, X, obs.mean, obs.cov, t) {
         if (is.null(obs.mean[[t]])) {
           f.2.y.ind <- seq_along(var.names)
         } else {
-          f.2.y.ind <- obs.mean[[t]] %>% 
-            purrr::map(\(x)which(var.names %in% names(x))) %>% 
-            unlist %>% 
+          f.2.y.ind <- obs.mean[[t]] %>%
+            purrr::map(\(x)which(var.names %in% names(x))) %>%
+            unlist %>%
             unique
         }
         seq.ind <- f.2.y.ind %>% purrr::map(function(start){
@@ -430,7 +430,7 @@ MCMC_block_function <- function(block) {
   conf$addSampler(target = samplerLists[[X.mod.ind]]$target, type = "ess",
                   control = list(propCov= block$data$pf, adaptScaleOnly = TRUE,
                                  latents = "X", pfOptimizeNparticles = TRUE))
-  
+
   #add toggle Y sampler.
   for (i in 1:block$constant$YN) {
     conf$addSampler(paste0("y.censored[", i, "]"), 'toggle', control=list(type='RW'))
@@ -451,7 +451,6 @@ MCMC_block_function <- function(block) {
 
   #run MCMC
   dat <- runMCMC(Cmcmc, niter = block$MCMC$niter, nburnin = block$MCMC$nburnin, thin = block$MCMC$nthin, nchains = block$MCMC$nchain)
-
   #update aq, bq, mua, and pa
   M <- colMeans(dat)
   block$update$aq <- block$Inits$q

modules/assim.sequential/R/Nimble_codes.R

@@ -187,14 +187,14 @@ GEF.MultiSite.Nimble <-  nimbleCode({
       # Sorting out qs
       q[1:YN, 1:YN] ~ dwish(R = aq[1:YN, 1:YN], df = bq) ## aq and bq are estimated over time
     }
-    
+
     for (i in 1:nH) {
       tmpX[i]  <- X.mod[H[i]]
       Xs[i] <- tmpX[i]
     }
     ## add process error to x model but just for the state variables that we have data and H knows who
     X[1:YN]  ~ dmnorm(Xs[1:YN], prec = q[1:YN, 1:YN])
-    
+
     ## Likelihood
     y.censored[1:YN] ~ dmnorm(X[1:YN], prec = r[1:YN, 1:YN])
 

modules/assim.sequential/R/SDA_OBS_Assembler.R

@@ -36,7 +36,7 @@ SDA_OBS_Assembler <- function(settings){
   }
 
   #prepare site_info offline, because we need to submit this to server remotely, which might not support the Bety connection.
-  site_info <- settings$run %>% 
+  site_info <- settings$run %>%
     purrr::map('site')%>% 
     purrr::map(function(site.list){
       #conversion from string to number
@@ -189,12 +189,13 @@ SDA_OBS_Assembler <- function(settings){
     for (j in seq_along(obs.mean[[i]])) {
       if (sum(is.na(obs.mean[[i]][[j]]))){
         na_ind <- which(is.na(obs.mean[[i]][[j]]))
-        obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
-        if(length(new_diag(obs.cov[[i]][[j]])) == 1){
+        #obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
+        if(length(obs.mean[[i]][[j]]) == 1){
           obs.cov[[i]][[j]] <- obs.cov[[i]][[j]][-na_ind]
         }else{
           obs.cov[[i]][[j]] <- obs.cov[[i]][[j]][-na_ind, -na_ind]
         }
+	obs.mean[[i]][[j]] <- obs.mean[[i]][[j]][-na_ind]
       }
       SoilC_ind <- which(names(obs.mean[[i]][[j]]) == "TotSoilCarb")
       if (length(SoilC_ind) > 0){
@@ -217,7 +218,7 @@ SDA_OBS_Assembler <- function(settings){
                          Obs_Prep[var_ind] %>% purrr::map(~.x$end.date)),
                     function(var_timestep, var_start_date, var_end_date){
                       obs_timestep2timepoint(var_start_date, var_end_date, var_timestep)
-                    }) %>% 
+                    }) %>%
         purrr::map(function(all_timepoints){
           all_timepoints[which(!all_timepoints %in% time_points)]
         }) %>% 

modules/assim.sequential/inst/MultiSite-Exs/SDA/Create_Multi_settings.R

@@ -17,6 +17,7 @@ ERA5_dir <- "/projectnb/dietzelab/dongchen/anchorSites/ERA5_2012_2021/"
 XML_out_dir <- "/projectnb/dietzelab/dongchen/anchorSites/SDA/pecan.xml"
 
 pft_csv_dir <- "/projectnb/dietzelab/dongchen/anchorSites/site_pft.csv"
+modis_phenology_dir <- "/projectnb/dietzelab/Cherry/pft_files/leaf_phenology.csv"
 
 #Obs_prep part
 #AGB
@@ -290,7 +291,8 @@ template <- PEcAn.settings::Settings(list(
       #                               )),
       # soilinitcond = structure(list(path = "/projectnb/dietzelab/ahelgeso/EFI_Forecast_Challenge/"
       #                               )),
-      pft.site = structure(list(path = pft_csv_dir))
+      pft.site = structure(list(path = pft_csv_dir)),
+      leaf_phenology = structure(list(path = modis_phenology_dir))
     ))
   ))
 ))