From 7f707ed1a7f813b51709d21cc44e12fa8dc0b144 Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Wed, 31 May 2023 11:26:14 +0200
Subject: [PATCH] evaluate +1m posterior

---
 dartwrf/evaluate_obs_space.py | 60 ++++++++++++++++++++---------------
 dartwrf/workflows.py          | 10 ++++++
 evaluate_plus1.py             | 23 ++++++++++++++
 3 files changed, 68 insertions(+), 25 deletions(-)
 create mode 100755 evaluate_plus1.py

diff --git a/dartwrf/evaluate_obs_space.py b/dartwrf/evaluate_obs_space.py
index 0436d92..2726910 100755
--- a/dartwrf/evaluate_obs_space.py
+++ b/dartwrf/evaluate_obs_space.py
@@ -1,49 +1,59 @@
 import os, sys, shutil, warnings
-import time as time_module
 import datetime as dt
-import numpy as np
 
 from dartwrf.exp_config import exp
 from dartwrf.server_config import cluster
 from dartwrf import assim_synth_obs as aso
 
-def get_previous_obsseq_file(time):
-     oso_input = cluster.archivedir+'/obs_seq_out' + time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out-beforeQC")
+def evaluate_one_time(init, valid):
+     """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
+
+     Returns:
+          None
+     """
+
+     # # 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)
+
 
-     if not os.path.isfile(oso_input):  # fallback
-          oso_input = cluster.archivedir+'/obs_seq_out' + time.strftime("/%Y-%m-%d_%H:%M_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")
 
-     return oso_input
+          from dartwrf.obs import obsseq
+          oso = obsseq.ObsSeq(f_obs_seq_out)
+          oso.df['observations'] = -9999
+          oso.to_dart(cluster.dart_rundir+'/obs_seq.out')
 
+     aso.evaluate(valid, output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate")
 
 
 if __name__ == "__main__":
      """Evaluate the ensemble forecast in observation space at a given time, apart from the analysis time.
 
      Note: Observations are not assimilated. This is only for evaluation purposes.
-     """
 
-     init = dt.datetime.strptime(sys.argv[1], "%Y-%m-%d_%H:%M")
-     time = dt.datetime.strptime(sys.argv[2], "%Y-%m-%d_%H:%M")
+     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
 
      # we need an existing run_DART folder
      aso.prepare_run_DART_folder()
 
-     # # prepare nature and prior ensemble
-     aso.prepare_nature_dart(time)
-     aso.prepare_prior_ensemble(time, prior_init_time=init, prior_valid_time=time, prior_path_exp=cluster.archivedir)
-
-     # tell DART to use the prior as input
-     aso.write_list_of_inputfiles_prior()
-
-     if use_other_obsseq:  # use a different obsseq file
-          oso_input = use_other_obsseq
-     else:  # from same exp
-
-          # use the last assimilation obsseq file for the observation locations (note: observed values are not valid)     
-          oso_input = get_previous_obsseq_file(init)
-          shutil.copy(oso_input, cluster.dart_rundir+'/obs_seq.out')
+     for (init, valid) in init_valid_tuples:
 
-     aso.evaluate(time, output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate")
\ No newline at end of file
+          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
diff --git a/dartwrf/workflows.py b/dartwrf/workflows.py
index f7914d2..10e7e40 100644
--- a/dartwrf/workflows.py
+++ b/dartwrf/workflows.py
@@ -322,6 +322,16 @@ class WorkFlows(object):
         id = self.cluster.run_job("obsseq_netcdf", cfg_update={"time": "10", "mail-type": "FAIL,END"}, 
                 depends_on=[depends_on])
         return id
+    
+    def evaluate_plus1(self, list_assim_times, depends_on=None):
+        list_of_tuples = [(init, (init+dt.timedelta(minutes=1))) for init in list_assim_times]
+        arg = ' '.join([ttuple[0].strftime('%Y-%m-%d_%H:%M,')+ttuple[1].strftime('%Y-%m-%d_%H:%M') for ttuple in list_of_tuples])
+
+        cmd = self.cluster.python+' '+self.cluster.scripts_rundir+'/evaluate_obs_space.py '+arg
+        id = self.cluster.run_job(cmd, 'eval+1'+self.exp.expname, cfg_update={"ntasks": "12", "mem": "50G", "ntasks-per-node": "12", "ntasks-per-core": "2", 
+                                                                              "time": "15", "mail-type": "FAIL"}, 
+                depends_on=[depends_on])
+        return id
 
 
     def verify_sat(self, depends_on=None):
diff --git a/evaluate_plus1.py b/evaluate_plus1.py
new file mode 100755
index 0000000..b128c3b
--- /dev/null
+++ b/evaluate_plus1.py
@@ -0,0 +1,23 @@
+#!/usr/bin/python3
+"""
+Evaluate the forecast in observation space one minute after the analysis time.
+"""
+import datetime as dt
+from dartwrf.workflows import WorkFlows
+
+w = WorkFlows(exp_config='exp_nonlin.py', server_config='jet.py')
+
+
+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)
+
+
+# evaluate the forecast at +1 minute after the assimilation time
+w.evaluate_plus1(assim_times)
\ No newline at end of file
-- 
GitLab