eval posterior

dartwrf/evaluate_obs_space.py

-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

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)