various fixes

scripts/apply_obs_op_dart.py

@@ -7,56 +7,13 @@ from config.cfg import exp, cluster
 from utils import symlink, copy, mkdir, sed_inplace, append_file, print
 import create_obsseq as osq
 import assim_synth_obs as aso
-from assim_synth_obs import read_prior_obs, set_DART_nml, generate_observations, assimilate
 import pre_assim
-
-def run_operator(obscfg, time):
-    """
-    time_for_dart (dt.datetime) : needs to be consistent with wrfout files!
-    """
-
-    # get observation file (obs not important, but their locations)
-    # this should correspond to configuration to have same locations as in real assim
-    os.chdir(cluster.dartrundir)
-
-    n_obs = obscfg['n_obs']
-    error_var = (obscfg['err_std'])**2
-    sat_channel = obscfg.get('sat_channel', False)
-    cov_loc = obscfg['cov_loc_radius_km']
-    dist_obs = obscfg.get('distance_between_obs_km', False)
-
-    obs_coords = osq.calc_obs_locations(n_obs, coords_from_domaincenter=False, 
-                                        distance_between_obs_km=dist_obs, 
-                                        fpath_obs_locations=None)
-    osq.sat(time, sat_channel, obs_coords, error_var,
-            output_path=cluster.dartrundir)
-    assert os.path.exists(cluster.dartrundir + '/obs_seq.in')
-
-    # prepare dummy nature
-    os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
 import wrfout_add_geo
-    wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', 
-                        cluster.dartrundir+'/wrfout_d01')
 
-    print('running perfect model obs')
-    os.system('mpirun -np 12 ./perfect_model_obs')
+"""Apply observation operator to some ensemble state
+i.e. wrfout files in an archive directory
 
-    # set namelist for filter (calc only forward op)
-    aso.set_DART_nml(sat_channel=sat_channel,
-                  just_prior_values=True)
-
-    # run filter
-    assert os.path.exists(cluster.dartrundir+'/obs_seq.out')
-    print('running filter')
-    os.system('mpirun -np 40 ./filter')
-
-    # copy output to archive
-    savedir = cluster.archivedir()+'/obs_seq_final_1min/'
-    mkdir(savedir)
-    obsname = obscfg['kind']
-    
-    copy(cluster.dartrundir+'/obs_seq.final', savedir+fout)
-    print('output of observation operator saved to', fout)
+"""
 
 
 if __name__ == '__main__':
@@ -66,6 +23,7 @@ if __name__ == '__main__':
     print(prev_forecast_init, time)
 
     # link ensemble states to run_DART directory
+    # we want the observation operator applied to these states!
     pre_assim.run(time, prev_forecast_init, exppath_firstguess)
 
     savedir = cluster.archivedir()+'/obs_seq_final_1min/'
@@ -85,9 +43,16 @@ if __name__ == '__main__':
                      just_prior_values=True)
 
         osq.create_obsseq_in(time, obscfg)
-        aso.generate_observations()
+
+        # prepare dummy nature (this Hx is irrelevant)
+        os.chdir(cluster.dartrundir)
+        os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
+        wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
+                           cluster.dartrundir+'/wrfout_d01')
+        aso.run_perfect_model_obs()
         aso.assimilate()
 
         archive_stage = savedir+kind
+        # only the prior state values are of interest in this file
         aso.archive_diagnostics(archive_stage, time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final'))
 

scripts/assim_synth_obs.py

@@ -67,27 +67,6 @@ def read_obsseqout(f):
             obs.append(observed)
     return true, obs
 
-# def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
-#     """
-#     overwrite observation errors in a obs_seq.out file
-#     according to the values in OEs
-#     """
-#     # write to txt (write whole obs_seq.out again)
-#     obsseq = open(fpath_obsseqin, 'r').readlines()
-#     obsseq_new = obsseq.copy()
-#     i_obs = 0
-#     for i, line in enumerate(obsseq):
-#         if 'kind\n' in line:
-#             i_line_oe = i+9  # 9 for satellite obs
-#             obsseq_new[i_line_oe] = ' '+str(OEs[i_obs])+'   \n'
-#             i_obs += 1
-
-#     os.rename(fpath_obsseqin, fpath_obsseqin+'-bak')  # backup
-#     # write cloud dependent errors (actually whole file)
-#     with open(fpath_obsseqin, 'w') as f:
-#         for line in obsseq_new:
-#             f.write(line)
-
 
 def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
                  just_prior_values=False):
@@ -95,7 +74,7 @@ def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
     cov_loc_radian = cov_loc_radius_km/earth_radius_km
     
     if just_prior_values:
-        template = cluster.scriptsdir+'/../templates/input.prioronly.nml'
+        template = cluster.scriptsdir+'/../templates/input.eval.nml'
     else:
         template = cluster.scriptsdir+'/../templates/input.nml'
     copy(template, cluster.dartrundir+'/input.nml')
