linear posterior

dartwrf/calc_linear_posterior.py

+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 assimilate as aso
+from dartwrf import dart_nml
+from dartwrf.obs import create_obsseq_out as osq_out
+
+pattern_osf_linear = cluster.archivedir + "/diagnostics/%Y-%m-%d_%H:%M_obs_seq.final-linear"  # how an obs_seq.final file is archived
+
+def calc_lin_posterior(assim_time):
+     """Run filter to get the linear, optimal-interpolation posterior in obs-space
+
+     Args:
+          prior_init_time (datetime): initialization of prior forecast
+          assim_time (datetime): assimilation time
+
+     Returns:
+          None
+     """
+     # what is the prior for this assim time?
+     # read link_to_prior
+     link_to_prior = cluster.archivedir + assim_time.strftime('/%Y-%m-%d_%H:%M/')+'link_to_prior.txt'
+     with open(link_to_prior) as f:
+          lines = [l.strip() for l in f.readlines()]
+          prior_exp = lines[0]
+          prior_init_time = dt.datetime.strptime(lines[1].strip('/'), '%Y-%m-%d_%H:%M')
+     
+     aso.prepare_prior_ensemble(assim_time, prior_init_time=prior_init_time, prior_valid_time=assim_time, prior_path_exp=prior_exp)
+     #aso.prepare_nature_dart(assim_time) 
+     
+     dart_nml.write_namelist()
+
+     # does an observation exist at this time?
+     #f_oso = assim_time.strftime(aso.pattern_obs_seq_out)
+     f_oso = cluster.archivedir+assim_time.strftime("/diagnostics/%Y-%m-%d_%H:%M/assim_stage0/obs_seq.out")
+
+     if os.path.exists(f_oso):
+          # use the existing file
+          shutil.copy(f_oso, cluster.dart_rundir+'/obs_seq.out')
+     else:
+          raise NotImplementedError(f_oso+' does not exist!')
+     
+     # link in the linear filter version
+     from dartwrf.utils import symlink
+     path_linear_filter = '/jetfs/home/lkugler/data/compile/DART/DART-10.8.3-kalmanposterior/models/wrf/work/filter'
+     symlink(path_linear_filter,
+             os.path.join(cluster.dart_rundir, 'filter'))
+
+     aso.filter(nproc=12)
+     aso.archive_filter_diagnostics(assim_time, pattern_osf_linear)
+
+
+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.
+
+     Usage: python3 evaluate_obs_space.py init1,valid1 init2,valid2 ...
+     """
+     args = sys.argv[1:]
+     #arg_tuples = [a.split(',') for a in args]
+
+     # we need an existing run_DART folder
+     aso.prepare_run_DART_folder()
+
+     for assim_time in args:
+          #prior_init_time = dt.datetime.strptime(prior_init_time, "%Y-%m-%d_%H:%M")
+          assim_time = dt.datetime.strptime(assim_time, "%Y-%m-%d_%H:%M")
+          calc_lin_posterior(assim_time)
+          
\ No newline at end of file

dartwrf/evaluate_obs_space.py

@@ -3,40 +3,49 @@ import datetime as dt
 
 from dartwrf.exp_config import exp
 from dartwrf.server_config import cluster
-from dartwrf import assim_synth_obs as aso
+from dartwrf import assimilate as aso
+from dartwrf import dart_nml
+from dartwrf.obs import create_obsseq_out as osq_out
 
 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:
           assim_time (datetime): initialization time of the forecast
-          valid_time (datetime): valid_time of the forecast
+          valid_time (datetime): valid_time of the forecast (e.g. +1 minute after assim_time)
 
      Returns:
           None
      """
-     aso.prepare_nature_dart(valid_time)
+     aso.prepare_run_DART_folder()
      aso.prepare_prior_ensemble(valid_time, prior_init_time=assim_time, prior_valid_time=valid_time, prior_path_exp=cluster.archivedir)
+     dart_nml.write_namelist()
 
      # does an observation exist at this time?
