diff --git a/config/cfg.py b/config/cfg.py
index c3f91b77bb11f5626796ee5dfe0c49dc27c78007..91c91dc7c34cdbb34b96e490e729274c30e6aea6 100755
--- a/config/cfg.py
+++ b/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,]
diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py
index cb7a161a3a82be81d6caaa6e7724979d8e3862ba..91b39dff9328115ca2cf8d2e2141eafe064158e8 100755
--- a/scripts/create_obsseq.py
+++ b/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)):
@@ -84,8 +90,8 @@ def calc_obs_locations(n_obs, distance_between_obs_meters,
for i in range(n_obs_x):
for j in range(n_obs_x):
- coords.append((lats[skip+i*dx,skip+j*dx],
- lons[skip+i*dx,skip+j*dx]))
+ coords.append((lats[skip+i*dx, skip+j*dx],
+ lons[skip+i*dx, skip+j*dx]))
try:
if fpath_obs_locations:
@@ -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)):
@@ -283,23 +274,19 @@ 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,
+ print('effective number of observations (with vertical levels):', n_obs,
' on each level:', len(coords))
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)
diff --git a/scripts/gen_synth_obs.py b/scripts/gen_synth_obs.py
index 25041296c12f096e15f4a22646fe7f37509a009e..c655b85558e3bc9fe3cff1c5d0298d9c10dd89b8 100755
--- a/scripts/gen_synth_obs.py
+++ b/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()
+ 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:
- 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)
- 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)