Skip to content

Commit

Permalink
Merge pull request #3161 from DongchenZ/develop
Browse files Browse the repository at this point in the history
Update SDA_runner script.
  • Loading branch information
mdietze authored Apr 17, 2023
2 parents 1865831 + 150c420 commit 1043bfd
Showing 1 changed file with 90 additions and 36 deletions.
126 changes: 90 additions & 36 deletions modules/assim.sequential/inst/SDA_runner.R
Original file line number Diff line number Diff line change
@@ -1,39 +1,93 @@
#query_pfts
setwd("/projectnb/dietzelab/dongchen/All_NEON_SDA")
settings <- PEcAn.settings::read.settings("test_NEON2.xml")


#prepare obs
load("Obs_2012_2021/obs_cov.Rdata")
load("Obs_2012_2021/obs_mean.Rdata")

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(Kendall)

setwd("/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA")
settings_dir <- "/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/IC/pecan.xml"
settings <- PEcAn.settings::read.settings(settings_dir)
#prepare samples
get.parameter.samples(settings, ens.sample.method = settings$ensemble$samplingspace$parameters$method)

# PEcAn.uncertainty::get.parameter.samples(settings, ens.sample.method = settings$ensemble$samplingspace$parameters$method)
#prepare settings
new.settings <- PEcAn.settings::prepare.settings(settings)


#debug on sda.enkf_Multisite.R
settings <- new.settings

#run sda main function
sda.enkf.multisite(settings,
obs.mean,
obs.cov,
Q = NULL,
restart = FALSE,
forceRun = TRUE,
keepNC = TRUE,
control=list(trace = TRUE,
FF = FALSE,
interactivePlot = FALSE,
TimeseriesPlot = FALSE,
BiasPlot = FALSE,
plot.title = NULL,
facet.plots = FALSE,
debug = FALSE,
pause = FALSE,
Profiling = FALSE,
OutlierDetection=FALSE))
settings <- PEcAn.settings::prepare.settings(settings)
All_obs_prep <- settings$state.data.assimilation$Obs_Prep
folder_prefix <- c("AGB", "LAI", "SMP", "SoilC")
comb <- combn(4, 3)
for (i in 1:4) {
settings$state.data.assimilation$Obs_Prep <- All_obs_prep[c(comb[,i], 5:8)]
folder_name <- paste0(folder_prefix[comb[,i]], collapse = '_')
#prep obs
obs <- PEcAnAssimSequential::SDA_OBS_Assembler(settings = settings)
#prepare obs
load("/projectnb/dietzelab/dongchen/All_NEON_SDA/test_OBS/Rdata/obs.mean.Rdata")
load("/projectnb/dietzelab/dongchen/All_NEON_SDA/test_OBS/Rdata/obs.cov.Rdata")

for (i in 1:length(obs.mean)) {
if(is.null(obs.mean[[i]][[1]])){
next
}
for (j in 1:length(obs.mean[[i]])) {
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
}else{
if(obs.cov[[i]][[j]] == 0){
obs.cov[[i]][[j]] <- 1
}
}
}
}
load("/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA_ALL/enkf.Rdata")
PEcAn.utils::sendmail("[email protected]", "[email protected]", "SDA progress report", paste("Started working on", folder_name))
if(is.character(try(PEcAnAssimSequential::sda.enkf.multisite(settings = settings,
obs.mean = obs.mean,
obs.cov = obs.cov,
Q = NULL,
restart = FALSE,
forceRun = TRUE,
keepNC = TRUE,
pre_enkf_params = enkf.params,
control=list(trace = TRUE,
FF = FALSE,
interactivePlot = FALSE,
TimeseriesPlot = F,
BiasPlot = FALSE,
plot.title = NULL,
facet.plots = FALSE,
debug = FALSE,
pause = FALSE,
Profiling = FALSE,
OutlierDetection=FALSE)), silent = T))){
PEcAn.utils::sendmail("[email protected]", "[email protected]", "SDA running error", "Error")
break
}else{
PEcAn.utils::sendmail("[email protected]", "[email protected]", "SDA progress report", paste0(folder_name, " has been completed!"))
}
file.rename('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA',
file.path('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42',
paste0("SDA_", folder_name)))
dir.create('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA')
dir.create('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA/out')
dir.create('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA/run')
file.copy(file.path('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42',
paste0("SDA_", folder_name),
"samples.Rdata"),
file.path('/projectnb/dietzelab/dongchen/All_NEON_SDA/NEON42/SDA',
"samples.Rdata"))
}

0 comments on commit 1043bfd

Please sign in to comment.