Adding SDA restart capability

CHANGELOG.md

@@ -54,6 +54,7 @@ convert data for a single PFT fixed (#1329, #2974, #2981)
   (#2989; @nanu1605)
 - Fixed a bug with ED2 where ED2IN tags supplied in `settings` that were not in the ED2IN template file were not getting added to ED2IN config files (#3034, #3033)
 - Fixed a bug where warnings were printed for file paths on remote servers even when they did exist (#3020)
+- Fixed bug in model2netcdf.SIPNET that caused LE to be overestimaed 10^3 (#3036)
 - Added an updated ED2IN template file, `models/ed/inst/ED2IN.r2.2.0.github`, to reflect new variables in the development version of ED2
 - `PEcAn.data.land::gSSURGO.Query` has been updated to work again after changes to the gSSURGO API.
 - `PEcAn.settings::read.settings()` now prints a warning when falling back on default `"pecan.xml"` if the named `inputfile` doesn't exist.

base/settings/R/check.all.settings.R

@@ -67,12 +67,17 @@ check.inputs <- function(settings) {
         }
       } else if ("path" %in% names(settings$run$inputs[[tag]])) {
         # can we find the file so we can set the tag.id
-        id <- PEcAn.DB::dbfile.id(
-          "Input",
-          settings$run$inputs[[tag]][["path"]],
-          dbcon,
-          hostname)
-        if (!is.na(id)) {
+        #adding for to loop over ensemble member filepaths 
+        id <- list()
+        path <- settings$run$inputs[[tag]][["path"]]
+        for (j in 1:length(path)){
+          id[j] <- PEcAn.DB::dbfile.id(
+            "Input",
+            path[[j]],
+            dbcon,
+            hostname)
+        }
+        if (any(!is.na(id))) {
           settings$run$inputs[[tag]][["id"]] <- id
         }
       }

models/sipnet/R/met2model.SIPNET.R

@@ -147,7 +147,9 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date
       }
 
       Rain <- ncdf4::ncvar_get(nc, "precipitation_flux")
-      pres <- ncdf4::ncvar_get(nc,'air_pressure') ## in pascal
+
+      press <- ncdf4::ncvar_get(nc,'air_pressure') ## in pascal
+
       SW <- ncdf4::ncvar_get(nc, "surface_downwelling_shortwave_flux_in_air")  ## in W/m2
 
       PAR <- try(ncdf4::ncvar_get(nc, "surface_downwelling_photosynthetic_photon_flux_in_air"))  ## in umol/m2/s
@@ -172,12 +174,14 @@ met2model.SIPNET <- function(in.path, in.prefix, outfolder, start_date, end_date
       SVP <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::get.es(Tair_C), "millibar", "Pa")  ## Saturation vapor pressure
       VPD <- try(ncdf4::ncvar_get(nc, "water_vapor_saturation_deficit"))  ## in Pa
       if (!is.numeric(VPD)) {
-        VPD <- SVP * (1 - PEcAn.data.atmosphere::qair2rh(Qair, Tair_C, pres/100))
+
+        VPD <- SVP * (1 - PEcAn.data.atmosphere::qair2rh(Qair, Tair_C, press = press/100))
+
         PEcAn.logger::logger.info("water_vapor_saturation_deficit absent; VPD calculated from Qair, Tair, and SVP (saturation vapor pressure) ")
       }
       e_a <- SVP - VPD
       VPDsoil <- PEcAn.utils::ud_convert(PEcAn.data.atmosphere::get.es(soilT), "millibar", "Pa") *
-        (1 - PEcAn.data.atmosphere::qair2rh(Qair, soilT, pres/100))
+        (1 - PEcAn.data.atmosphere::qair2rh(Qair, soilT, press/100))
 
       ncdf4::nc_close(nc)
     } else {

models/sipnet/R/model2netcdf.SIPNET.R

@@ -196,15 +196,15 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     output[[10]] <- (sub.sipnet.output$plantWoodC * 0.001) + (sub.sipnet.output$plantLeafC * 0.001) + 
       (sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001 # Total living C kgC/m2
     output[[11]] <- (sub.sipnet.output$soil * 0.001) + (sub.sipnet.output$litter * 0.001)  # Total soil C kgC/m2
-    if (revision == "r136") {
-      output[[12]] <- (sub.sipnet.output$evapotranspiration * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s  # Qle W/m2
-    } else {
+    if (revision == "unk") {
       ## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : ***
       ## *** it says "npp" in the header (written by L774) but the values being written are trackers.evapotranspiration (L806) ***
       ## evapotranspiration in SIPNET is cm^3 water per cm^2 of area, to convert it to latent heat units W/m2 multiply with :
       ## 0.01 (cm2m) * 1000 (water density, kg m-3) * latent heat of vaporization (J kg-1)
       ## latent heat of vaporization is not constant and it varies slightly with temperature, get.lv() returns 2.5e6 J kg-1 by default
       output[[12]] <- (sub.sipnet.output$npp * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s  # Qle W/m2
+    } else {
+     output[[12]] <- (sub.sipnet.output$evapotranspiration * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s  # Qle W/m2
     }
     output[[13]] <- (sub.sipnet.output$fluxestranspiration * 10) / timestep.s  # Transpiration kgW/m2/s
     output[[14]] <- (sub.sipnet.output$soilWater * 10)  # Soil moisture kgW/m2

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

@@ -491,9 +491,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       param[which(param[, 1] == "microbeInit"), 2] <- IC$microbe
     }
   }
+
   else if (length(settings$run$inputs$poolinitcond$path)>0) {
     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))
 
     if(!is.null(IC.pools)){
@@ -503,7 +505,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
         param[which(param[, 1] == "plantWoodInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$wood, "g m-2", "g m-2")
       }
       ## laiInit m2/m2
-      lai <- try(ncdf4::ncvar_get(IC.nc,"LAI"),silent = TRUE)
+      lai <- IC.pools$LAI
       if (!is.na(lai) && is.numeric(lai)) {
         param[which(param[, 1] == "laiInit"), 2] <- lai
       }

models/sipnet/R/write_restart.SIPNET.R

@@ -42,7 +42,7 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
                 file.path(outdir, runid, paste0("sipnet.", as.Date(start.time), ".out")))
     system(paste("rm", file.path(rundir, runid, "sipnet.clim")))
   } else {
-    print(paste("Files not renamed -- Need to rerun year", start.time, "before next time step"))
+    print(paste("Files not renamed -- Need to rerun timestep", start.time, "before next time step"))
   }
 
   settings$run$start.date <- start.time

modules/assim.sequential/DESCRIPTION

@@ -53,3 +53,4 @@ License: BSD_3_clause + file LICENSE
 Copyright: Authors
 Encoding: UTF-8
 RoxygenNote: 7.1.2
+

modules/assim.sequential/NAMESPACE

@@ -22,6 +22,7 @@ export(alltocs)
 export(alr)
 export(assessParams)
 export(block_matrix)
+export(build_X)
 export(conj_wt_wishart_sampler)
 export(dwtmnorm)
 export(generate_colors_sda)
@@ -45,6 +46,7 @@ export(sampler_toggle)
 export(sda.enkf)
 export(sda.enkf.multisite)
 export(sda.enkf.original)
+export(sda_matchparam)
 export(simple.local)
 export(tobit.model)
 export(tobit2space.model)

modules/assim.sequential/R/Analysis_sda_multiSite.R

@@ -377,7 +377,7 @@ GEF.MultiSite<-function(setting, Forecast, Observed, H, extraArg,...){
     dimensions.tobit = list(X = length(elements.W.Data),
                             X.mod = ncol(X),
                             Q = c(nrow(aqq), ncol(aqq))
-                            )
+    )
 
     # Contants defined in the model
     constants.tobit <-
@@ -402,13 +402,35 @@ GEF.MultiSite<-function(setting, Forecast, Observed, H, extraArg,...){
         r = solve(R)
       )
 
-  # This is the first step in making the nimble model - Nimble does some preliminary checks on the code    
-    model_pred <- nimbleModel(GEF.MultiSite.Nimble,
-                              data = data.tobit,
-                              dimensions = dimensions.tobit,
-                              constants = constants.tobit,
-                              inits = inits.pred,
-                              name = 'base')
+  # This is the first step in making the nimble model - Nimble does some preliminary checks on the code 
+    #special case for YN == 1 to run nimble model without for loops around nH
+    if(constants.tobit$YN == 1){
+      #add error message if trying to run SDA with 1 obs and 1 state variable no model currently exists to handle this case, need to remove for loop from GEF_singleobs_nimble for this case and save new model
+      if(constants.tobit$N == 1){
+        PEcAn.logger::logger.error("No model exists for assimilating 1 observation and 1 state variable, add more state variables or edit GEF_singleobs_nimble to work with 1 state variable")
+      }
+      #slight adjustment to inputs for nimble function when running with 1 obs
+      inits.pred$qq <- 0.368
+      dimensions.tobit$y.censored <- 1
+      dimensions.tobit$y.ind <- 1
+      constants.tobit$q.type <- NULL
+
+      model_pred <- nimbleModel(GEF_singleobs_nimble,
+                                data = data.tobit,
+                                dimensions = dimensions.tobit,
+                                constants = constants.tobit,
+                                inits = inits.pred,
+                                name = 'base')
+    }else{
+      model_pred <- nimbleModel(GEF.MultiSite.Nimble,
+                                data = data.tobit,
+                                dimensions = dimensions.tobit,
+                                constants = constants.tobit,
+                                inits = inits.pred,
+                                name = 'base')
+    }
+
+
 
 
     model_pred$initializeInfo()
@@ -418,20 +440,23 @@ GEF.MultiSite<-function(setting, Forecast, Observed, H, extraArg,...){
     conf$addMonitors(c("X","Xall","q","Xs")) 
     samplerNumberOffset <<- length(conf$getSamplers())
 
-    for(i in 1:length(y.ind)) {
-      node <- paste0('y.censored[',i,']')
-      conf$addSampler(node, 'toggle', control=list(type='RW'))
-    }
-
+      for(i in 1:length(y.ind)) {
+        node <- paste0('y.censored[',i,']')
+        conf$addSampler(node, 'toggle', control=list(type='RW'))
+      }
+
+
     conf$printSamplers()
 
     Rmcmc <<- buildMCMC(conf)
     Cmodel <<- compileNimble(model_pred)
-    Cmcmc <<- compileNimble(Rmcmc, project = model_pred)
+    Cmcmc <<- compileNimble(Rmcmc, project = model_pred, showCompilerOutput = TRUE)
 
     for(i in 1:length(y.ind)) {
-      valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 1-y.ind[i])
-    }
+        valueInCompiledNimbleFunction(Cmcmc$samplerFunctions[[samplerNumberOffset+i]], 'toggle', 1-y.ind[i])
+      }
+
+
 
     save(
       inits.pred,

modules/assim.sequential/R/Nimble_codes.R

@@ -331,3 +331,37 @@ conj_wt_wishart_sampler <-  nimbleFunction(
     }
   )
 )
+
+GEF_singleobs_nimble <-  nimbleCode({
+
+  # Sorting out qs
+  qq ~ dgamma(aq, bq) ## aq and bq are estimated over time
+  q[1, 1] <- qq * diag(YN)
+
+  # # X model
+  X.mod[1:N] ~ dmnorm(mean = muf[1:N], cov = pf[1:N, 1:N])
+  # # got rid of for loop no need when nH = 1
+  Xs[1] <- X.mod[H]
+
+  ## add process error to x model but just for the state variables that we have data and H knows who
+  #changed model from dmnorm to dnorm to accomodate when only assimilating 1 obs
+  X[1]  ~ dnorm(Xs[1], q[1, 1])
+
+  ## Likelihood
+  #changed model from dmnorm to dnorm to accomodate when only assimilating 1 obs
+  y.censored[1] ~ dnorm(X[1], r[1, 1])
+
+  #puting the ones that they don't have q in Xall - They come from X.model
+  # If I don't have data on then then their q is not identifiable, so we use the same Xs as Xmodel
+  for (j in 1:nNotH) {
+    tmpXmod[j]  <- X.mod[NotH[j]]
+    Xall[NotH[j]] <- tmpXmod[j]
+  }
+  # # # #These are the one that they have data and their q can be estimated
+  #got rid of for loop no need when nH = 1
+  Xall[H]  <- X[1]
+
+  y.ind[1] ~ dinterval(y.censored[1], 0)
+
+
+})

modules/assim.sequential/R/build_X.R

@@ -0,0 +1,95 @@
+#' build_X
+#' 
+#' @name build_X
+#' @author Alexis Helgeson
+#' 
+#' @description builds X matrix for SDA
+#'
+#' @param new.params object created from sda_matchparam, passed from sda.enkf_MultiSite
+#' @param nens number of ensemble members i.e. runs
+#' @param read_restart_times passed from sda.enkf_MultiSite
+#' @param settings settings object, passed from sda.enkf_MultiSite
+#' @param outdir location of previous run output folder containing .nc files
+#' @param out.configs object created for build_X passed from sda.enkf_MultiSite
+#' @param t Default t=1, for function to work within time loop
+#' @param var.names list of state variables taken from settings object
+#' @param my.read_restart object that points to the model restart function i.e. read_restart.SIPNET
+#'
+#' @return X ready to be passed to SDA Analysis code
+#' @export
+#'
+#' @examples
+build_X <- function(out.configs, settings, new.params, nens, read_restart_times, outdir, t = 1, var.names, my.read_restart){
+
+  if(t == 1){
+    reads <-
+      furrr::future_pmap(list(out.configs %>% `class<-`(c("list")), settings, new.params),function(configs,settings,siteparams) {
+        # Loading the model package - this is required bc of the furrr
+        #library(paste0("PEcAn.",settings$model$type), character.only = TRUE)
+        #source("~/pecan/models/sipnet/R/read_restart.SIPNET.R")
+
+        X_tmp <- vector("list", 2)
+
+        for (i in seq_len(nens)) {
+          X_tmp[[i]] <- do.call( my.read_restart,
+                                 args = list(
+                                   outdir = outdir,
+                                   runid = settings$runs$id[i] %>% as.character(),
+                                   stop.time = read_restart_times[t+1],
+                                   settings = settings,
+                                   var.names = var.names,
+                                   params = siteparams[[i]]
+                                 )
+          )
+
+        }
+        return(X_tmp)
+      })
+
+  }else{
+    reads <-
+      furrr::future_pmap(list(out.configs %>% `class<-`(c("list")), settings, new.params),function(configs,settings,siteparams) {
+        # Loading the model package - this is required bc of the furrr
+        #library(paste0("PEcAn.",settings$model$type), character.only = TRUE)
+        #source("~/pecan/models/sipnet/R/read_restart.SIPNET.R")
+
+        X_tmp <- vector("list", 2)
+
+        for (i in seq_len(nens)) {
+          X_tmp[[i]] <- do.call( my.read_restart,
+                                 args = list(
+                                   outdir = outdir,
+                                   runid = configs$runs$id[i] %>% as.character(),
+                                   stop.time = read_restart_times[t+1],
+                                   settings = settings,
+                                   var.names = var.names,
+                                   params = siteparams[[i]]
+                                 )
+          )
+
+        }
+        return(X_tmp)
+      })
+  }
+
+  #commented out text below describes future work - Alexis
+  # #let's read the parameters of each site/ens
+  # params.list <- reads %>% map(~.x %>% map("params"))
+  # # Now let's read the state variables of site/ens
+  # X <- reads %>% map(~.x %>% map_df(~.x[["X"]] %>% t %>% as.data.frame))
+  # 
+  # # Now we have a matrix that columns are state variables and rows are ensembles.
+  # # this matrix looks like this
+  # #         GWBI    AbvGrndWood   GWBI    AbvGrndWood
+  # #[1,]  3.872521     37.2581  3.872521     37.2581
+  # # But therer is an attribute called `Site` which tells yout what column is for what site id - check out attr (X,"Site")
+  # if (multi.site.flag){
+  #   X <- X %>%
+  #     map_dfc(~.x) %>% 
+  #     as.matrix() %>%
+  #     `colnames<-`(c(rep(var.names, length(X)))) %>%
+  #     `attr<-`('Site',c(rep(site.ids, each=length(var.names))))
+  # }
+
+  return(reads)
+}