From 32a1267403ebfbac971a019983ebe7a76c8e22aa Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Tue, 20 Jun 2023 11:04:46 +0200
Subject: [PATCH] eval posterior
---
dartwrf/evaluate_obs_space.py | 50 +++++++++++++++++------------------
evaluate_plus1.py | 21 ++++++++-------
2 files changed, 36 insertions(+), 35 deletions(-)
diff --git a/dartwrf/evaluate_obs_space.py b/dartwrf/evaluate_obs_space.py
index 2726910..1d11971 100755
--- a/dartwrf/evaluate_obs_space.py
+++ b/dartwrf/evaluate_obs_space.py
@@ -1,40 +1,42 @@
-import os, sys, shutil, warnings
+import os, sys, shutil, warnings, glob
import datetime as dt
from dartwrf.exp_config import exp
from dartwrf.server_config import cluster
from dartwrf import assim_synth_obs as aso
-def evaluate_one_time(init, valid):
+def evaluate_one_time(assim_time, valid_time, use_other_obsseq=False):
"""Evaluate the ensemble forecast in observation space at a given time, apart from the analysis time.
Args:
- init (datetime): initialization time of the forecast
- valid (datetime): valid time of the forecast
+ assim_time (datetime): initialization time of the forecast
+ valid_time (datetime): valid_time of the forecast
Returns:
None
"""
+ aso.prepare_nature_dart(valid_time)
+ aso.prepare_prior_ensemble(valid_time, prior_init_time=assim_time, prior_valid_time=valid_time, prior_path_exp=cluster.archivedir)
- # # prepare nature and prior ensemble
- aso.prepare_nature_dart(valid)
- aso.prepare_prior_ensemble(valid, prior_init_time=init, prior_valid_time=valid, prior_path_exp=cluster.archivedir)
+ # does an observation exist at this time?
+ f_oso = cluster.archivedir+valid_time.strftime("/obs_seq_out/%Y-%m-%d_%H:%M:%S_obs_seq.out")
-
- if use_other_obsseq: # use a different obsseq file
- f_obs_seq_out = use_other_obsseq
- shutil.copy(f_obs_seq_out, cluster.dart_rundir+'/obs_seq.out')
-
- else: # from same exp
- # use the last assimilation obsseq file for the observation locations (note: observed values are not valid)
- f_obs_seq_out = cluster.archivedir+valid.strftime("/obs_seq_out/%Y-%m-%d_%H:%M_obs_seq.out")
+ if os.path.exists(f_oso):
+ # use the existing file
+ shutil.copy(f_oso, cluster.dart_rundir+'/obs_seq.out')
+ else:
+ # use an old obsseq file and
+ f_oso = cluster.archivedir+assim_time.strftime("/obs_seq_out/%Y-%m-%d_%H:%M:%S_obs_seq.out")
from dartwrf.obs import obsseq
- oso = obsseq.ObsSeq(f_obs_seq_out)
- oso.df['observations'] = -9999
+ oso = obsseq.ObsSeq(f_oso)
+
+ # overwrite obs/truth values with "missing value"
+ oso.df['observations'] = -888888.0
+ oso.df['truth'] = -888888.0
oso.to_dart(cluster.dart_rundir+'/obs_seq.out')
- aso.evaluate(valid, output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate")
+ aso.evaluate(valid_time, output_format="%Y-%m-%d_%H:%M:%S_obs_seq.final-evaluate")
if __name__ == "__main__":
@@ -45,15 +47,13 @@ if __name__ == "__main__":
Usage: python3 evaluate_obs_space.py init1,valid1 init2,valid2 ...
"""
args = sys.argv[1:]
- init_valid_tuples = [a.split(',') for a in args]
-
- use_other_obsseq = False
+ arg_tuples = [a.split(',') for a in args]
# we need an existing run_DART folder
aso.prepare_run_DART_folder()
- for (init, valid) in init_valid_tuples:
+ for (assim_time, valid_time) in arg_tuples:
- init = dt.datetime.strptime(init, "%Y-%m-%d_%H:%M")
- valid = dt.datetime.strptime(valid, "%Y-%m-%d_%H:%M")
- evaluate_one_time(init, valid)
\ No newline at end of file
+ assim_time = dt.datetime.strptime(assim_time, "%Y-%m-%d_%H:%M")
+ valid_time = dt.datetime.strptime(valid_time, "%Y-%m-%d_%H:%M:%S")
+ evaluate_one_time(assim_time, valid_time)
\ No newline at end of file
diff --git a/evaluate_plus1.py b/evaluate_plus1.py
index b128c3b..f7f772f 100755
--- a/evaluate_plus1.py
+++ b/evaluate_plus1.py
@@ -6,18 +6,19 @@ import datetime as dt
from dartwrf.workflows import WorkFlows
w = WorkFlows(exp_config='exp_nonlin.py', server_config='jet.py')
+id = None
-
-assim_times = [dt.datetime(2008,7,30,12),
- dt.datetime(2008,7,30,12,30),
- dt.datetime(2008,7,30,13),
- dt.datetime(2008,7,30,13,30),
- dt.datetime(2008,7,30,14),]
+assim_times = [dt.datetime(2008,7,30,12), ]
+ # dt.datetime(2008,7,30,12,30),
+ # dt.datetime(2008,7,30,13),
+ # dt.datetime(2008,7,30,13,30),
+ # dt.datetime(2008,7,30,14),]
# generate observations at +1 minute after the assimilation time
-obs_times = [each+dt.timedelta(minutes=1) for each in assim_times]
-w.generate_obsseq_out(obs_times)
+tuples = []
+for init in assim_times:
+ for s in range(30,3*60+1,30):
+ tuples.append((init, init+dt.timedelta(seconds=s)))
-# evaluate the forecast at +1 minute after the assimilation time
-w.evaluate_plus1(assim_times)
\ No newline at end of file
+w.evaluate_obs_posterior_after_analysis(tuples, depends_on=id)
--
GitLab