@@ -128,6 +107,7 @@ def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
         append_file(cluster.dartrundir+'/input.nml', rttov_nml)
 
 def obs_operator_ensemble():
+    # assumes that prior ensemble is already linked to advance_temp<i>/wrfout_d01
     print('running obs operator on ensemble forecast')
     os.chdir(cluster.dartrundir)
 
@@ -139,12 +119,12 @@ def obs_operator_ensemble():
             # ens members are already linked to advance_temp<i>/wrfout_d01
             copy(cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01',
                  cluster.dartrundir+'/wrfout_d01')
-            
-            wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
-
             # DART may need a wrfinput file as well, which serves as a template for dimension sizes
             symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
             
+            # add geodata
+            wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
+
             # run perfect_model obs (forward operator)
             os.system('mpirun -np 12 ./perfect_model_obs > /dev/null')
 
@@ -160,25 +140,26 @@ def obs_operator_ensemble():
     else:
         raise NotImplementedError()
 
-def obs_operator_nature():
+def obs_operator_nature(time):
     print('running obs operator on nature run')
-    prepare_nature_dart()
-
-    os.chdir(cluster.dartrundir)
-    os.remove(cluster.dartrundir+'/obs_seq.out')
-    os.system('mpirun -np 12 ./perfect_model_obs')
+    prepare_nature_dart(time)
+    run_perfect_model_obs()
     true, _ = read_obsseqout(cluster.dartrundir+'/obs_seq.out')
     return true
 
-def prepare_nature_dart():
+
+def link_nature_to_dart_truth(time):
     # get wrfout_d01 from nature run
     shutil.copy(time.strftime(cluster.nature_wrfout),
                 cluster.dartrundir+'/wrfout_d01')
+    # DART may need a wrfinput file as well, which serves as a template for dimension sizes
+    symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
 
+
+def prepare_nature_dart(time):
+    link_nature_to_dart_truth(time)
     wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
 
-    # DART may need a wrfinput file as well, which serves as a template for dimension sizes
-    symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
 
 def calc_obserr_WV73(Hx_nature, Hx_prior):
 
@@ -196,8 +177,7 @@ def calc_obserr_WV73(Hx_nature, Hx_prior):
         OEs[iobs] = oe_nature
     return OEs
 
-def generate_observations():
-    print('generate actual observations')
+def run_perfect_model_obs():
     os.chdir(cluster.dartrundir)
     try_remove(cluster.dartrundir+'/obs_seq.out')
     if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
@@ -212,16 +192,17 @@ def assimilate():
         raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir)
     os.system('mpirun -np 48 ./filter')
 
-def archive_diagnostics(archive_stage, fname_final):
+def archive_diagnostics(archive_dir, time):
     print('archive obs space diagnostics')
