use different values of obserr in obs generation and assimilation

config/cfg.py

@@ -9,42 +9,49 @@ class ExperimentConfiguration(object):
 
 
 exp = ExperimentConfiguration()
-exp.expname = "exp_v1.16_P1_40mem"
+exp.expname = "exp_v1.16_P1-1_Radar"
 exp.model_dx = 2000
-exp.timestep = 8
 exp.n_ens = 40
 exp.n_nodes = 10
 
 n_obs = 1600  # 50km res: 64:8x8; 10km res: 1600:40x40  # radar: n_obs for each observation height level
 
 vis = dict(plotname='VIS 0.6µm',  plotunits='[1]',
-           kind='MSG_4_SEVIRI_BDRF',
-           sat_channel=1, n_obs=n_obs, err_std=0.03,
+           kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs, 
+           err_std=0.03,
+           error_generate=0.03, error_assimilate=0.06,
            cov_loc_radius_km=32)
+
 wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
-            kind='MSG_4_SEVIRI_TB',
-            sat_channel=6, n_obs=n_obs, err_std=False,
+            kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs, 
+            err_std=False,
+            error_generate=False, error_assimilate=False,
             cov_loc_radius_km=32)
+
 ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
-             kind='MSG_4_SEVIRI_TB',
-             sat_channel=9, n_obs=n_obs, err_std=5.,
+             kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs, 
+             err_std=5.,
+             error_generate=5., error_assimilate=10.,
              cov_loc_radius_km=32)
 
 radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
-             kind='RADAR_REFLECTIVITY', 
-             n_obs=n_obs, err_std=5.,
+             kind='RADAR_REFLECTIVITY', n_obs=n_obs, 
+             error_generate=2.5, error_assimilate=5.,
              heights=np.arange(1000, 15001, 1000),
-             cov_loc_radius_km=30, cov_loc_vert_km=4)
+             cov_loc_radius_km=32, cov_loc_vert_km=4)
 
 t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
-           kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=0.1, 
+           kind='SYNOP_TEMPERATURE', n_obs=n_obs, 
+           error_generate=0.1, error_assimilate=1.,
            cov_loc_radius_km=20, cov_loc_vert_km=3)
+
 psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
-            kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
+            kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, 
+            error_generate=50., error_assimilate=100.,
             cov_loc_radius_km=32, cov_loc_vert_km=5)
 
 
-exp.observations = [ ] #ir108, vis, wv73]  # ir108, wv73, vis]
+exp.observations = [radar, ] # 108, wv73, vis]
 
 # directory paths depend on the name of the experiment
 cluster.expname = exp.expname

scripts/assim_synth_obs.py

@@ -53,7 +53,8 @@ def read_prior_obs(f_obs_prior):
             OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens)))
     return OBSs
 
-def read_obsseqout(f):
+def read_truth_obs_obsseq(f):
+    """Reads observed and true values from obs_seq.out/final files."""
     obsseq = open(f, 'r').readlines()
     true = []
     obs = []
@@ -66,6 +67,44 @@ def read_obsseqout(f):
             obs.append(observed)
     return true, obs
 
+def replace_errors_obsseqout(f, new_errors):
+    """Replaces existing observation errors in obs_seq.final files
+    
+    new_errors (np.array) : standard deviation,
+                            shape must match the number of observations
+    """
+    obsseq = open(f, 'r').readlines()
+
+    # find number of lines between two ' OBS   ' lines:
+    first_obs = None
+    for i, line in enumerate(obsseq):
+        if ' OBS ' in line:
+            if first_obs is None:
+                first_obs = i
+            else:
+                second_obs = i
+                break
+    lines_between = second_obs - first_obs
+    lines_obserr_after_obsnr = lines_between - 1  # obserr line is one before ' OBS   ' line
+
+    # replace values in list obsseq
+    i_obs = 0
+    for i, line in enumerate(obsseq):
+        if ' OBS ' in line:
+            line_error_obs_i = i+lines_obserr_after_obsnr
+
+            previous_err_var = obsseq[line_error_obs_i]
+            new_err_obs_i = new_errors[i_obs]**2  # variance in obs_seq.out
+            #print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
+            obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n'
+
+            i_obs += 1  # next iteration
+
+    with open(f, 'w') as file:
+        for line in obsseq:
+            file.write(line)
+        print(f, 'saved.')
+
 
 def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
                  just_prior_values=False):
