Update NEE and LE gap-fill functions

base/utils/R/read.output.R

@@ -77,9 +77,9 @@ read.output <- function(runid, outdir,
   }
 
   # create list of *.nc years - look only for files formatted as
-  # YYYY.nc, the default pecan output file name standard
+  # YYYY.nc or YYYY-MM-DD.nc
   if (is.null(ncfiles)) {
-    ncfiles_sub <- list.files(path = outdir, pattern = "^-?[[:digit:]]{4}\\.nc$", full.names = FALSE)
+    ncfiles_sub <- list.files(path = outdir, pattern = "^-?[[:digit:]]{4}(-[[:digit:]]{2}-[[:digit:]]{2})?\\.nc$", full.names = FALSE)
     ncfiles <- file.path(outdir, ncfiles_sub)
   } else {
     # Assume the NetCDF files follow the PEcAn standard format
@@ -128,7 +128,7 @@ read.output <- function(runid, outdir,
 
   # Deduce file years from their names
   nc_years <- suppressWarnings(
-    as.numeric(gsub("^(-?[[:digit:]]{4})\\.nc", "\\1", ncfiles_sub))
+    as.numeric(gsub("^(-?[[:digit:]]{4})(-[[:digit:]]{2}-[[:digit:]]{2})?\\.nc$", "\\1", ncfiles_sub))
   )
   if (any(is.na(nc_years))) {
     PEcAn.logger::logger.debug(

docker/depends/pecan_package_dependencies.csv

@@ -33,6 +33,7 @@
 "data.table","*","base/utils","Suggests",FALSE
 "data.table","*","base/visualization","Imports",FALSE
 "data.table","*","models/biocro","Imports",FALSE
+"data.table","*","modules/data.land","Suggests",FALSE
 "data.table","*","modules/data.remote","Suggests",FALSE
 "dataone","*","modules/data.land","Suggests",FALSE
 "datapack","*","modules/data.land","Suggests",FALSE
@@ -41,6 +42,7 @@
 "dbplyr",">= 2.4.0","base/db","Imports",FALSE
 "devtools","*","models/ed","Suggests",FALSE
 "doParallel","*","modules/data.atmosphere","Suggests",FALSE
+"doParallel","*","modules/data.land","Suggests",FALSE
 "doParallel","*","modules/data.remote","Imports",FALSE
 "doSNOW","*","base/remote","Suggests",FALSE
 "doSNOW","*","base/utils","Suggests",FALSE
@@ -668,6 +670,7 @@
 "tibble","*","models/fates","Imports",FALSE
 "tibble","*","models/lpjguess","Imports",FALSE
 "tibble","*","modules/data.atmosphere","Imports",FALSE
+"tibble","*","modules/data.land","Suggests",FALSE
 "tibble","*","modules/data.remote","Suggests",FALSE
 "tibble","*","modules/meta.analysis","Suggests",FALSE
 "tictoc","*","modules/assim.sequential","Suggests",FALSE
@@ -723,6 +726,7 @@
 "withr","*","modules/data.land","Suggests",FALSE
 "withr","*","modules/uncertainty","Suggests",FALSE
 "xgboost","*","modules/assim.sequential","Suggests",FALSE
+"xgboost","*","modules/data.land","Suggests",FALSE
 "XML","*","base/workflow","Imports",FALSE
 "XML","*","models/biocro","Imports",FALSE
 "XML","*","models/maat","Imports",FALSE

local/sda_runner.R

@@ -0,0 +1,186 @@
+# loading libraries
+library(dplyr)
+library(xts)
+library(PEcAn.all)
+library(purrr)
+library(furrr)
+library(lubridate)
+library(nimble)
+library(ncdf4)
+library(PEcAnAssimSequential)
+library(dplyr)
+library(sp)
+library(raster)
+library(zoo)
+library(ggplot2)
+library(mnormt)
+library(sjmisc)
+library(stringr)
+library(doParallel)
+library(doSNOW)
+library(data.table)
+library(Kendall)
+library(lgarch)
+library(parallel)
+library(foreach)
+library(terra)
+setwd("/projectnb/dietzelab/guYANG/pecan/")
+## read settings xml file.
+load("/projectnb/dietzelab/guYANG/pecan/pecan.Rdata")
+
+### PEcAn Settings
+# settings_dir <- "/projectnb/dietzelab/guYANG/pecan/output/pecan_flux.xml"
+# settings <- PEcAn.settings::read.settings(settings_dir)
+
+# Change dir name and settings
+settings$outdir      <- "/projectnb/dietzelab/guYANG/pecan/output/5obs_monthly/"
+settings$rundir      <- file.path(settings$outdir, "run")
+settings$modeloutdir <- file.path(settings$outdir, "out")
+settings$host$rundir <- file.path(settings$outdir, "run")
+settings$host$outdir <- file.path(settings$outdir, "out")
+settings$host$folder <- file.path(settings$outdir, "out")
+settings$ensemble$size <- 25
+settings$state.data.assimilation$adjustment <- "FALSE"
+settings$host$prerun <- "module load R/4.4.0"
+###### Change Q type
+settings$state.data.assimilation$q.type <- "vector"
+settings$state.data.assimilation$aqq.Init <- "1"
+settings$state.data.assimilation$bqq.Init <- "1"
+## Fix the multi output in one timestep bug
+settings$model$jobtemplate <- "/projectnb/dietzelab/guYANG/pecan/runners/test7/sipnet_template.job"
+
+# Load the selected sites
+load("/projectnb/dietzelab/guYANG/pecan/runners/wishart_sda/sda_idx.Rdata")
+all_ids  <- vapply(settings, \(s) as.character(s$run$site$id), "")
+# Sub settings for this run
+settings <- settings[all_ids %in% keep_ids]
+settings <- PEcAn.settings::as.MultiSettings(settings)
+
+# setup the batch job settings.
+general.job <- list(cores = 28, folder.num = 80)
+batch.settings = structure(list(
+  general.job = general.job,
+  qsub.cmd = "qsub -l h_rt=24:00:00 -l mem_per_core=4G -l buyin -pe omp @CORES@ -V -N @NAME@ -o @STDOUT@ -e @STDERR@ -S /bin/bash"
+))
+settings$state.data.assimilation$batch.settings <- batch.settings
+
+# update settings with the actual PFTs.
+settings <- PEcAn.settings::prepare.settings(settings)
+
+# load obs
+# load("/projectnb/dietzelab/guYANG/pecan/runners/test10/obs.mean.RData")
+# load("/projectnb/dietzelab/guYANG/pecan/runners/test10/obs.cov.RData")
+
+load("/projectnb/dietzelab/guYANG/pecan/output/obs/Rdata/obs.mean.monthly.noSM.Rdata")
+load("/projectnb/dietzelab/guYANG/pecan/output/obs/Rdata/obs.cov.monthly.noSM.Rdata")
+
+sub_obs <- function(L, keep) setNames(lapply(L, \(l) l[names(l) %in% keep]), names(L))
+obs.mean <- sub_obs(obs.mean, keep_ids)
+obs.cov  <- sub_obs(obs.cov,  keep_ids)
+
+# replace zero observations and variances with small numbers.
+for (i in 1:length(obs.mean)) {
+  if(is.null(obs.mean[[i]][[1]])){
+    next
+  }
+  for (j in 1:length(obs.mean[[i]])) {
+    if (length(obs.mean[[i]][[j]])==0) {
+      next
+    }
+    obs.mean[[i]][[j]][which(obs.mean[[i]][[j]]==0)] <- 0.01
+    if(length(obs.cov[[i]][[j]]) > 1){
+      diag(obs.cov[[i]][[j]])[which(diag(obs.cov[[i]][[j]]<=0.1))] <- 0.1
+    }else{
+      if(obs.cov[[i]][[j]] <= 0.1){
+        obs.cov[[i]][[j]] <- 0.1
+      }
+    }
+  }
+}
+
+if (length(obs.cov[[i]][[j]]) > 1) {
+  d <- diag(obs.cov[[i]][[j]])
+  d[d <= 0.1] <- 0.1
+  diag(obs.cov[[i]][[j]]) <- d
+}
+
+# load PFT parameter file.
+# samples_src <- "/projectnb/dietzelab/guYANG/pecan/cali_samples.Rdata"
+samples_src <- "/projectnb/dietzelab/dongchen/anchorSites/NA_runs/SDA_8k_site/samples.Rdata"
+samples_dst <- file.path(settings$outdir, "samples.Rdata")
+dir.create(settings$outdir, recursive = TRUE, showWarnings = FALSE)
+if (!file.exists(samples_dst)) file.copy(samples_src, samples_dst, overwrite = TRUE)
+
+control = list(
+  TimeseriesPlot = FALSE,
+  OutlierDetection = FALSE,
+  send_email = NULL,
+  keepNC = TRUE,
+  forceRun = TRUE,
+  run_parallel = FALSE,
+  MCMC.args = list(
+    niter = 200000,   
+    nburnin = 100000,
+    nthin = 5,
+    nchain = 3
+  ),
+  merge_nc = FALSE,
+  execution = "qsub_parallel"   # or "qsub" / "local"
+)
+
+res <- PEcAnAssimSequential:::sda.enkf_local(
+  settings = settings,
+  obs.mean = obs.mean,
+  obs.cov = obs.cov,
+  Q = NULL,
+  pre_enkf_params = NULL,
+  ensemble.samples = NULL,
+  control = list(
+    TimeseriesPlot = FALSE,
+    OutlierDetection = FALSE,
+    send_email = NULL,
+    keepNC = TRUE,
+    forceRun = TRUE,
+    run_parallel = FALSE,
+    MCMC.args = list(
+      niter = 200000,   
+      nburnin = 100000,
+      nthin = 5,
+      nchain = 3
+    ),
+    merge_nc = FALSE,
+    execution = "qsub_parallel"   # or "qsub" / "local"
+  )
+)
+
+#### qsub
+# PEcAnAssimSequential:::qsub_sda(
+#   settings = settings,
+#   obs.mean = obs.mean,
+#   obs.cov = obs.cov,
+#   Q = NULL,
+#   pre_enkf_params = NULL,
+#   ensemble.samples = NULL,
+#   outdir = "/projectnb/dietzelab/guYANG/pecan/output/qsub",
+#   control = control
+# )
+
+# job_lines <- c(
+#   "#!/bin/bash",
+#   "module load R/4.4.0",
+#   "Rscript /projectnb/dietzelab/guYANG/pecan/pecan/local/sda_runner.R"
+# )
+# writeLines(job_lines, "/projectnb/dietzelab/guYANG/pecan/output/log/sda_runner_1.sh")
+
+# qsub -l h_rt=10:00:00 \
+# -l buyin \
+# -l mem_per_core=8G \
+# -pe omp 28 \
+# -V \
+# -N sda_runner_1 \
+# -o /projectnb/dietzelab/guYANG/pecan/output/log/sda_runner_1.out \
+# -e /projectnb/dietzelab/guYANG/pecan/output/log/sda_runner_1.err \
+# -M yanggu@bu.edu \
+# -m abe \
+# -S /bin/bash \
+# /projectnb/dietzelab/guYANG/pecan/output/log/sda_runner_1.sh