+     f_oso = valid_time.strftime(aso.pattern_obs_seq_out)
      f_oso = cluster.archivedir+valid_time.strftime("/obs_seq_out/%Y-%m-%d_%H:%M:%S_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_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_time, output_format="%Y-%m-%d_%H:%M:%S_obs_seq.final-evaluate")
+          if False:
+               # use an old obsseq file and overwrite obs/truth values with "missing value"
+               f_oso = cluster.archivedir+valid_time.strftime("/obs_seq_out/%Y-%m-%d_%H:%M_obs_seq.out")
+          
+               from dartwrf.obs import obsseq
+               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')
+          else:
+               # generate the observations for the specified valid_time
+               aso.prepare_nature_dart(valid_time)
+               osq_out.generate_obsseq_out(valid_time)
+     
+     aso.evaluate(valid_time, f_out_pattern=cluster.archivedir + "/diagnostics/%Y-%m-%d_%H:%M:%S_obs_seq.final-evaluate")
 
 
 if __name__ == "__main__":

tests/test_assim.py

 import os, shutil
-import datetime as dt
 
-from dartwrf import assim_synth_obs
-from dartwrf.obs import obsseq
-from dartwrf.server_config import cluster
+def test_dartwrf_cycled_da():
+    import datetime as dt
+    
+    from dartwrf.workflows import WorkFlows
+    w = WorkFlows(exp_config='exp_template.py', server_config='jet.py')
 
-class ExperimentConfiguration(object):
-    def __init__(self):
-        pass
+    timedelta_integrate = dt.timedelta(minutes=15)
+    timedelta_btw_assim = dt.timedelta(minutes=15)
 
-exp = ExperimentConfiguration()
-exp.expname = "test"
+    id = None
 
-obscfg_wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
-                kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=4,
-                error_generate=1., error_assimilate=False,
-                cov_loc_radius_km=32)
-exp.observations = [obscfg_wv73]
+    if True:  # random
+        # need some ensemble data to test
+        prior_path_exp = '/jetfs/scratch/lkugler/data/sim_archive/exp_v1.19_P2_noDA+1'
 
-time = dt.datetime(2008,7,30,12)
+        init_time = dt.datetime(2008, 7, 30, 12, 30)
+        time = dt.datetime(2008, 7, 30, 13)
+        last_assim_time = dt.datetime(2008, 7, 30, 13)
+        forecast_until = dt.datetime(2008, 7, 30, 13,10)
 
+        w.prepare_WRFrundir(init_time)
+        # id = w.run_ideal(depends_on=id)
+
+    prior_init_time = init_time
+
+    while time <= last_assim_time:
+
+        # usually we take the prior from the current time
+        # but one could use a prior from a different time from another run
+        # i.e. 13z as a prior to assimilate 12z observations
+        prior_valid_time = time
+
+        id = w.assimilate(time, prior_init_time, prior_valid_time, prior_path_exp, depends_on=id)
+        
+        # 1) Set posterior = prior
+        id = w.prepare_IC_from_prior(prior_path_exp, prior_init_time, prior_valid_time, depends_on=id)
+
+        # 2) Update posterior += updates from assimilation
+        id = w.update_IC_from_DA(time, depends_on=id)
+
+        # How long shall we integrate?
+        timedelta_integrate = timedelta_btw_assim
+        output_restart_interval = timedelta_btw_assim.total_seconds()/60
+        if time == last_assim_time:  # this_forecast_init.minute in [0,]:  # longer forecast every full hour
+            timedelta_integrate = forecast_until - last_assim_time  # dt.timedelta(hours=4)
+            output_restart_interval = 9999  # no restart file after last assim
+
+        # 3) Run WRF ensemble
+        id = w.run_ENS(begin=time,  # start integration from here
+                    end=time + timedelta_integrate,  # integrate until here
+                    output_restart_interval=output_restart_interval,
+                    first_minutes=False,
+                    depends_on=id)
+        
+        # as we have WRF output, we can use own exp path as prior
+        prior_path_exp = w.cluster.archivedir
+
+        time += timedelta_btw_assim
+        prior_init_time = time - timedelta_btw_assim
+    
 
 def test_overwrite_OE_assim():
+    import datetime as dt
+    from dartwrf import assimilate as aso
+    from dartwrf.obs import obsseq
+    from dartwrf.server_config import cluster
+
     # checks if modified entries are correctly written to DART files
-    input_template = cluster.scriptsdir + "/../tests/obs_seq.orig.out"
-    input_temporary = cluster.scriptsdir + "/../tests/obs_seq.out"
+    input_template = "./obs_seq.orig.out"
+    input_temporary = "./obs_seq.out"
     # TODO: MISSING FILE
-    # input_osf = cluster.scriptsdir + "/../tests/obs_seq.final"
-    output = cluster.scriptsdir + "/../tests/obs_seq.test.out"
+    # input_osf = "./obs_seq.final"
+    output = "./obs_seq.test.out"
     shutil.copy(input_template, input_temporary)
 
     oso = obsseq.ObsSeq(input_temporary)
-    osf = obsseq.ObsSeq(input_osf)
+    #osf = obsseq.ObsSeq(input_osf)
 
-    assim_synth_obs.set_obserr_assimilate_in_obsseqout(oso, osf, outfile=output)
+    aso.set_obserr_assimilate_in_obsseqout(oso, outfile=output)
 
     var_orig = oso.df['variance']
 
@@ -43,4 +88,6 @@ def test_overwrite_OE_assim():
     os.remove(input_temporary)
 
 if __name__ == '__main__':
-    test_overwrite_OE_assim()
+    # test_dartwrf_cycled_da()
+    # test_overwrite_OE_assim()
+    pass