-    mkdir(archive_stage)
-    copy(cluster.dartrundir+'/obs_seq.final', archive_stage+'/'+fname_final)
+    mkdir(archive_dir)
+    copy(cluster.dartrundir+'/obs_seq.final', 
+         archive_dir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final'))
 
-    try:
-        print('archive regression diagnostics')
-        copy(cluster.dartrundir+'/reg_diagnostics', archive_stage+'/reg_diagnostics')
-    except Exception as e:
-        warnings.warn(str(e))
+    # try:  # what are regression diagnostics?!
+    #     print('archive regression diagnostics')
+    #     copy(cluster.dartrundir+'/reg_diagnostics', archive_dir+'/reg_diagnostics')
+    # except Exception as e:
+    #     warnings.warn(str(e))
 
 def recycle_output():
     print('move output to input')
@@ -289,7 +270,7 @@ if __name__ == "__main__":
 
             osq.create_obsseq_in(time, obscfg, zero_error=True)  # zero error to get truth vals
 
-            Hx_nat = obs_operator_nature() 
+            Hx_nat = obs_operator_nature(time) 
             Hx_prior = obs_operator_ensemble()  # files are already linked to DART directory
             
             obscfg['err_std'] = calc_obserr_WV73(Hx_nat, Hx_prior)
@@ -297,10 +278,12 @@ if __name__ == "__main__":
             obscfg['err_std'] = np.ones(n_obs) * obscfg['err_std']
 
         osq.create_obsseq_in(time, obscfg)  # now with correct errors
-        generate_observations()
+        prepare_nature_dart(time)
+        run_perfect_model_obs()
 
         assimilate()
-        archive_diagnostics(archive_stage, '/obs_seq.final')
+        dir_obsseq = cluster.archivedir()+'/obs_seq_final/assim_stage'+str(istage)+'_'+kind
+        archive_diagnostics(dir_obsseq, time)
 
         if istage < n_stages-1:
             # recirculation: filter output -> input
@@ -311,7 +294,7 @@ if __name__ == "__main__":
         elif istage == n_stages-1:
             # last assimilation, continue integration now
             copy(cluster.dartrundir+'/input.nml', archive_stage+'/input.nml')
-            pass  # call update wrfinput from filteroutput later
+            archive_output_mean(archive_stage)
 
         else:
             RuntimeError('this should never occur?!')

scripts/prepare_namelist.py

 import os, sys, shutil, warnings
 import datetime as dt
+from config.cfg import exp, cluster
 
-sys.path.append(os.getcwd())
+sys.path.append(cluster.scriptsdir)
 from config.cfg import exp, cluster
 from utils import sed_inplace, copy, symlink, mkdir
 
-def run(cluster, iens, begin, end, hist_interval=5, radt=5):
+def run(iens, begin, end, hist_interval=5, radt=5, archive=True):
+    """
+    Args:
+    archive (bool): if True, write to archivedir of experiment
+        if False, write to WRF run directory
+    """
     rundir = cluster.wrf_rundir(iens)
     print(rundir)
     copy(cluster.namelist, rundir+'/namelist.input')
@@ -15,10 +21,14 @@ def run(cluster, iens, begin, end, hist_interval=5, radt=5):
     sed_inplace(rundir+'/namelist.input', '<hist_interval>', str(int(hist_interval)))
     sed_inplace(rundir+'/namelist.input', '<radt>', str(int(radt)))
 
+    if archive:
         archdir = cluster.archivedir()+begin.strftime('/%Y-%m-%d_%H:%M/'+str(iens)+'/')
+        os.makedirs(archdir, exist_ok=True)
+    else:
+        archdir = './'
     print('namelist for run from', begin, end, 'output to', archdir)
     sed_inplace(rundir+'/namelist.input', '<archivedir>', archdir)
-    os.makedirs(archdir, exist_ok=True)
+    
 
     # set times
     for k, v in {'<y1>': '%Y', '<m1>': '%m', '<d1>': '%d',
@@ -29,6 +39,7 @@ def run(cluster, iens, begin, end, hist_interval=5, radt=5):
         sed_inplace(rundir+'/namelist.input', k, end.strftime(v))
 
     #########################
+    if archive:
         try:
             print('copy wrfinput of this run to archive')
             wrfin_old = rundir+'/wrfinput_d01'
@@ -48,4 +59,4 @@ if __name__ == '__main__':
 
     print('prepare namelists for all ens members')
     for iens in range(1, exp.n_ens+1):
-        run(cluster, iens, begin, end, hist_interval=intv, radt=radt)
+        run(iens, begin, end, hist_interval=intv, radt=radt)

scripts/run_obs_diag.py

@@ -2,31 +2,40 @@ import os, sys, shutil, glob
 from config.cfg import exp, cluster
 from utils import symlink, copy, sed_inplace, append_file
 
-def run(folder_obs_seq_final):
 rundir_program = '/home/fs71386/lkugler/data/DART-WRF/rundir/'
 
-    files = sorted(glob.glob(folder_obs_seq_final+'/*.final'))  # input for obs_diag program
+
+def prepare(obserr_iszero='.true.'):
+    copy(cluster.scriptsdir+'/../templates/input.eval.nml',
+            rundir_program+'/input.nml')
+    sed_inplace(rundir_program+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
+    sed_inplace(rundir_program+'/input.nml', '<zero_error_obs>', obserr_iszero)
+    sed_inplace(rundir_program+'/input.nml', '<horiz_dist_only>', '.false.')  # dummy
+    sed_inplace(rundir_program+'/input.nml', '<vert_norm_hgt>', '5000.0')  # dummy
+    
+    append_file(rundir_program+'/input.nml', cluster.scriptsdir+'/../templates/obs_def_rttov.VIS.nml')
+
+
+def write_input_filelist(filepaths):
     fpath = rundir_program+'/obsdiag_inputlist.txt'
     print('writing', fpath)
-
     if os.path.exists(fpath):
         os.remove(fpath)
 
     with open(fpath, 'w') as f:
-        for fin in files:
+        for fin in filepaths:
             f.write(fin)
             f.write('\n')
 
+
+def run_obsdiag(filepaths, f_out='./obsdiag.nc'):
+    write_input_filelist(filepaths)
+
     for obserr_iszero in ['.true.', '.false.']:
-        print('ensure correct input.nml')
-        copy(cluster.scriptsdir+'/../templates/input.prioronly.nml',
-             rundir_program+'/input.nml')
-        sed_inplace(rundir_program+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
-        sed_inplace(rundir_program+'/input.nml', '<zero_error_obs>', obserr_iszero)
-        append_file(rundir_program+'/input.nml', cluster.scriptsdir+'/../templates/obs_def_rttov.VIS.nml')
+        prepare(obserr_iszero=obserr_iszero)
 
-        # run obs_diag
-        print('running obs_diag program')
+        # run_allinoneplace obs_diag
+        print('------ running obs_diag program')
         os.chdir(rundir_program)
         symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
         try:
@@ -36,24 +45,29 @@ def run(folder_obs_seq_final):
         os.system('./obs_diag >& obs_diag.log')  # caution, this overwrites obs_seq_to_netcdf
 
         # move output to archive
-        outdir = '/'.join(folder_obs_seq_final.split('/')[:-1])
+        #outdir = outdir #'/'.join(folder_obs_seq_final.split('/')[:-1])
         if obserr_iszero == '.true.':
-            fout = '/obs_diag_wrt_truth.nc'   
+            fout = outdir+'/'+f_out[:-3]+'_wrt_truth.nc'   
         elif obserr_iszero == '.false.':
-            fout = '/obs_diag_wrt_obs.nc'
-        print('moving output to', outdir+fout)
-        copy(rundir_program+'/obs_diag_output.nc', outdir+fout)
+            fout = outdir+'/'+f_out[:-3]+'_wrt_obs.nc' 
+        shutil.move(rundir_program+'/obs_diag_output.nc', fout)
+        print(fout, 'saved.')
+
 
+def run_obs_seq_to_netcdf(filepaths, f_out='./obs_epoch.nc'):
 
-    print('running obs_seq_to_netcdf program')
-    shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', cluster.dart_srcdir+'/obs_seq_to_netcdf')
-    symlink(cluster.dart_srcdir+'/obs_seq_to_netcdf', rundir_program+'/obs_seq_to_netcdf')
+    write_input_filelist(filepaths)
+    print('------ running obs_seq_to_netcdf program')
+    shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', rundir_program+'/obs_seq_to_netcdf')
     os.system('./obs_seq_to_netcdf  >& obs_seq_to_netcdf.log')  # caution, overwrites its own binary?!
-    print('moving output to', outdir+'/obs_seq...')
-    os.system('mv '+rundir_program+'/obs_epoch_*.nc '+outdir+'/')
+    shutil.move(rundir_program+'/obs_epoch_001.nc', f_out)
+    print(f_out, 'saved.')
 
 
 if __name__ == '__main__':
-    #folder_obs_seq_final = '/home/fs71386/lkugler/data/sim_archive/exp_v1.11_LMU_filter2/obs_seq_final/'
+    #folder_obs_seq_final = '/home/fs71386/lkugler/data/DART-WRF/rundir/test' 
     folder_obs_seq_final = str(sys.argv[1])
-    run(folder_obs_seq_final)
+    files = sorted(glob.glob(folder_obs_seq_final+'/*.final'))  # input for obs_diag program
+    
+    run_obsdiag(files, outdir=folder_obs_seq_final)  # input must be files with posterior data!!
+    run_obs_seq_to_netcdf(files, outdir=folder_obs_seq_final)  # input can be files without posterior data

templates/input.prioronly.nml

@@ -31,7 +31,7 @@
    async                    = 0,
    adv_ens_command          = "../shell_scripts/advance_model.csh",
    ens_size                 = <n_ens>,
-   obs_sequence_in_name     = "obs_seq_all.out",
+   obs_sequence_in_name     = "obs_seq.out",
    obs_sequence_out_name    = "obs_seq.final",
    input_state_file_list    = "input_list.txt"
    output_state_file_list   = "output_list.txt"
@@ -143,7 +143,8 @@
    /
 
 &obs_kind_nml
-   assimilate_these_obs_types = 'RADIOSONDE_TEMPERATURE',
+   assimilate_these_obs_types = ,
+   evaluate_these_obs_types   = 'RADIOSONDE_TEMPERATURE',
                                 'RADIOSONDE_U_WIND_COMPONENT',
                                 'RADIOSONDE_V_WIND_COMPONENT',
                                 'SAT_U_WIND_COMPONENT',
@@ -160,7 +161,6 @@
                                 'MSG_4_SEVIRI_RADIANCE',
                                 'MSG_4_SEVIRI_TB',
                                 'MSG_4_SEVIRI_BDRF'
-   evaluate_these_obs_types   = 'RADIOSONDE_SPECIFIC_HUMIDITY',
    /
 
 # Notes for obs_def_radar_mod_nml:
@@ -246,7 +246,7 @@
 &location_nml
    horiz_dist_only             = .false.,
    vert_normalization_pressure = 6666666.7,
-   vert_normalization_height   = 5000000.0,
+   vert_normalization_height   = 500.0,
    vert_normalization_level    = 2666.7,
    vert_normalization_scale_height = 10.0,
    approximate_distance        = .false.,