minor

scripts/create_obsseq.py

@@ -27,12 +27,8 @@ def add_timezone_UTC(t):
 
 def get_dart_date(time_dt):
     # assumes input is UTC!
-    #try:
-    #    time_dt = add_timezone_UTC(time_dt)
-    #except:
-    #    time_dt = time_dt.replace(tzinfo=dt.timezone.utc)  # overwrites existing timezone!
     days_since_010120 = 145731
-    ref_day = dt.datetime(2000,1,1, tzinfo=dt.timezone.utc)
+    ref_day = dt.datetime(2000, 1, 1, tzinfo=dt.timezone.utc)
     dart_date_day = str((time_dt - ref_day).days + days_since_010120)
     secs_thatday = str(int((time_dt - round_to_day(time_dt)).total_seconds()))
     return dart_date_day, secs_thatday
@@ -102,6 +98,65 @@ def calc_obs_locations(n_obs, distance_between_obs_meters,
     assert len(coords) == n_obs, (len(coords), n_obs)
     return coords
 
+def write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords,
+                         dart_date_day, secs_thatday, output_path):
+    n_obs_str = str(int(n_obs))
+    error_var = str(error_var)
+    line_obstypedef = obs_kind_nr+' '+obs_name
+
+    msg = """
+ obs_sequence
+obs_kind_definitions
+           1
+         """+line_obstypedef+"""
+  num_copies:            0  num_qc:            0
+  num_obs:            """+n_obs_str+"  max_num_obs:            "+n_obs_str+"""
+  first:            1  last:            """+n_obs_str
+
+    for i_obs in range(1, int(n_obs)+1):
+
+        lon = coords[i_obs-1][1]
+        lat = coords[i_obs-1][0]
+        hgt_m = str(coords[i_obs-1][2])+'.000'
+
+        lon_rad = str(degr_to_rad(lon))
+        lat_rad = str(degr_to_rad(lat))
+
+        # compile text
+        if i_obs < int(n_obs):
+            msg += """
+ OBS            """+str(i_obs)+"""
+          -1           """+str(i_obs+1)+"""          -1
+obdef
+loc3d
+     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
+kind
+         """+obs_kind_nr+"""
+ """+secs_thatday+"""     """+dart_date_day+"""
+  """+error_var
+        if i_obs == int(n_obs):  # last_observation
+            # compile text
+            msg += """
+ OBS            """+str(i_obs)+"""
+          """+str(i_obs-1)+"""           -1          -1
+obdef
+loc3d
+     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
+kind
+         """+obs_kind_nr+"""
+ """+secs_thatday+"""     """+dart_date_day+"""
+  """+error_var
+
+    fpath = output_path+'/obs_seq.in'
+    try:
+        os.remove(fpath)
+    except OSError:
+        pass
+
+    with open(fpath, 'w') as f:
+        f.write(msg)
+        print(fpath, 'saved.')
+
 def sat(time_dt, channel_id, n_obs, error_var, distance_between_obs_meters,
         output_path='./', fpath_obs_locations=False):
     """Create obs_seq.in
@@ -135,8 +190,8 @@ def sat(time_dt, channel_id, n_obs, error_var, distance_between_obs_meters,
     channel_id = str(channel_id)
 
     coords = calc_obs_locations(n_obs, distance_between_obs_meters,
-                coords_from_domaincenter=False,
-                fpath_obs_locations=False)
+                                coords_from_domaincenter=False,
+                                fpath_obs_locations=fpath_obs_locations)
 
     time_dt = add_timezone_UTC(time_dt)
     sun_az = str(get_azimuth(lat0, lon0, time_dt))
@@ -218,7 +273,8 @@ def gen_coords(n_obs, distance_between_obs_meters, heights,
                coords_from_domaincenter=True, fpath_obs_locations=False):
 
     coords = calc_obs_locations(n_obs, distance_between_obs_meters,
-                coords_from_domaincenter=True, fpath_obs_locations=False)
+                coords_from_domaincenter=coords_from_domaincenter,
+                fpath_obs_locations=fpath_obs_locations)
 
     # append height
     coords2 = []
@@ -227,7 +283,8 @@ def gen_coords(n_obs, distance_between_obs_meters, heights,
             coords2.append(coords[i] + (hgt_m,))
 
     n_obs = len(coords2)
-    print('effective number of observations (with vertical levels):', n_obs, ' on each level:', len(coords))
+    print('effective number of observations (with vertical levels):', n_obs, 
+          ' on each level:', len(coords))
     return coords2
 
 
@@ -236,13 +293,13 @@ def generic_obs(obs_type, time_dt, n_obs, error_var, distance_between_obs_meters
 
     obs_codes = {'RASO_T': {'name': 'RADIOSONDE_TEMPERATURE', 'nr': '5'},
                  'RADAR': {'name': 'RADAR_REFLECTIVITY', 'nr': '37'}
-                }
+                 }
 
-    coords_from_domaincenter = True
     heights = np.arange(1000, 15001, 1000)
 
     coords = gen_coords(n_obs, distance_between_obs_meters, heights,
-                        coords_from_domaincenter=True, fpath_obs_locations=False)
+                        coords_from_domaincenter=False, 
+                        fpath_obs_locations=fpath_obs_locations)
 
     dart_date_day, secs_thatday = get_dart_date(time_dt)
     print('secs, days:', secs_thatday, dart_date_day)
@@ -251,68 +308,9 @@ def generic_obs(obs_type, time_dt, n_obs, error_var, distance_between_obs_meters
     obs_kind_nr = obs_codes[obs_type]['nr']
 
     write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords,
-                            dart_date_day, secs_thatday, output_path)
+                         dart_date_day, secs_thatday, output_path)
 
 
-def write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords,
-                         dart_date_day, secs_thatday, output_path):
-    n_obs_str = str(int(n_obs))
-    error_var = str(error_var)
-    line_obstypedef = obs_kind_nr+' '+obs_name
-
-    msg = """
- obs_sequence
-obs_kind_definitions
-           1
-         """+line_obstypedef+"""
-  num_copies:            0  num_qc:            0
-  num_obs:            """+n_obs_str+"  max_num_obs:            "+n_obs_str+"""
-  first:            1  last:            """+n_obs_str
-
-    for i_obs in range(1, int(n_obs)+1):
-
-        lon = coords[i_obs-1][1]
-        lat = coords[i_obs-1][0]
-        hgt_m = str(coords[i_obs-1][2])+'.000'
-
-        lon_rad = str(degr_to_rad(lon))
-        lat_rad = str(degr_to_rad(lat))
-
-        # compile text
-        if i_obs < int(n_obs):
-            msg += """
- OBS            """+str(i_obs)+"""
-          -1           """+str(i_obs+1)+"""          -1
-obdef
-loc3d
-     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
-kind
-         """+obs_kind_nr+"""
- """+secs_thatday+"""     """+dart_date_day+"""
-  """+error_var
-        if i_obs == int(n_obs):  # last_observation
-            # compile text
-            msg += """
- OBS            """+str(i_obs)+"""
-          """+str(i_obs-1)+"""           -1          -1
-obdef
-loc3d
-     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
-kind
-         """+obs_kind_nr+"""
- """+secs_thatday+"""     """+dart_date_day+"""
-  """+error_var
-
-    fpath = output_path+'/obs_seq.in'
-    try:
-        os.remove(fpath)
-    except OSError:
-        pass
-
-    with open(fpath, 'w') as f:
-        f.write(msg)
-        print(fpath, 'saved.')
-
 if __name__ == '__main__':
     time_dt = dt.datetime(2008, 7, 30, 15, 30)
     n_obs = 100
@@ -326,7 +324,7 @@ if __name__ == '__main__':
 
     error_var = (5.)**2
     generic_obs('RADAR', time_dt, n_obs, error_var, distance_between_obs_meters,
-                    output_path='./', fpath_obs_locations='./domain.pkl')
+                output_path='./', fpath_obs_locations='./domain.pkl')
 
     # error_var = (0.5)**2
     # generic_obs('RASO_T', time_dt, n_obs, error_var, distance_between_obs_meters,

scripts/gen_synth_obs.py

@@ -46,8 +46,7 @@ def read_prior_obs(f_obs_prior):
             for j in range(5, 5+exp.n_ens):
                 prior_ens.append(float(obsseq[i+j].strip()))
 
-            OBSs.append(dict(observed=observed, truth=truth,
-                            prior_ens=np.array(prior_ens)))
+            OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens)))
     return OBSs
 
 def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
@@ -88,6 +87,7 @@ def set_input_nml(sat_channel=False, just_prior_values=False):
             rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
         append_file(cluster.dartrundir+'/input.nml', rttov_nml)
 
+
 if __name__ == "__main__":
 
     time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
@@ -103,16 +103,16 @@ if __name__ == "__main__":
 
         # generate obs_seq.in
         if not obscfg['sat']:
-            create_obsseq.generic_obs(obscfg['kind'], time, n_obs, error_var, 
+            create_obsseq.generic_obs(obscfg['kind'], time, n_obs, error_var,
                         distance_between_obs_meters,
                         output_path=cluster.dartrundir,
                         fpath_obs_locations=cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M')
                                             +'/obs_coords.pkl')
         else:
-            create_obsseq.sat(time, obscfg['channel'], n_obs, error_var, 
-                            distance_between_obs_meters,
-                            output_path=cluster.dartrundir, 
-                            fpath_obs_locations='./domain.pkl')
+            create_obsseq.sat(time, obscfg['channel'], n_obs, error_var,
+                                distance_between_obs_meters,
+                                output_path=cluster.dartrundir, 
+                                fpath_obs_locations='./domain.pkl')
 
         if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
             raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)

scripts/run_obs_diag.py

@@ -17,26 +17,33 @@ def run(folder_obs_seq_final):
             f.write(fin)
             f.write('\n')
 
+    for obserr_iszero in ['.true.', '.false.']:
+        print('ensure correct input.nml')
+        copy(cluster.scriptsdir+'/../templates/input.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')
+
+        # run obs_diag
+        print('running obs_diag program')
+        os.chdir(rundir_program)
+        symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
+        try:
+            os.remove(rundir_program+'/obs_seq_to_netcdf')
+        except:
+            pass
+        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])
+        if obserr_iszero == '.true.':
+            fout = '/obs_diag_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)
 
-    print('ensure correct input.nml')
-    copy(cluster.scriptsdir+'/../templates/input.nml',
-         rundir_program+'/input.nml') #cluster.dartrundir+'/input.nml')
-    sed_inplace(rundir_program+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
-    append_file(rundir_program+'/input.nml', cluster.scriptsdir+'/../templates/obs_def_rttov.VIS.nml')
-
-    # run obs_diag
-    print('running obs_diag program')
-    os.chdir(rundir_program)
-    symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
-    try:
-        os.remove(rundir_program+'/obs_seq_to_netcdf')
-    except:
-        pass
-    os.system('./obs_diag >& obs_diag.log')  # caution, this overwrites obs_seq_to_netcdf
-
-    outdir = '/'.join(folder_obs_seq_final.split('/')[:-1])
-    print('moving output to', outdir+'/obs_diag_output.nc')
-    copy(rundir_program+'/obs_diag_output.nc', outdir+'/obs_diag_output.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')

templates/input.nml

@@ -316,7 +316,7 @@
    print_mismatched_locs = .false.,
    create_rank_histogram = .true.,
    outliers_in_histogram = .true.,
-   use_zero_error_obs    = .true.,
+   use_zero_error_obs    = <zero_error_obs>,
    verbose               = .false.
    /