@@ -129,7 +168,7 @@ def obs_operator_ensemble(istage):
             os.system('mpirun -np 12 ./perfect_model_obs > /dev/null')
 
             # truth values in obs_seq.out are H(x) values
-            true, _ = read_obsseqout(cluster.dartrundir+'/obs_seq.out')
+            true, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out')
             list_ensemble_truths.append(true)
         
         n_obs = len(list_ensemble_truths[0])
@@ -144,7 +183,7 @@ def obs_operator_nature(time):
     print('running obs operator on nature run')
     prepare_nature_dart(time)
     run_perfect_model_obs()
-    true, _ = read_obsseqout(cluster.dartrundir+'/obs_seq.out')
+    true, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out')
     return true
 
 
@@ -286,32 +325,47 @@ if __name__ == "__main__":
 
         archive_stage = archive_time + '/assim_stage'+str(istage)+'/'
         n_obs = obscfg['n_obs']
+        n_obs_3d = n_obs * len(obscfg['heights'])
         sat_channel = obscfg.get('sat_channel', False)
+        error_generate = obscfg['error_generate']
+        error_assimilate = obscfg['error_assimilate']
 
         set_DART_nml(sat_channel=sat_channel, 
                      cov_loc_radius_km=obscfg['cov_loc_radius_km'],
                      cov_loc_vert_km=obscfg.get('cov_loc_vert_km', False))
 
-        use_error_parametrization = obscfg['err_std'] == False
+        use_error_parametrization = error_generate == False
         if use_error_parametrization:
+
             if sat_channel != 6:
                 raise NotImplementedError('sat channel '+str(sat_channel))
 
-            osq.create_obsseq_in(time, obscfg, zero_error=True)  # zero error to get truth vals
+            osq.create_obsseq_in(time, obscfg, obs_errors=0)  # zero error to get truth vals
 
             Hx_nat = obs_operator_nature(time) 
             Hx_prior = obs_operator_ensemble(istage)  # files are already linked to DART directory
-            
-            obscfg['err_std'] = calc_obserr_WV73(Hx_nat, Hx_prior)
-
+            error_generate = calc_obserr_WV73(Hx_nat, Hx_prior)
+ 
         # create obs template file, now with correct errors
-        osq.create_obsseq_in(time, obscfg, archive_obs_coords=archive_stage+'/obs_coords.pkl')
+        osq.create_obsseq_in(time, obscfg, obs_errors=error_generate,
+                             archive_obs_coords=archive_stage+'/obs_coords.pkl')
         prepare_nature_dart(time)  # link WRF files to DART directory
         run_perfect_model_obs()  # actually create observations that are used to assimilate
 
+        # debug option: overwrite time in prior files
         # for iens in range(1,41):
         #    os.system('ncks -A -v Times '+cluster.dartrundir+'/wrfout_d01 '+cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01')
 
+        if error_assimilate is not False:
+            # overwrite error values in obs_seq.out
+
+            if isinstance(error_assimilate, int):
+                error_assimilate = np.zeros(n_obs_3d) + error_assimilate
+            else:
+                assert len(error_assimilate) == n_obs_3d
+
+            replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
+
         assimilate()
         dir_obsseq = cluster.archivedir()+'/obs_seq_final/assim_stage'+str(istage)
         archive_diagnostics(dir_obsseq, time)

scripts/create_obsseq.py

