Skip to content
Snippets Groups Projects
Commit 8ed445f0 authored by Lukas Kugler's avatar Lukas Kugler
Browse files

complete rewrite create obsseq

parent ded54c86
No related branches found
No related tags found
No related merge requests found
...@@ -17,11 +17,14 @@ exp.n_nodes = 10 ...@@ -17,11 +17,14 @@ exp.n_nodes = 10
n_obs = 64 # radar: n_obs for each observation height level n_obs = 64 # radar: n_obs for each observation height level
vis = dict(sat_channel=1, n_obs=n_obs, err_std=0.03, vis = dict(kind='MSG_4_SEVIRI_BDRF',
sat_channel=1, n_obs=n_obs, err_std=0.03,
cov_loc_radius_km=10) cov_loc_radius_km=10)
wv = dict(sat_channel=6, n_obs=n_obs, err_std=False, wv = dict(kind='MSG_4_SEVIRI_TB',
sat_channel=6, n_obs=n_obs, err_std=False,
cov_loc_radius_km=10) cov_loc_radius_km=10)
ir108 = dict(sat_channel=9, n_obs=n_obs, err_std=5., ir108 = dict(kind='MSG_4_SEVIRI_TB',
sat_channel=9, n_obs=n_obs, err_std=5.,
cov_loc_radius_km=10) cov_loc_radius_km=10)
radar = dict(kind='RADAR', n_obs=n_obs, err_std=5., radar = dict(kind='RADAR', n_obs=n_obs, err_std=5.,
......
...@@ -88,17 +88,14 @@ def edit_obserr_in_obsseq(fpath_obsseqin, OEs): ...@@ -88,17 +88,14 @@ def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
f.write(line) f.write(line)
def set_input_nml(sat_channel=False, just_prior_values=False, def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False):
cov_loc_radius_km=32, cov_loc_vert_km=False):
"""descr""" """descr"""
cov_loc_radian = cov_loc_radius_km/earth_radius_km cov_loc_radian = cov_loc_radius_km/earth_radius_km
if just_prior_values: copy(cluster.scriptsdir+'/../templates/input.nml',
template = cluster.scriptsdir+'/../templates/input.prioronly.nml' cluster.dartrundir+'/input.nml')
else:
template = cluster.scriptsdir+'/../templates/input.nml'
copy(template, cluster.dartrundir+'/input.nml')
# options are overwritten with settings
options = {'<n_ens>': str(int(exp.n_ens)), options = {'<n_ens>': str(int(exp.n_ens)),
'<cov_loc_radian>': str(cov_loc_radian)} '<cov_loc_radian>': str(cov_loc_radian)}
...@@ -160,7 +157,6 @@ def obs_operator_ensemble(): ...@@ -160,7 +157,6 @@ def obs_operator_ensemble():
def obs_operator_nature(): def obs_operator_nature():
prepare_nature_dart() prepare_nature_dart()
set_input_nml(sat_channel=sat_channel)
os.chdir(cluster.dartrundir) os.chdir(cluster.dartrundir)
os.remove(cluster.dartrundir+'/obs_seq.out') os.remove(cluster.dartrundir+'/obs_seq.out')
os.system('mpirun -np 12 ./perfect_model_obs') os.system('mpirun -np 12 ./perfect_model_obs')
...@@ -194,11 +190,7 @@ def calc_obserr_WV73(Hx_nature, Hx_prior): ...@@ -194,11 +190,7 @@ def calc_obserr_WV73(Hx_nature, Hx_prior):
OEs[iobs] = oe_nature OEs[iobs] = oe_nature
return OEs return OEs
def generate_observations(): def generate_observations():
set_input_nml(sat_channel=sat_channel, cov_loc_radius_km=cov_loc,
cov_loc_vert_km=cov_loc_vert_km)
# generate actual observations (with correct error) # generate actual observations (with correct error)
os.chdir(cluster.dartrundir) os.chdir(cluster.dartrundir)
os.remove(cluster.dartrundir+'/obs_seq.out') os.remove(cluster.dartrundir+'/obs_seq.out')
...@@ -206,7 +198,6 @@ def generate_observations(): ...@@ -206,7 +198,6 @@ def generate_observations():
raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir) raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
os.system('mpirun -np 12 ./perfect_model_obs') os.system('mpirun -np 12 ./perfect_model_obs')
def assimilate(): def assimilate():
os.chdir(cluster.dartrundir) os.chdir(cluster.dartrundir)
os.remove(cluster.dartrundir+'/obs_seq.final') os.remove(cluster.dartrundir+'/obs_seq.final')
...@@ -214,17 +205,14 @@ def assimilate(): ...@@ -214,17 +205,14 @@ def assimilate():
raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir) raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir)
os.system('mpirun -np 48 ./filter') os.system('mpirun -np 48 ./filter')
def archive_diagnostics(time): def archive_diagnostics(archive_stage):
print('archive obs space diagnostics') print('archive obs space diagnostics')
savedir = cluster.archivedir()+'/obs_seq_final/' mkdir(archive_stage)
mkdir(savedir) copy(cluster.dartrundir+'/obs_seq.final', archive_stage+'/obs_seq.final')
copy(cluster.dartrundir+'/obs_seq.final', savedir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final'))
try: try:
print('archive regression diagnostics') print('archive regression diagnostics')
savedir = cluster.archivedir()+'/reg_factor/' copy(cluster.dartrundir+'/reg_diagnostics', archive_stage+'/reg_diagnostics')
mkdir(savedir)
copy(cluster.dartrundir+'/reg_diagnostics', savedir+time.strftime('/%Y-%m-%d_%H:%M_reg_diagnostics'))
except Exception as e: except Exception as e:
warnings.warn(str(e)) warnings.warn(str(e))
...@@ -234,26 +222,19 @@ def recycle_output(): ...@@ -234,26 +222,19 @@ def recycle_output():
os.rename(cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4), os.rename(cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4),
cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01') cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01')
def archive_output_mean(archive_stage):
def archive_output_mean(time, istage):
istage = str(istage)
print('copy prior posterior members to archive')
for iens in range(1, exp.n_ens+1): for iens in range(1, exp.n_ens+1):
savedir = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/')+str(iens) savedir = archive_stage+'/'+str(iens)
mkdir(savedir) mkdir(savedir)
copy(cluster.dartrundir+'/input.nml', copy(cluster.dartrundir+'/input.nml', archive_stage+'/input.nml')
cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/input '+istage+'.nml'))
# filter_in = cluster.dartrundir+'/preassim_member_'+str(iens).zfill(4)+'.nc' # filter_in = cluster.dartrundir+'/preassim_member_'+str(iens).zfill(4)+'.nc'
# filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4) # filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
# copy mean and sd to archive # copy mean and sd to archive
print('copy preassim, postassim mean and sd')
for f in ['output_mean.nc', 'output_sd.nc']: for f in ['output_mean.nc', 'output_sd.nc']:
copy(cluster.dartrundir+'/'+f, copy(cluster.dartrundir+'/'+f, archive_stage+'/'+f)
cluster.archivedir()+'/'+f[:-3]+time.strftime('_'+istage+'_%Y-%m-%d_%H:%M:%S'))
if __name__ == "__main__": if __name__ == "__main__":
...@@ -279,20 +260,19 @@ if __name__ == "__main__": ...@@ -279,20 +260,19 @@ if __name__ == "__main__":
""" """
time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M') time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
folder_obs_coords = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/') archive_time = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/')
os.chdir(cluster.dartrundir) os.chdir(cluster.dartrundir)
os.system('rm -f obs_seq.out obs_seq.in obs_seq.final') # remove any existing observation files os.system('rm -f obs_seq.in obs_seq.out obs_seq.final') # remove any existing observation files
n_stages = len(exp.observations) n_stages = len(exp.observations)
for istage, obscfg in enumerate(exp.observations): for istage, obscfg in enumerate(exp.observations):
kind = obscfg['kind']
archive_stage = archive_time + '/stage'+str(istage)+'_'+kind
n_obs = obscfg['n_obs'] n_obs = obscfg['n_obs']
sat_channel = obscfg.get('sat_channel', False) sat_channel = obscfg.get('sat_channel', False)
obs_coords = osq.calc_obs_locations(n_obs, obscfg['folder_obs_coords'] = archive_stage+'/obs_coords.pkl'
coords_from_domaincenter=False,
distance_between_obs_km=obscfg.get('distance_between_obs_km', False),
folder_obs_coords=folder_obs_coords)
set_DART_nml(sat_channel=sat_channel, set_DART_nml(sat_channel=sat_channel,
cov_loc=obscfg['cov_loc_radius_km'], cov_loc=obscfg['cov_loc_radius_km'],
...@@ -303,25 +283,25 @@ if __name__ == "__main__": ...@@ -303,25 +283,25 @@ if __name__ == "__main__":
if sat_channel != 6: if sat_channel != 6:
raise NotImplementedError('sat channel '+str(sat_channel)) raise NotImplementedError('sat channel '+str(sat_channel))
osq.create_obsseq_in(obscfg, np.zeros(n_obs)) osq.create_obsseq_in(time, obscfg, zero_error=True) # zero error to get truth vals
Hx_nat = obs_operator_nature() Hx_nat = obs_operator_nature()
Hx_prior = obs_operator_ensemble() # files are already linked to DART directory Hx_prior = obs_operator_ensemble() # files are already linked to DART directory
obserr_var = calc_obserr_WV73(Hx_nat, Hx_prior) obscfg['err_std'] = calc_obserr_WV73(Hx_nat, Hx_prior)
else: else:
obserr_var = np.ones(n_obs) * obscfg['err_std']**2 obscfg['err_std'] = np.ones(n_obs) * obscfg['err_std']
osq.create_obsseq_in(obscfg, obserr_var) # now with correct errors osq.create_obsseq_in(time, obscfg) # now with correct errors
generate_observations() generate_observations()
assimilate() assimilate()
archive_diagnostics(time) archive_diagnostics(archive_stage)
if istage < n_stages-1: if istage < n_stages-1:
# recirculation: filter output -> input # recirculation: filter output -> input
recycle_output() recycle_output()
archive_output_mean(time, istage) archive_output_mean(archive_stage)
elif istage == n_stages-1: elif istage == n_stages-1:
# last assimilation, continue integration now # last assimilation, continue integration now
......
...@@ -4,8 +4,45 @@ ...@@ -4,8 +4,45 @@
import os, sys, warnings import os, sys, warnings
import numpy as np import numpy as np
import datetime as dt import datetime as dt
from config.cfg import exp, cluster
from pysolar.solar import get_altitude, get_azimuth from pysolar.solar import get_altitude, get_azimuth
def obskind_read():
raw_obskind_dart = """
5 RADIOSONDE_TEMPERATURE
6 RADIOSONDE_SPECIFIC_HUMIDITY
12 AIRCRAFT_U_WIND_COMPONENT
13 AIRCRAFT_V_WIND_COMPONENT
14 AIRCRAFT_TEMPERATURE
16 ACARS_U_WIND_COMPONENT
17 ACARS_V_WIND_COMPONENT
18 ACARS_TEMPERATURE
29 LAND_SFC_PRESSURE
30 SAT_U_WIND_COMPONENT
31 SAT_V_WIND_COMPONENT
36 DOPPLER_RADIAL_VELOCITY
37 RADAR_REFLECTIVITY
83 GPSRO_REFRACTIVITY
94 SYNOP_SURFACE_PRESSURE
95 SYNOP_SPECIFIC_HUMIDITY
96 SYNOP_TEMPERATURE
254 MSG_4_SEVIRI_RADIANCE
255 MSG_4_SEVIRI_TB
256 MSG_4_SEVIRI_BDRF"""
# lookup table for kind nr
alist = raw_obskind_dart.split()
assert len(alist) % 2 == 0, alist
obskind_nrs = {}
for i in range(0, len(alist)-1, 2):
obskind_nrs[alist[i+1]] = alist[i]
return obskind_nrs
#####################
# Global variables
# position on earth to calculate solar angles # position on earth to calculate solar angles
lat0 = 45. lat0 = 45.
lon0 = 0. lon0 = 0.
...@@ -13,18 +50,24 @@ lon0 = 0. ...@@ -13,18 +50,24 @@ lon0 = 0.
sat_az = "180.0" sat_az = "180.0"
sat_zen = "45.0" sat_zen = "45.0"
obs_kind_nrs = obskind_read()
def degr_to_rad(degr): def degr_to_rad(degr):
"""Convert to DART convention = radians""" """Convert to DART convention = radians"""
if degr < 0: if degr < 0:
degr += 360 degr += 360
return degr/360*2*np.pi return degr/360*2*np.pi
def round_to_day(dtobj): def round_to_day(dtobj):
return dtobj.replace(second=0, minute=0, hour=0) return dtobj.replace(second=0, minute=0, hour=0)
def add_timezone_UTC(t): def add_timezone_UTC(t):
return dt.datetime(t.year, t.month, t.day, t.hour, t.minute, tzinfo=dt.timezone.utc) return dt.datetime(t.year, t.month, t.day, t.hour, t.minute, tzinfo=dt.timezone.utc)
def get_dart_date(time_dt): def get_dart_date(time_dt):
# assumes input is UTC! # assumes input is UTC!
days_since_010120 = 145731 days_since_010120 = 145731
...@@ -33,8 +76,9 @@ def get_dart_date(time_dt): ...@@ -33,8 +76,9 @@ def get_dart_date(time_dt):
secs_thatday = str(int((time_dt - round_to_day(time_dt)).total_seconds())) secs_thatday = str(int((time_dt - round_to_day(time_dt)).total_seconds()))
return dart_date_day, secs_thatday return dart_date_day, secs_thatday
def calc_obs_locations(n_obs, coords_from_domaincenter=True, def calc_obs_locations(n_obs, coords_from_domaincenter=True,
distance_between_obs_km=9999, folder_obs_coords=False): distance_between_obs_km=9999, fpath_coords=False):
"""Calculate coordinate pairs for locations """Calculate coordinate pairs for locations
Args: Args:
...@@ -94,62 +138,44 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True, ...@@ -94,62 +138,44 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True,
lons[skip+i*dx, skip+j*dx])) lons[skip+i*dx, skip+j*dx]))
try: try:
if fpath_obs_locations: if fpath_coords:
import pickle import pickle
os.makedirs(os.path.dirname(fpath_obs_locations), exist_ok=True) os.makedirs(os.path.dirname(fpath_coords), exist_ok=True)
with open(fpath_obs_locations, 'wb') as f: with open(fpath_coords, 'wb') as f:
pickle.dump(coords, f); print(fpath_obs_locations, 'saved.') pickle.dump(coords, f)
print(fpath_coords, 'saved.')
except Exception as e: except Exception as e:
warnings.warn(str(e)) warnings.warn(str(e))
assert len(coords) == n_obs, (len(coords), n_obs) assert len(coords) == n_obs, (len(coords), n_obs)
return coords return coords
def main_obspart(obs, last=False):
"""
Args:
obs (object)
last (bool): True if this is the last observation in the obs_seq file
"""
if last:
line_link = " "+str(obs.i-1)+" -1 -1"
else:
line_link = " -1 "+str(obs.i+1)+" -1"
return """ def write_file(msg, output_path='./'):
OBS """+str(obs.i)+""" try:
"""+line_link+""" os.remove(output_path)
obdef except OSError:
loc3d pass
"""+obs.lon_rad+" "+obs.lat_rad+" "+obs.vert+" "+obs.vert_coord_sys+"""
kind
"""+obs.kind_nr+"""
"""+obs.secs_thatday+""" """+obs.dart_date_day+"""
"""+obs.error_var
def write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords, with open(output_path, 'w') as f:
dart_date_day, secs_thatday, output_path, f.write(msg)
vert_coord_sfc=False): print(output_path, 'saved.')
"""
Args:
dart_date_day (str): DART internal time formatted date
secs_thatday (str): DART internal time of day (seconds since 0 UTC)
vert_coord_sfc (bool):
if True, then vertical coordinate is height above ground, i.e. "surface observation"
if False, then vertical is hgt_AMSL
"""
vert_coord_sys = 3 # meters AMSL
if vert_coord_sfc:
vert_coord_sys = -1
n_obs = len(coords) def append_hgt_to_coords(coords, heights):
n_obs_str = str(n_obs) coords2 = []
error_var = str(error_var) for i in range(len(coords)):
line_obstypedef = obs_kind_nr+' '+obs_name for hgt_m in heights:
vert_coord_sys = str(vert_coord_sys) coords2.append(coords[i] + (hgt_m,))
msg = """ n_obs = len(coords2)
obs_sequence print('effective number of observations (with vertical levels):', n_obs,
' on each level:', len(coords))
return coords2
def preamble(n_obs, line_obstypedef):
n_obs_str = str(n_obs)
return """ obs_sequence
obs_kind_definitions obs_kind_definitions
1 1
"""+line_obstypedef+""" """+line_obstypedef+"""
...@@ -157,283 +183,145 @@ obs_kind_definitions ...@@ -157,283 +183,145 @@ obs_kind_definitions
num_obs: """+n_obs_str+" max_num_obs: "+n_obs_str+""" num_obs: """+n_obs_str+" max_num_obs: "+n_obs_str+"""
first: 1 last: """+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] def write_section(obs, last=False):
lat = coords[i_obs-1][0] """
hgt_m = str(coords[i_obs-1][2]) Args:
obs (object)
last (bool): True if this is the last observation in the obs_seq file
"""
lon_rad = str(degr_to_rad(obs['lon']))
lat_rad = str(degr_to_rad(obs['lat']))
lon_rad = str(degr_to_rad(lon)) if last:
lat_rad = str(degr_to_rad(lat)) line_link = " "+str(obs['i']-1)+" -1 -1"
else:
line_link = " -1 "+str(obs['i']+1)+" -1"
# compile text return """
if i_obs < int(n_obs): OBS """+str(obs['i'])+"""
msg += """ """+line_link+"""
OBS """+str(i_obs)+"""
-1 """+str(i_obs+1)+""" -1
obdef
loc3d
"""+lon_rad+" "+lat_rad+" "+hgt_m+" "+vert_coord_sys+"""
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 obdef
loc3d loc3d
"""+lon_rad+" "+lat_rad+" "+hgt_m+" "+vert_coord_sys+""" """+lon_rad+" "+lat_rad+" "+str(obs['vert_coord'])+" "+obs['vert_coord_sys']+"""
kind kind
"""+obs_kind_nr+""" """+obs['kind_nr']+"""
"""+secs_thatday+""" """+dart_date_day+""" """+obs['appendix']+"""
"""+error_var """+obs['secs_thatday']+""" """+obs['dart_date_day']+"""
"""+str(obs['obserr_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, coords, error_var, output_path='./'): def create_obsseq_in(time_dt, obscfg, zero_error=False,
archive_obs_coords=False):
"""Create obs_seq.in """Create obs_seq.in
Args: Args:
time_dt (dt.datetime): time of observation time_dt (dt.datetime): time of observation
obscfg (dict)
archive_obs_coords (str, False): path to folder
channel_id (int): SEVIRI channel number channel_id (int): SEVIRI channel number
see https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov12/ir_srf/rtcoef_msg_4_seviri_srf.html see https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov12/ir_srf/rtcoef_msg_4_seviri_srf.html
coords (list of 2-tuples with (lat,lon)) coords (list of 2-tuples with (lat,lon))
error_var (float): obserr_std (np.array): shape (n_obs,), one value for each observation,
gaussian error with this variance is added to the truth at observation time gaussian error with this std-dev is added to the truth at observation time
output_path (str): directory where `obs_seq.in` will be saved output_path (str): directory where `obs_seq.in` will be saved
""" """
# time_dt = add_timezone_UTC(time_dt)
# time_dt = dt.datetime(2008, 7, 30, 15, 30, tzinfo=dt.timezone.utc)
error_var = str(error_var)
n_obs = len(coords)
# Brightness temperature or Reflectance?
channel_id = int(channel_id)
if channel_id in [1, 2, 3, 12]:
line_obstypedef = ' 256 MSG_4_SEVIRI_BDRF'
code = '256'
else:
line_obstypedef = ' 255 MSG_4_SEVIRI_TB'
code = '255'
channel_id = str(channel_id)
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))
print('sunzen', sun_zen, 'sunazi', sun_az)
dart_date_day, secs_thatday = get_dart_date(time_dt)
print('secs, days:', secs_thatday, dart_date_day)
n_obs_str = str(int(n_obs))
error_var = str(error_var)
msg = """ n_obs = obscfg['n_obs']
obs_sequence coords = calc_obs_locations(n_obs,
obs_kind_definitions coords_from_domaincenter=False,
1 distance_between_obs_km=obscfg.get('distance_between_obs_km', False),
"""+line_obstypedef+""" fpath_coords=archive_obs_coords)
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]
lon_rad = str(degr_to_rad(lon))
lat_rad = str(degr_to_rad(lat))
# compile text kind = obscfg['kind']
if i_obs < int(n_obs): sat_channel = obscfg.get('sat_channel', False)
msg += """
OBS """+str(i_obs)+"""
-1 """+str(i_obs+1)+""" -1
obdef
loc3d
"""+lon_rad+""" """+lat_rad+""" -888888.0000000000 -2
kind
"""+code+"""
visir
"""+sat_az+""" """+sat_zen+""" """+sun_az+"""
"""+sun_zen+"""
12 4 21 """+channel_id+"""
-888888.000000000
1
"""+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+""" -888888.0000000000 -2
kind
"""+code+"""
visir
"""+sat_az+""" """+sat_zen+""" """+sun_az+"""
"""+sun_zen+"""
12 4 21 """+channel_id+"""
-888888.000000000
1
"""+secs_thatday+""" """+dart_date_day+"""
"""+error_var
fpath = output_path+'/obs_seq.in' # obs error has to be len(n_obs)
obserr_std = np.zeros(n_obs) if zero_error else obscfg['err_std']
try: try:
os.remove(fpath) assert len(obserr_std) == n_obs # fails for scalars
except OSError: except TypeError:
pass obserr_std = np.ones(n_obs) * obserr_std
with open(fpath, 'w') as f: # determine vertical coordinates
f.write(msg) if not sat_channel:
print(fpath, 'saved.') if 'SYNOP' in kind:
vert_coord_sys = "-1" # meters AGL
vert_coords = [2, ]
def calc_obs_locations_3d(coords, heights):
# append height
coords2 = []
for i in range(len(coords)):
for hgt_m in 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))
return coords2
def generic_obs(obs_kind, time_dt, coords, error_var, heights=False, output_path='./'):
obs_kind_nrs = {'RADIOSONDE_TEMPERATURE': '5',
'RADAR_REFLECTIVITY': '37',
'SYNOP_SURFACE_PRESSURE': '94',
'SYNOP_SPECIFIC_HUMIDITY': '95',
'SYNOP_TEMPERATURE': '96',
}
if 'SYNOP' in obs_kind:
is_sfc_obs = True
heights = [2,]
else: else:
is_sfc_obs = False vert_coord_sys = "3" # meters AMSL
vert_coords = obscfg['heights']
else:
vert_coord_sys = "-2" # undefined height
vert_coords = ["-888888.0000000000", ]
if not heights: coords = append_hgt_to_coords(coords, vert_coords)
heights = [5000., ] obs_kind_nr = obs_kind_nrs[kind]
coords = calc_obs_locations_3d(coords, heights) line_obstypedef = ' '+obs_kind_nr+' '+kind
dart_date_day, secs_thatday = get_dart_date(add_timezone_UTC(time_dt)) time_dt = add_timezone_UTC(time_dt)
dart_date_day, secs_thatday = get_dart_date(time_dt)
print('secs, days:', secs_thatday, dart_date_day) print('secs, days:', secs_thatday, dart_date_day)
obs_kind_nr = obs_kind_nrs[obs_kind] if sat_channel:
write_generic_obsseq(obs_kind, obs_kind_nr, error_var, coords, sun_az = str(get_azimuth(lat0, lon0, time_dt))
dart_date_day, secs_thatday, output_path, sun_zen = str(90. - get_altitude(lat0, lon0, time_dt))
vert_coord_sfc=is_sfc_obs) print('sunzen', sun_zen, 'sunazi', sun_az)
def obskind_read():
raw_obskind_dart = """
5 RADIOSONDE_TEMPERATURE
6 RADIOSONDE_SPECIFIC_HUMIDITY
12 AIRCRAFT_U_WIND_COMPONENT
13 AIRCRAFT_V_WIND_COMPONENT
14 AIRCRAFT_TEMPERATURE
16 ACARS_U_WIND_COMPONENT
17 ACARS_V_WIND_COMPONENT
18 ACARS_TEMPERATURE
29 LAND_SFC_PRESSURE
30 SAT_U_WIND_COMPONENT
31 SAT_V_WIND_COMPONENT
36 DOPPLER_RADIAL_VELOCITY
37 RADAR_REFLECTIVITY
83 GPSRO_REFRACTIVITY
94 SYNOP_SURFACE_PRESSURE
95 SYNOP_SPECIFIC_HUMIDITY
96 SYNOP_TEMPERATURE
254 MSG_4_SEVIRI_RADIANCE
255 MSG_4_SEVIRI_TB
256 MSG_4_SEVIRI_BDRF"""
# lookup table for kind nr
alist = raw_obskind_dart.split()
assert len(alist) % 2 == 0, alist
obskind_nrs = {}
for i in range(0, len(alist)-1, 2):
obskind_nrs[alist[i+1]] = alist[i]
return obskind_nrs
obskind_nrs = obskind_read()
def create_obsseq_in(obscfg, obserr_var):
"""
Args:
obserr_var (np.array): observation error variance
shape (n_obs,), one value for each observation,
"""
self.coords = osq.calc_obs_locations(obscfg['n_obs'],
coords_from_domaincenter=False,
distance_between_obs_km=obscfg.get('distance_between_obs_km', False),
fpath_obs_locations=folder_obs_coords+'/obs_coords_'+obscfg['kind']+'.pkl')
# for i_obs in obscfg['n_obs']:
# instruction = dict(kind_nr = obskind_nrs[obscfg['kind']],
# sat_channel = obscfg.get('sat_channel', False),
# heights = obscfg.get('heights', False),
# obs_kinds, time_dt, coords, error_var, heights=False, output_path='./'):
if 'SYNOP' in obs_kind: appendix = """visir
is_sfc_obs = True """+sat_az+""" """+sat_zen+""" """+sun_az+"""
heights = [2,] """+sun_zen+"""
12 4 21 """+str(sat_channel)+"""
-888888.000000000
1"""
else: else:
is_sfc_obs = False appendix = ''
if not heights: txt = preamble(n_obs, line_obstypedef)
heights = [5000., ]
coords = calc_obs_locations_3d(coords, heights)
dart_date_day, secs_thatday = get_dart_date(add_timezone_UTC(time_dt)) for i_obs in range(int(n_obs)):
print('secs, days:', secs_thatday, dart_date_day) last = False
if i_obs == int(n_obs)-1:
last = True # last_observation
obs_kind_nr = obs_kind_nrs[obs_kind] txt += write_section(dict(i=i_obs+1,
kind_nr=obs_kind_nr,
dart_date_day=dart_date_day,
secs_thatday=secs_thatday,
lon=coords[i_obs][1],
lat=coords[i_obs][0],
vert_coord=coords[i_obs][2],
vert_coord_sys=vert_coord_sys,
obserr_var=obserr_std[i_obs]**2,
appendix=appendix),
last=last)
for obs_kind in obs_kinds: write_file(txt, output_path=cluster.dartrundir+'/obs_seq.in')
write_generic_obsseq2(obs_kind, obs_kind_nr, error_var, coords,
dart_date_day, secs_thatday, output_path,
vert_coord_sfc=is_sfc_obs)
if __name__ == '__main__': if __name__ == '__main__':
time_dt = dt.datetime(2008, 7, 30, 10, 0) time_dt = dt.datetime(2008, 7, 30, 10, 0)
n_obs = 64 n_obs = 64 # radar: n_obs for each observation height level
sat_channel = 1
vis = dict(kind='MSG_4_SEVIRI_BDRF',
distance_between_obs_meters = 10000 sat_channel=1, n_obs=n_obs, err_std=0.03,
#error_var = 0.001 cov_loc_radius_km=10)
obs_coords = calc_obs_locations(n_obs, coords_from_domaincenter=False, wv = dict(kind='MSG_4_SEVIRI_TB',
distance_between_obs_km=distance_between_obs_meters, sat_channel=6, n_obs=n_obs, err_std=False,
fpath_obs_locations=None) cov_loc_radius_km=10)
#sat(time_dt, sat_channel, obs_coords, error_var, output_path='./') ir108 = dict(kind='MSG_4_SEVIRI_TB',
sat_channel=9, n_obs=n_obs, err_std=5.,
# error_var = (5.)**2 cov_loc_radius_km=10)
# generic_obs('RADAR_REFLECTIVITY', time_dt, obs_coords, error_var, heights=[5000.,], output_path='./')
radar = dict(kind='RADAR', n_obs=n_obs, err_std=5.,
error_var = (0.5)**2 heights=np.arange(1000, 15001, 1000),
generic_obs('RADIOSONDE_TEMPERATURE', time_dt, obs_coords, error_var, heights=[5000.,]) cov_loc_radius_km=10, cov_loc_vert_km=2)
t2m = dict(kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=1.0,
cov_loc_radius_km=32, cov_loc_vert_km=1)
psfc = dict(kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
cov_loc_radius_km=32, cov_loc_vert_km=5)
create_obsseq_in(time_dt, vis, archive_obs_coords='./coords_stage1.pkl')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment