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

models/sipnet/R/model2netcdf.SIPNET.R

@@ -97,13 +97,14 @@ sipnet2datetime <- function(sipnet_tval, base_year, base_month = 1,
 ##' @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.
+##' conflict is set to TRUE to enable the assimilation of monthly data.
 ##' @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) {
 
   ### Read in model output in SIPNET format
   sipnet_out_file <- file.path(outdir, prefix)
@@ -296,13 +297,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 subannual 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 +330,4 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
   }
 } # model2netcdf.SIPNET
 #--------------------------------------------------------------------------------------------------#
-### EOF
+### EOF

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

@@ -14,7 +14,7 @@
 ##' @export
 ##' @author Michael Dietze
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
-                                restart = NULL, spinup = NULL) {
+                                restart = NULL, spinup = NULL,obs_time=NULL,update_phenology=NULL) {
   ### WRITE sipnet.in
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
   config.text <- readLines(con = template.in, n = -1)
@@ -437,7 +437,26 @@ 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 (update_phenology){
+     leaf_pheno_outdir <- settings$model$leaf_phenology$outdir  ## read from settings
+     if (!is.null(leaf_pheno_outdir)) {
+    ##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]
+       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 modules/data.remote/R/phenology_MODIS_extract.R 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){
@@ -466,7 +485,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
@@ -522,7 +541,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     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))
+    IC.pools <- PEcAn.data.land::prepare_pools(IC.path, settings, constants = list(sla = SLA))
 
     if(!is.null(IC.pools)){
       IC.nc <- ncdf4::nc_open(IC.path) #for additional variables specific to SIPNET
@@ -630,4 +649,4 @@ remove.config.SIPNET <- function(main.outdir, settings) {
   } else {
     print("*** WARNING: Removal of files on remote host not yet implemented ***")
   }
-} # remove.config.SIPNET 
+} # remove.config.SIPNET

models/sipnet/R/write_restart.SIPNET.R

@@ -20,15 +20,17 @@
 ##' @param RENAME flag to either rename output file or not
 ##' @param new.params list of parameters to convert between different states 
 ##' @param inputs list of model inputs to use in write.configs.SIPNET
+##' @param obs_time obervation timepoints
+##' @param update_phenology TRUE if we want to update the phenological data (i.e. leaf-on and leaf-off dates) for each restart run during SDA
 ##' @param verbose decide if we want to print the outputs.
 ##'
 ##' @description Write restart files for SIPNET. WARNING: Some variables produce illegal values < 0 and have been hardcoded to correct these values!!
 ##' 
 ##' @return NONE
 ##' @export
 write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings, new.state,
-                                 RENAME = TRUE, new.params = FALSE, inputs, verbose = FALSE) {
-  
+                                 RENAME = TRUE, new.params = FALSE, inputs, obs_time=NULL, update_phenology=NULL, 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 +61,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 +108,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 +122,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)
@@ -130,6 +132,8 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                                            settings = settings,
                                            run.id = runid,
                                            inputs = inputs,
-                                           IC = analysis.save.mat))
+                                           IC = analysis.save.mat,
+                                           obs_time=obs_time,
+					                                 update_phenology=update_phenology))
   print(runid)
 } # write_restart.SIPNET

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/R/sda.enkf_MultiSite.R

@@ -50,7 +50,8 @@ sda.enkf.multisite <- function(settings,
                                             keepNC = TRUE,
                                             forceRun = TRUE,
                                             run_parallel = TRUE,
-                                            MCMC.args = NULL),
+                                            MCMC.args = NULL,
+                                            update_phenology = TRUE),
                                ...) {
   #add if/else for when restart points to folder instead if T/F set restart as T
   if(is.list(restart)){
@@ -296,7 +297,8 @@ sda.enkf.multisite <- function(settings,
             settings = settings,
             model = settings$model$type,
             write.to.db = settings$database$bety$write,
-            restart = restart.arg
+            restart = restart.arg,
+            update_phenology=control$update_phenology
           )
         }) %>%
         stats::setNames(site.ids)
@@ -420,7 +422,9 @@ sda.enkf.multisite <- function(settings,
               model = settings$model$type,
               write.to.db = settings$database$bety$write,
               restart = restart.arg,
-              rename = TRUE
+              rename = TRUE,
+              time = obs.year,
+              update_phenology=control$update_phenology
             )
           }) %>%
           stats::setNames(site.ids)
@@ -530,7 +534,7 @@ sda.enkf.multisite <- function(settings,
           if (is.null(control$MCMC.args)) {
             MCMC.args <- list(niter = 1e5,
                               nthin = 10,
-                              nchain = 3,
+                              nchain = 1,
                               nburnin = 5e4)
           } else {
             MCMC.args <- control$MCMC.args

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

@@ -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",
+                         jobtemplate = "~/sipnet_geo.job",
+                         leaf_phenology=structure(list(outdir="USER_DEFINED"))
   )),
 
   ###########################################################################