restructure obs_loc

config/cfg.py

@@ -12,12 +12,15 @@ exp.model_dx = 2000
 exp.timestep = 10
 exp.n_ens = 40
 exp.n_nodes = 10
-exp.n_obs = 81  # radar: n_obs for each observation height level
+
+n_obs = 81  # radar: n_obs for each observation height level
 
 vis = dict(sat=True, channel=1, n_obs=n_obs, err_std=0.03, 
-           distance_between_obs_meters=50000,)
+           distance_between_obs_km=50, 
+           cov_loc_radius_km=10)
 radar = dict(sat=False, kind='radar', n_obs=n_obs, err_std=5.,
-             distance_between_obs_meters=50000,)
+             distance_between_obs_km=50, 
+             cov_loc_radius_km=10)
 
 exp.observations = [vis,]
 

scripts/create_obsseq.py

@@ -33,13 +33,18 @@ def get_dart_date(time_dt):
     secs_thatday = str(int((time_dt - round_to_day(time_dt)).total_seconds()))
     return dart_date_day, secs_thatday
 
-def calc_obs_locations(n_obs, distance_between_obs_meters,
-                       coords_from_domaincenter=True, fpath_obs_locations=False):
+def calc_obs_locations(n_obs, coords_from_domaincenter=True,
+                       distance_between_obs_km=9999, fpath_obs_locations=False):
     """Calculate coordinate pairs for locations
 
     Args:
         n_obs (int):
             number of observations (must be a square of an integer: 4, 9, 1000, ...)
+        coords_from_domaincenter (bool):
+            if False: spread observations evenly
+            if True: spread from domaincenter, `distance_between_obs_km` apart
+        distance_between_obs_km (int):
+            only used when coords_from_domaincenter=True
         fpath_obs_locations (False or str):
             write an obs_coords.pkl, can be used to check observation locations
             if str, write to file
@@ -60,6 +65,7 @@ def calc_obs_locations(n_obs, distance_between_obs_meters,
         nx, ny = int(np.sqrt(n_obs)), int(np.sqrt(n_obs))
 
         m_per_degree = 2*np.pi*radius_earth_meters/360
+        distance_between_obs_meters = distance_between_obs_km*1000.
         dy_4km_in_degree = distance_between_obs_meters/m_per_degree
 
         for iy in range(int(-ny/2), int(ny/2+1)):
@@ -157,27 +163,22 @@ kind
         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):
+def sat(time_dt, channel_id, coords, error_var, output_path='./'):
     """Create obs_seq.in
 
     Args:
         time_dt (dt.datetime): time of observation
         channel_id (int): SEVIRI channel number
             see https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov12/ir_srf/rtcoef_msg_4_seviri_srf.html
-        n_obs (int):
-            number of observations squared (must be a square of an integer: 4, 9, 1000, ...)
+        coords (list of 2-tuples with (lat,lon))
         error_var (float):
             gaussian error with this variance is added to the truth at observation time
         output_path (str): directory where `obs_seq.in` will be saved
