diff --git a/config/cfg.py b/config/cfg.py index 8451eba854eeb7d897515d7dc8f53d8418a679ab..2802e07aee382e1555e3b5f7a0e6896b1aa8ab1c 100755 --- a/config/cfg.py +++ b/config/cfg.py @@ -17,11 +17,14 @@ exp.n_nodes = 10 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) -wv = dict(sat_channel=6, n_obs=n_obs, err_std=False, - cov_loc_radius_km=10) -ir108 = dict(sat_channel=9, n_obs=n_obs, err_std=5., +wv = dict(kind='MSG_4_SEVIRI_TB', + sat_channel=6, n_obs=n_obs, err_std=False, + cov_loc_radius_km=10) +ir108 = dict(kind='MSG_4_SEVIRI_TB', + sat_channel=9, n_obs=n_obs, err_std=5., cov_loc_radius_km=10) radar = dict(kind='RADAR', n_obs=n_obs, err_std=5., diff --git a/scripts/assim_synth_obs.py b/scripts/assim_synth_obs.py index e1e292b7cc99c080421125ef7a0bda505cc63cc1..a41e8c5c5f6f1e8bdfe9035daa0e6f1b8b4b7dbc 100755 --- a/scripts/assim_synth_obs.py +++ b/scripts/assim_synth_obs.py @@ -88,17 +88,14 @@ def edit_obserr_in_obsseq(fpath_obsseqin, OEs): f.write(line) -def set_input_nml(sat_channel=False, just_prior_values=False, - cov_loc_radius_km=32, cov_loc_vert_km=False): +def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False): """descr""" cov_loc_radian = cov_loc_radius_km/earth_radius_km - if just_prior_values: - template = cluster.scriptsdir+'/../templates/input.prioronly.nml' - else: - template = cluster.scriptsdir+'/../templates/input.nml' - copy(template, cluster.dartrundir+'/input.nml') + copy(cluster.scriptsdir+'/../templates/input.nml', + cluster.dartrundir+'/input.nml') + # options are overwritten with settings options = {'<n_ens>': str(int(exp.n_ens)), '<cov_loc_radian>': str(cov_loc_radian)} @@ -160,7 +157,6 @@ def obs_operator_ensemble(): def obs_operator_nature(): prepare_nature_dart() - set_input_nml(sat_channel=sat_channel) os.chdir(cluster.dartrundir) os.remove(cluster.dartrundir+'/obs_seq.out') os.system('mpirun -np 12 ./perfect_model_obs') @@ -194,11 +190,7 @@ def calc_obserr_WV73(Hx_nature, Hx_prior): OEs[iobs] = oe_nature return OEs - 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) os.chdir(cluster.dartrundir) os.remove(cluster.dartrundir+'/obs_seq.out') @@ -206,7 +198,6 @@ def generate_observations(): raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir) os.system('mpirun -np 12 ./perfect_model_obs') - def assimilate(): os.chdir(cluster.dartrundir) os.remove(cluster.dartrundir+'/obs_seq.final') @@ -214,17 +205,14 @@ def assimilate(): raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir) os.system('mpirun -np 48 ./filter') -def archive_diagnostics(time): +def archive_diagnostics(archive_stage): print('archive obs space diagnostics') - savedir = cluster.archivedir()+'/obs_seq_final/' - mkdir(savedir) - copy(cluster.dartrundir+'/obs_seq.final', savedir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final')) + mkdir(archive_stage) + copy(cluster.dartrundir+'/obs_seq.final', archive_stage+'/obs_seq.final') try: print('archive regression diagnostics') - savedir = cluster.archivedir()+'/reg_factor/' - mkdir(savedir) - copy(cluster.dartrundir+'/reg_diagnostics', savedir+time.strftime('/%Y-%m-%d_%H:%M_reg_diagnostics')) + copy(cluster.dartrundir+'/reg_diagnostics', archive_stage+'/reg_diagnostics') except Exception as e: warnings.warn(str(e)) @@ -233,27 +221,20 @@ def recycle_output(): for iens in range(1, exp.n_ens+1): os.rename(cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4), cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01') - - -def archive_output_mean(time, istage): - istage = str(istage) - print('copy prior posterior members to archive') +def archive_output_mean(archive_stage): 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) - copy(cluster.dartrundir+'/input.nml', - cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/input '+istage+'.nml')) + copy(cluster.dartrundir+'/input.nml', archive_stage+'/input.nml') # filter_in = cluster.dartrundir+'/preassim_member_'+str(iens).zfill(4)+'.nc' # filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4) # copy mean and sd to archive - print('copy preassim, postassim mean and sd') for f in ['output_mean.nc', 'output_sd.nc']: - copy(cluster.dartrundir+'/'+f, - cluster.archivedir()+'/'+f[:-3]+time.strftime('_'+istage+'_%Y-%m-%d_%H:%M:%S')) + copy(cluster.dartrundir+'/'+f, archive_stage+'/'+f) if __name__ == "__main__": @@ -279,20 +260,19 @@ if __name__ == "__main__": """ 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.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) for istage, obscfg in enumerate(exp.observations): + kind = obscfg['kind'] + archive_stage = archive_time + '/stage'+str(istage)+'_'+kind n_obs = obscfg['n_obs'] sat_channel = obscfg.get('sat_channel', False) - obs_coords = osq.calc_obs_locations(n_obs, - coords_from_domaincenter=False, - distance_between_obs_km=obscfg.get('distance_between_obs_km', False), - folder_obs_coords=folder_obs_coords) + obscfg['folder_obs_coords'] = archive_stage+'/obs_coords.pkl' set_DART_nml(sat_channel=sat_channel, cov_loc=obscfg['cov_loc_radius_km'], @@ -303,25 +283,25 @@ if __name__ == "__main__": if sat_channel != 6: 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_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: - 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() assimilate() - archive_diagnostics(time) + archive_diagnostics(archive_stage) if istage < n_stages-1: # recirculation: filter output -> input recycle_output() - archive_output_mean(time, istage) + archive_output_mean(archive_stage) elif istage == n_stages-1: # last assimilation, continue integration now diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py index 1d4050496c3f5d36afda45ed9728ba63a475c64e..ec0496b2dbc7e1f03ea4238bd462d940a78314d7 100755 --- a/scripts/create_obsseq.py +++ b/scripts/create_obsseq.py @@ -4,8 +4,45 @@ import os, sys, warnings import numpy as np import datetime as dt +from config.cfg import exp, cluster 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 lat0 = 45. lon0 = 0. @@ -13,18 +50,24 @@ lon0 = 0. sat_az = "180.0" sat_zen = "45.0" +obs_kind_nrs = obskind_read() + + def degr_to_rad(degr): """Convert to DART convention = radians""" if degr < 0: degr += 360 return degr/360*2*np.pi + def round_to_day(dtobj): return dtobj.replace(second=0, minute=0, hour=0) + def add_timezone_UTC(t): return dt.datetime(t.year, t.month, t.day, t.hour, t.minute, tzinfo=dt.timezone.utc) + def get_dart_date(time_dt): # assumes input is UTC! days_since_010120 = 145731 @@ -33,8 +76,9 @@ 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, 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 Args: @@ -94,154 +138,44 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True, lons[skip+i*dx, skip+j*dx])) try: - if fpath_obs_locations: + if fpath_coords: import pickle - os.makedirs(os.path.dirname(fpath_obs_locations), exist_ok=True) - with open(fpath_obs_locations, 'wb') as f: - pickle.dump(coords, f); print(fpath_obs_locations, 'saved.') + os.makedirs(os.path.dirname(fpath_coords), exist_ok=True) + with open(fpath_coords, 'wb') as f: + pickle.dump(coords, f) + print(fpath_coords, 'saved.') except Exception as e: warnings.warn(str(e)) assert len(coords) == n_obs, (len(coords), n_obs) 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 """ - OBS """+str(obs.i)+""" -"""+line_link+""" -obdef -loc3d - """+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, - dart_date_day, secs_thatday, output_path, - vert_coord_sfc=False): - """ - 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) - n_obs_str = str(n_obs) - error_var = str(error_var) - line_obstypedef = obs_kind_nr+' '+obs_name - vert_coord_sys = str(vert_coord_sys) - - 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]) - - 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+" "+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 -loc3d - """+lon_rad+" "+lat_rad+" "+hgt_m+" "+vert_coord_sys+""" -kind - """+obs_kind_nr+""" - """+secs_thatday+""" """+dart_date_day+""" - """+error_var - fpath = output_path+'/obs_seq.in' +def write_file(msg, output_path='./'): try: - os.remove(fpath) + os.remove(output_path) except OSError: pass - with open(fpath, 'w') as f: + with open(output_path, 'w') as f: f.write(msg) - print(fpath, 'saved.') + print(output_path, 'saved.') -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 - 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 - """ - # 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) +def append_hgt_to_coords(coords, heights): + coords2 = [] + for i in range(len(coords)): + for hgt_m in heights: + coords2.append(coords[i] + (hgt_m,)) - dart_date_day, secs_thatday = get_dart_date(time_dt) - print('secs, days:', secs_thatday, dart_date_day) + n_obs = len(coords2) + print('effective number of observations (with vertical levels):', n_obs, + ' on each level:', len(coords)) + return coords2 - n_obs_str = str(int(n_obs)) - error_var = str(error_var) - msg = """ - obs_sequence +def preamble(n_obs, line_obstypedef): + n_obs_str = str(n_obs) + return """ obs_sequence obs_kind_definitions 1 """+line_obstypedef+""" @@ -249,191 +183,145 @@ obs_kind_definitions 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)) +def write_section(obs, last=False): + """ + 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'])) - # 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+""" -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 + if last: + line_link = " "+str(obs['i']-1)+" -1 -1" + else: + line_link = " -1 "+str(obs['i']+1)+" -1" + + return """ + OBS """+str(obs['i'])+""" +"""+line_link+""" obdef loc3d - """+lon_rad+""" """+lat_rad+""" -888888.0000000000 -2 + """+lon_rad+" "+lat_rad+" "+str(obs['vert_coord'])+" "+obs['vert_coord_sys']+""" 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' - try: - os.remove(fpath) - except OSError: - pass - - with open(fpath, 'w') as f: - f.write(msg) - print(fpath, 'saved.') + """+obs['kind_nr']+""" +"""+obs['appendix']+""" +"""+obs['secs_thatday']+""" """+obs['dart_date_day']+""" +"""+str(obs['obserr_var']) -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 create_obsseq_in(time_dt, obscfg, zero_error=False, + archive_obs_coords=False): + """Create obs_seq.in + Args: + time_dt (dt.datetime): time of observation + obscfg (dict) + archive_obs_coords (str, False): path to folder + + channel_id (int): SEVIRI channel number + 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)) + obserr_std (np.array): shape (n_obs,), one value for each observation, + 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 + """ -def generic_obs(obs_kind, time_dt, coords, error_var, heights=False, output_path='./'): + n_obs = obscfg['n_obs'] + coords = calc_obs_locations(n_obs, + coords_from_domaincenter=False, + distance_between_obs_km=obscfg.get('distance_between_obs_km', False), + fpath_coords=archive_obs_coords) - obs_kind_nrs = {'RADIOSONDE_TEMPERATURE': '5', - 'RADAR_REFLECTIVITY': '37', - 'SYNOP_SURFACE_PRESSURE': '94', - 'SYNOP_SPECIFIC_HUMIDITY': '95', - 'SYNOP_TEMPERATURE': '96', - } + kind = obscfg['kind'] + sat_channel = obscfg.get('sat_channel', False) - if 'SYNOP' in obs_kind: - is_sfc_obs = True - heights = [2,] + # obs error has to be len(n_obs) + obserr_std = np.zeros(n_obs) if zero_error else obscfg['err_std'] + try: + assert len(obserr_std) == n_obs # fails for scalars + except TypeError: + obserr_std = np.ones(n_obs) * obserr_std + + # determine vertical coordinates + if not sat_channel: + if 'SYNOP' in kind: + vert_coord_sys = "-1" # meters AGL + vert_coords = [2, ] + else: + vert_coord_sys = "3" # meters AMSL + vert_coords = obscfg['heights'] else: - is_sfc_obs = False + vert_coord_sys = "-2" # undefined height + vert_coords = ["-888888.0000000000", ] - if not heights: - heights = [5000., ] - coords = calc_obs_locations_3d(coords, heights) + coords = append_hgt_to_coords(coords, vert_coords) + obs_kind_nr = obs_kind_nrs[kind] + 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) - obs_kind_nr = obs_kind_nrs[obs_kind] - write_generic_obsseq(obs_kind, obs_kind_nr, error_var, coords, - dart_date_day, secs_thatday, output_path, - vert_coord_sfc=is_sfc_obs) + if sat_channel: + 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) -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: - is_sfc_obs = True - heights = [2,] + appendix = """visir + """+sat_az+""" """+sat_zen+""" """+sun_az+""" + """+sun_zen+""" + 12 4 21 """+str(sat_channel)+""" + -888888.000000000 + 1""" else: - is_sfc_obs = False + appendix = '' - if not heights: - heights = [5000., ] - coords = calc_obs_locations_3d(coords, heights) + txt = preamble(n_obs, line_obstypedef) - dart_date_day, secs_thatday = get_dart_date(add_timezone_UTC(time_dt)) - print('secs, days:', secs_thatday, dart_date_day) + for i_obs in range(int(n_obs)): + 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__': time_dt = dt.datetime(2008, 7, 30, 10, 0) - n_obs = 64 - sat_channel = 1 - - distance_between_obs_meters = 10000 - #error_var = 0.001 - obs_coords = calc_obs_locations(n_obs, coords_from_domaincenter=False, - distance_between_obs_km=distance_between_obs_meters, - fpath_obs_locations=None) - #sat(time_dt, sat_channel, obs_coords, error_var, output_path='./') - - # error_var = (5.)**2 - # generic_obs('RADAR_REFLECTIVITY', time_dt, obs_coords, error_var, heights=[5000.,], output_path='./') - - error_var = (0.5)**2 - generic_obs('RADIOSONDE_TEMPERATURE', time_dt, obs_coords, error_var, heights=[5000.,]) + n_obs = 64 # radar: n_obs for each observation height level + + vis = dict(kind='MSG_4_SEVIRI_BDRF', + sat_channel=1, n_obs=n_obs, err_std=0.03, + cov_loc_radius_km=10) + wv = dict(kind='MSG_4_SEVIRI_TB', + sat_channel=6, n_obs=n_obs, err_std=False, + cov_loc_radius_km=10) + ir108 = dict(kind='MSG_4_SEVIRI_TB', + sat_channel=9, n_obs=n_obs, err_std=5., + cov_loc_radius_km=10) + + radar = dict(kind='RADAR', n_obs=n_obs, err_std=5., + heights=np.arange(1000, 15001, 1000), + 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')