@@ -211,13 +211,15 @@ kind
 """+str(obs['obserr_var'])
 
 
-def create_obsseq_in(time_dt, obscfg, zero_error=False,
+def create_obsseq_in(time_dt, obscfg, obs_errors=False,
                      archive_obs_coords=False):
     """Create obs_seq.in
 
     Args:
         time_dt (dt.datetime): time of observation
         obscfg (dict)
+        obs_errors (int, np.array) : values of observation errors (standard deviations)
+            e.g. 0 = use zero error
         archive_obs_coords (str, False): path to folder
 
     channel_id (int): SEVIRI channel number
@@ -255,8 +257,8 @@ def create_obsseq_in(time_dt, obscfg, zero_error=False,
 
     # define obs error
     obserr_std = np.zeros(n_obs_3d) 
-    if not zero_error:
-        obserr_std += obscfg['err_std']
+    if obs_errors:
+        obserr_std += obs_errors
 
     # other stuff for obsseq.in
     obs_kind_nr = obs_kind_nrs[kind]
@@ -305,25 +307,44 @@ def create_obsseq_in(time_dt, obscfg, zero_error=False,
 if __name__ == '__main__':
     # for testing
     time_dt = dt.datetime(2008, 7, 30, 9, 0)
-    n_obs = 16  # radar: n_obs for each observation height level
-
-    vis = dict(kind='MSG_4_SEVIRI_BDRF',
-            sat_channel=1, n_obs=n_obs, err_std=0.03,
-            cov_loc_radius_km=10)
-    wv = dict(kind='MSG_4_SEVIRI_TB',
-            sat_channel=6, n_obs=n_obs, err_std=False,
-            cov_loc_radius_km=10)
-    ir108 = dict(kind='MSG_4_SEVIRI_TB',
-                sat_channel=9, n_obs=n_obs, err_std=5.,
-                cov_loc_radius_km=10)
-
-    radar = dict(kind='RADAR_REFLECTIVITY', n_obs=n_obs, err_std=5.,
-                heights=np.arange(1000, 2001, 1000),
-                cov_loc_radius_km=10, cov_loc_vert_km=2)
-
-    t2m = dict(kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=1.0, 
-            cov_loc_radius_km=32, cov_loc_vert_km=1)
-    psfc = dict(kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
-            cov_loc_radius_km=32, cov_loc_vert_km=5)
-
-    create_obsseq_in(time_dt, radar, archive_obs_coords='./coords_stage1.pkl')
+    n_obs = 1600  # radar: n_obs for each observation height level
+
+    vis = dict(plotname='VIS 0.6µm',  plotunits='[1]',
+            kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs, 
+            err_std=0.03,
+            error_generate=0.03, error_assimilate=0.06,
+            cov_loc_radius_km=32)
+
+    wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
+                kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs, 
+                err_std=False,
+                error_generate=False, error_assimilate=False,
+                cov_loc_radius_km=32)
+
+    ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
+                kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs, 
+                err_std=5.,
+                error_generate=5., error_assimilate=10.,
+                cov_loc_radius_km=32)
+
+    radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
+                kind='RADAR_REFLECTIVITY', n_obs=n_obs, 
+                error_generate=2.5, error_assimilate=5.,
+                heights=np.arange(1000, 15001, 1000),
+                cov_loc_radius_km=30, cov_loc_vert_km=4)
+
+    t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
+            kind='SYNOP_TEMPERATURE', n_obs=n_obs, 
+            error_generate=0.1, error_assimilate=1.,
+            cov_loc_radius_km=20, cov_loc_vert_km=3)
+
+    psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
+                kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, 
+                error_generate=50., error_assimilate=100.,
+                cov_loc_radius_km=32, cov_loc_vert_km=5)
+
+    #create_obsseq_in(time_dt, radar, archive_obs_coords=False) #'./coords_stage1.pkl')
+
+    error_assimilate = 5.*np.ones(n_obs*len(radar['heights']))
+    import assim_synth_obs as aso
+    aso.replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)