-        fpath_obs_locations (False or str):
-            write an obs_coords.pkl, can be used to check observation locations
-            if str, write to file
     """
     # time_dt = add_timezone_UTC(time_dt)
     # time_dt = dt.datetime(2008, 7, 30, 15, 30, tzinfo=dt.timezone.utc)
-    assert n_obs == int(n_obs)
     error_var = str(error_var)
+    n_obs = len(coords)
 
     # Brightness temperature or Reflectance?
     channel_id = int(channel_id)
@@ -189,10 +190,6 @@ def sat(time_dt, channel_id, n_obs, error_var, distance_between_obs_meters,
         code = '255'
     channel_id = str(channel_id)
 
-    coords = calc_obs_locations(n_obs, distance_between_obs_meters,
-                                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))
     sun_zen = str(90. - get_altitude(lat0, lon0, time_dt))
@@ -269,13 +266,7 @@ kind
         print(fpath, 'saved.')
 
 
-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=coords_from_domaincenter,
-                fpath_obs_locations=fpath_obs_locations)
-
+def calc_obs_locations_3d(coords, heights):
     # append height
     coords2 = []
     for i in range(len(coords)):
@@ -288,18 +279,14 @@ def gen_coords(n_obs, distance_between_obs_meters, heights,
     return coords2
 
 
-def generic_obs(obs_type, time_dt, n_obs, error_var, distance_between_obs_meters,
-                output_path='./', fpath_obs_locations=False):
+def generic_obs(obs_type, time_dt, coords, error_var, output_path='./'):
 
     obs_codes = {'RASO_T': {'name': 'RADIOSONDE_TEMPERATURE', 'nr': '5'},
                  'RADAR': {'name': 'RADAR_REFLECTIVITY', 'nr': '37'}
                  }
 
     heights = np.arange(1000, 15001, 1000)
-
-    coords = gen_coords(n_obs, distance_between_obs_meters, heights,
-                        coords_from_domaincenter=False, 
-                        fpath_obs_locations=fpath_obs_locations)
+    coords = calc_obs_locations_3d(coords, heights)
 
     dart_date_day, secs_thatday = get_dart_date(time_dt)
     print('secs, days:', secs_thatday, dart_date_day)

scripts/gen_synth_obs.py

@@ -4,7 +4,9 @@ import numpy as np
 from scipy.interpolate import interp1d
 from config.cfg import exp, cluster
 from utils import symlink, copy, sed_inplace, append_file
-import create_obsseq
+import create_obsseq as osq
+
+earth_radius_km = 6370
 
 # fit of Fig 7, Harnisch 2016
 x_ci = [0,   5, 10.5, 13, 16]
@@ -70,7 +72,8 @@ def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
         for line in obsseq_new:
             f.write(line)
 
-def set_input_nml(sat_channel=False, just_prior_values=False):
+def set_input_nml(sat_channel=False, just_prior_values=False,
+                  cov_loc_radius_km=32):
     """descr"""
     if just_prior_values:
         template = cluster.scriptsdir+'/../templates/input.prioronly.nml'
@@ -78,6 +81,8 @@ def set_input_nml(sat_channel=False, just_prior_values=False):
         template = cluster.scriptsdir+'/../templates/input.nml'
     copy(template, cluster.dartrundir+'/input.nml')
     sed_inplace(cluster.dartrundir+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
+    cov_loc_radian = cov_loc_radius_km/earth_radius_km
+    sed_inplace(cluster.dartrundir+'/input.nml', '<cov_loc_radian>', str(cov_loc_radian))
 
     # input.nml for RTTOV
     if sat_channel > 0:
@@ -91,6 +96,7 @@ def set_input_nml(sat_channel=False, just_prior_values=False):
 if __name__ == "__main__":
 
     time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
+    fpath_obs_coords = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/obs_coords.pkl')
 
     # remove any existing observation files
     os.system('rm -f '+cluster.dartrundir+'/obs_seq_*.out')
@@ -100,32 +106,27 @@ if __name__ == "__main__":
 
         n_obs = obscfg['n_obs']
         error_var = (obscfg['err_std'])**2
-        distance_between_obs_meters = obscfg['distance_between_obs_meters']
+        sat_channel = obscfg.get('channel', False)
+        cov_loc = obscfg.get['cov_loc_radius_km']
+        dist_obs = obscfg.get('distance_between_obs_km', False)
 
         # generate obs_seq.in
-        if not obscfg['sat']:
-            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')
-            set_input_nml()
-        else:
-            create_obsseq.sat(time, obscfg['channel'], n_obs, error_var,
-                                distance_between_obs_meters,
-                                output_path=cluster.dartrundir, 
-                                fpath_obs_locations='./domain.pkl')
-            
-            if obscfg['channel'] == 6:
-                # for cloud dependent error
-                set_input_nml(sat_channel=obscfg['channel'], just_prior_values=True)
+        obs_coords = osq.calc_obs_locations(n_obs, coords_from_domaincenter=False, 
+                                            distance_between_obs_km=dist_obs, 
+                                            fpath_obs_locations=fpath_obs_coords)
+
+        if obscfg['sat']:
+            osq.sat(time, sat_channel, obs_coords, error_var,
+                    output_path=cluster.dartrundir)
         else:
-                set_input_nml(sat_channel=obscfg['channel'])
+            osq.generic_obs(obscfg['kind'], time, obs_coords, error_var,
+                            output_path=cluster.dartrundir)
 
         if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
             raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
 
         # generate observations (obs_seq.out)
+        set_input_nml(sat_channel=sat_channel, cov_loc_radius_km=cov_loc)
         os.chdir(cluster.dartrundir)
         t = dt.datetime.now()
         os.system('mpirun -np 12 ./perfect_model_obs')
@@ -133,9 +134,9 @@ if __name__ == "__main__":
 
         # cloud dependent observation error
         if obscfg['sat']:
-            if obscfg['channel'] == 6:
+            if sat_channel == 6:
                 # run ./filter to have prior observation estimates from model state
-                set_input_nml(sat_channel=obscfg['channel'], just_prior_values=True)
+                set_input_nml(sat_channel=sat_channel, just_prior_values=True)
                 t = dt.datetime.now()
                 os.system('mv obs_seq.out obs_seq_all.out; mpirun -np 20 ./filter')
                 print('1st filter', (dt.datetime.now()-t).total_seconds())
@@ -166,10 +167,12 @@ if __name__ == "__main__":
                 print('real obs gen', (dt.datetime.now()-t).total_seconds())
 
                 # correct input.nml for actual assimilation later on
-                set_input_nml(sat_channel=obscfg['channel'], just_prior_values=False)
+                set_input_nml(sat_channel=sat_channel,
+                              cov_loc_radius_km=cov_loc)
 
         # rename according to i_obs
-        os.rename(cluster.dartrundir+'/obs_seq.out', cluster.dartrundir+'/obs_seq_'+str(i_obs)+'.out')
+        os.rename(cluster.dartrundir+'/obs_seq.out', 
+                  cluster.dartrundir+'/obs_seq_'+str(i_obs)+'.out')
 
     # concatenate the created obs_seq_*.out files
     os.chdir(cluster.dartrundir)