diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py
index 58a79384e1e6bf59a8ef76123df43072cabf77ec..cb7a161a3a82be81d6caaa6e7724979d8e3862ba 100755
--- a/scripts/create_obsseq.py
+++ b/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,
diff --git a/scripts/gen_synth_obs.py b/scripts/gen_synth_obs.py
index 3cd791d60d8506490b8f9f3c33ecf0fed26a6dc7..46a4ab3b45743485b860a2f497db3aad13feb30e 100755
--- a/scripts/gen_synth_obs.py
+++ b/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)
diff --git a/scripts/run_obs_diag.py b/scripts/run_obs_diag.py
index 5fd07f90752f0a62f16e59eb23bbb25832663b43..41a8c260174425e80b1c727871ce71b1c66ae077 100644
--- a/scripts/run_obs_diag.py
+++ b/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')
diff --git a/templates/input.nml b/templates/input.nml
index 101c3fdc534fd8ef5e8e2a9a8762cfe5e8193174..8787dd0e75cc1de87f2b881775697dc48178c901 100644
--- a/templates/input.nml
+++ b/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.
/