diff --git a/scripts/gen_synth_obs.py b/scripts/gen_synth_obs.py
index 7a7b9ee3428aab9254a77c82fd5d4e4cdf053a7e..3cd791d60d8506490b8f9f3c33ecf0fed26a6dc7 100755
--- a/scripts/gen_synth_obs.py
+++ b/scripts/gen_synth_obs.py
@@ -91,81 +91,81 @@ 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')
- print(dt.datetime.now())
- # error_var = (5.)**2
- # kind = 'RADAR'
+ # remove any existing observation files
+ os.system('rm -f 'cluster.dartrundir+'/obs_seq_*.out')
- # error_var = (0.5)**2
- # kind = 'RASO_T'
+ for i_obs, obscfg in enumerate(exp.observations):
- #obs.generic_obs(kind, time, exp.n_obs, exp.radar_err**2, exp.distance_between_obs_meters,
- # output_path=cluster.dartrundir,
- # fpath_obs_locations=cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M')
- # +'/obs_coords.pkl')
+ n_obs = obscfg['n_obs']
+ error_var = (obscfg['err_std'])**2
+ distance_between_obs_meters = obscfg['distance_between_obs_meters']
-
- for channel_id in exp.sat_channels:
-
- n_obs = 100
- channel_id = 6
- distance_between_obs_meters = 10000
-
- error_var = (1.)**2
- create_obsseq.sat(time, channel_id, n_obs, error_var,
- distance_between_obs_meters,
- output_path=cluster.dartrundir,
- fpath_obs_locations='./domain.pkl')
- print(dt.datetime.now())
+ # 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')
+ 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 not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
- # generate observations
+ # generate observations (obs_seq.out)
set_input_nml(sat_channel=channel_id, just_prior_values=True)
os.chdir(cluster.dartrundir)
t = dt.datetime.now()
os.system('mpirun -np 12 ./perfect_model_obs')
print('1st perfect_model_obs', (dt.datetime.now()-t).total_seconds())
- if channel_id == 6:
- # run ./filter to have prior observation estimates from model state
- set_input_nml(sat_channel=channel_id, 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())
-
- # find the observation error for a pair of (H(x_nature), H(x_background))
- f_obs_prior = cluster.dartrundir+'/obs_seq.final'
- OBSs = read_prior_obs(f_obs_prior)
-
- # compute the observation error necessary
- # to achieve a certain operational FGD distribution
- OEs = []
- for obs in OBSs:
- bt_y = obs['truth']
- bt_x_ens = obs['prior_ens']
- CIs = [cloudimpact_73(bt_x, bt_y) for bt_x in bt_x_ens]
-
- oe_nature = oe_73(np.mean(CIs))
- OEs.append(oe_nature)
-
- # write obs_seq.out
- fpath_obsseqout = cluster.dartrundir+'/obs_seq.in'
- edit_obserr_in_obsseq(fpath_obsseqout, OEs)
- print('after editing oe', dt.datetime.now())
-
- # generate actual observations (with correct error)
- os.chdir(cluster.dartrundir)
- t = dt.datetime.now()
- os.system('mpirun -np 12 ./perfect_model_obs')
- print('real obs gen', (dt.datetime.now()-t).total_seconds())
-
- # correct input.nml for actual assimilation later on
- set_input_nml(sat_channel=channel_id, just_prior_values=False)
-
- # FIXME: missing is the concatenation of
+ # cloud dependent observation error
+ if obscfg['sat']:
+ if obscfg['channel'] == 6:
+ # run ./filter to have prior observation estimates from model state
+ set_input_nml(sat_channel=obscfg['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())
+
+ # find the observation error for a pair of (H(x_nature), H(x_background))
+ f_obs_prior = cluster.dartrundir+'/obs_seq.final'
+ OBSs = read_prior_obs(f_obs_prior)
+
+ # compute the observation error necessary
+ # to achieve a certain operational FGD distribution
+ OEs = []
+ for obs in OBSs:
+ bt_y = obs['truth']
+ bt_x_ens = obs['prior_ens']
+ CIs = [cloudimpact_73(bt_x, bt_y) for bt_x in bt_x_ens]
+
+ oe_nature = oe_73(np.mean(CIs))
+ OEs.append(oe_nature)
+
+ # write obs_seq.out
+ fpath_obsseqout = cluster.dartrundir+'/obs_seq.in'
+ edit_obserr_in_obsseq(fpath_obsseqout, OEs)
+
+ # generate actual observations (with correct error)
+ os.chdir(cluster.dartrundir)
+ t = dt.datetime.now()
+ os.system('mpirun -np 12 ./perfect_model_obs')
+ 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)
+
+ # rename according to i_obs
+ 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)
- os.system('cat obs_seq.out >> obs_seq_all.out')
- print(dt.datetime.now())
- # FIXME: what if different channels in one obs_seq.out -> need different input.nml for different channels
+ os.system('cat obs_seq_*.out >> obs_seq_all.out')
+
+ print(dt.datetime.now())
\ No newline at end of file