From 0da787c84dbf909247bf68a6506ade971bf5006c Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Mon, 13 Sep 2021 19:05:49 +0200
Subject: [PATCH] bugfix
---
scripts/assim_synth_obs.py | 99 +++++++++++++++++++++-----------------
scripts/create_obsseq.py | 16 +++---
2 files changed, 65 insertions(+), 50 deletions(-)
diff --git a/scripts/assim_synth_obs.py b/scripts/assim_synth_obs.py
index b15eec7..df29afc 100755
--- a/scripts/assim_synth_obs.py
+++ b/scripts/assim_synth_obs.py
@@ -74,6 +74,7 @@ def replace_errors_obsseqout(f, new_errors):
new_errors (np.array) : standard deviation,
shape must match the number of observations
"""
+ debug = True
obsseq = open(f, 'r').readlines()
# find number of lines between two ' OBS ' lines:
@@ -96,7 +97,7 @@ def replace_errors_obsseqout(f, new_errors):
previous_err_var = obsseq[line_error_obs_i]
new_err_obs_i = new_errors[i_obs]**2 # variance in obs_seq.out
- #print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
+ if debug: print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n'
i_obs += 1 # next iteration
@@ -189,9 +190,10 @@ def obs_operator_ensemble(istage):
# DART may need a wrfinput file as well, which serves as a template for dimension sizes
symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
+ # I dont think this is necessary, we do this already in pre_assim.py
# add geodata, if istage>0, wrfout is DART output (has coords)
- if istage == 0:
- wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
+ #if istage == 0:
+ # wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
# run perfect_model obs (forward operator)
os.system('mpirun -np 12 ./perfect_model_obs > /dev/null')
@@ -221,6 +223,7 @@ def link_nature_to_dart_truth(time):
cluster.dartrundir+'/wrfout_d01')
# DART may need a wrfinput file as well, which serves as a template for dimension sizes
symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
+ print('linked', time.strftime(cluster.nature_wrfout), 'to', cluster.dartrundir+'/wrfout_d01')
def prepare_nature_dart(time):
@@ -254,7 +257,7 @@ def run_perfect_model_obs():
os.system('mpirun -np 12 ./perfect_model_obs > log.perfect_model_obs')
if not os.path.exists(cluster.dartrundir+'/obs_seq.out'):
raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir,
- '\n look for '+cluster.dartrundir+'log.perfect_model_obs')
+ '\n look for '+cluster.dartrundir+'/log.perfect_model_obs')
def assimilate(nproc=96):
print('time now', dt.datetime.now())
@@ -279,6 +282,8 @@ def archive_diagnostics(archive_dir, time):
# warnings.warn(str(e))
def recycle_output():
+ """Use output of assimilation (./filter) as input for another assimilation (with ./filter)
+ Specifically, this copies the state fields from filter_restart_d01.000x to the wrfout files in advance_temp folders"""
update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'QRAIN', 'U10', 'V10', 'T2', 'Q2', 'TSK', 'PSFC', 'CLDFRA']
updates = ','.join(update_vars)
@@ -344,6 +349,9 @@ if __name__ == "__main__":
Assumptions:
- x_ensemble is already linked for DART to advance_temp<iens>/wrfout_d01
+
+ Example call:
+ python assim.py 2008-08-07_12:00
"""
time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
@@ -360,20 +368,21 @@ if __name__ == "__main__":
################################################
print(' 1) get the assimilation errors in a single vector ')
- error_assimilate = []
+ error_generate = []
for i, obscfg in enumerate(exp.observations):
n_obs = obscfg['n_obs']
- n_obs_z = obscfg.get('heights', 1)
+ n_obs_z = len(obscfg.get('heights', [1,]))
n_obs_3d = n_obs * n_obs_z
parametrized = obscfg.get('sat_channel') == 6
+
if not parametrized:
err_this_type = np.zeros(n_obs_3d) + obscfg['error_generate']
else: # error parametrization for WV73
# get observations for sat 6
- osq.create_obsseqin_alltypes(time, [obscfg,], obs_errors='error_generate')
+ osq.create_obsseqin_alltypes(time, [obscfg,], np.zeros(n_obs_3d))
run_perfect_model_obs()
# depends on obs_seq.out produced before by run_perfect_model_obs()
@@ -382,11 +391,11 @@ if __name__ == "__main__":
Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
err_this_type = calc_obserr_WV73(Hx_nat, Hx_prior)
- error_assimilate.extend(err_this_type)
+ error_generate.extend(err_this_type) # the obs-error we assume for generating observations
################################################
print(' 2) generate observations ')
- osq.create_obsseqin_alltypes(time, exp.observations, obs_errors='error_generate',
+ osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=error_generate,
archive_obs_coords=archive_stage+'/obs_coords.pkl')
first_obstype = exp.observations[0]
@@ -397,58 +406,60 @@ if __name__ == "__main__":
################################################
print(' 3) assimilate with adapted observation-errors ')
+ error_assimilate = np.zeros(n_obs_3d) + obscfg['error_assimilate'] # the obs-error we assume for assimilation
replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
- t = time.time()
+ t = time_module.time()
assimilate()
+ print('filter took', time_module.time()-t, 'seconds')
dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage)
archive_diagnostics(dir_obsseq, time)
archive_output(archive_stage)
- sys.exit() # below is the code for separate assimilation of different obs types
+ #sys.exit() # below is the code for separate assimilation of different obs types
- for istage, obscfg in enumerate(exp.observations):
- print('running observation stage', istage, obscfg)
+ #for istage, obscfg in enumerate(exp.observations):
+ # print('running observation stage', istage, obscfg)
- archive_stage = archive_time + '/assim_stage'+str(istage)+'/'
- n_obs = obscfg['n_obs']
- n_obs_3d = n_obs * len(obscfg['heights'])
- sat_channel = obscfg.get('sat_channel', False)
+ # archive_stage = archive_time + '/assim_stage'+str(istage)+'/'
+ # n_obs = obscfg['n_obs']
+ # n_obs_3d = n_obs * len(obscfg['heights'])
+ # sat_channel = obscfg.get('sat_channel', False)
- # debug option: overwrite time in prior files
- # for iens in range(1,41):
- # os.system('ncks -A -v Times '+cluster.dartrundir+'/wrfout_d01 '+cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01')
+ # # debug option: overwrite time in prior files
+ # # for iens in range(1,41):
+ # # os.system('ncks -A -v Times '+cluster.dartrundir+'/wrfout_d01 '+cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01')
- if error_assimilate == False:
- # use error parametrization for assimilation-obs.errors
+ # if error_assimilate == False:
+ # # use error parametrization for assimilation-obs.errors
- if sat_channel != 6:
- raise NotImplementedError('sat channel '+str(sat_channel))
+ # if sat_channel != 6:
+ # raise NotImplementedError('sat channel '+str(sat_channel))
- # depends on obs_seq.out produced before by run_perfect_model_obs()
- Hx_nat, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out')
+ # # depends on obs_seq.out produced before by run_perfect_model_obs()
+ # Hx_nat, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out')
- Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
- error_assimilate = calc_obserr_WV73(Hx_nat, Hx_prior)
+ # Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
+ # error_assimilate = calc_obserr_WV73(Hx_nat, Hx_prior)
- else:
- # overwrite error values in obs_seq.out
- error_assimilate = np.zeros(n_obs_3d) + error_assimilate # ensure shape
+ # else:
+ # # overwrite error values in obs_seq.out
+ # error_assimilate = np.zeros(n_obs_3d) + error_assimilate # ensure shape
- replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
- assimilate()
- dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage)
- archive_diagnostics(dir_obsseq, time)
+ # replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
+ # assimilate()
+ # dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage)
+ # archive_diagnostics(dir_obsseq, time)
- if istage < n_stages-1:
- # recirculation: filter output -> input
- archive_output(archive_stage)
- recycle_output()
+ # if istage < n_stages-1:
+ # # recirculation: filter output -> input
+ # archive_output(archive_stage)
+ # recycle_output()
- elif istage == n_stages-1:
- # last assimilation, continue integration now
- archive_output(archive_stage)
+ # elif istage == n_stages-1:
+ # # last assimilation, continue integration now
+ # archive_output(archive_stage)
- else:
- RuntimeError('this should never occur?!')
+ # else:
+ # RuntimeError('this should never occur?!')
diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py
index 2b54b79..1d06ac7 100755
--- a/scripts/create_obsseq.py
+++ b/scripts/create_obsseq.py
@@ -143,7 +143,7 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True,
omit_covloc_radius_on_boundary = True
if omit_covloc_radius_on_boundary: # in order to avoid an increment step on the boundary
- skip_km = cov_loc_radius_km*1.5
+ skip_km = 50 # cov_loc_radius_km*1.5
skip_gridpoints = int(skip_km/model_dx_km) # skip this many gridpoints on each side
gridpoints_left = nx - 2*skip_gridpoints
@@ -376,13 +376,14 @@ def create_obsseq_in_separate_obs(time_dt, obscfg, obs_errors=False,
write_file(txt, output_path=cluster.dartrundir+'/obs_seq.in')
-def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors='error_assimilate', archive_obs_coords=False):
+def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors, archive_obs_coords=False):
"""Create obs_seq.in with multiple obs types in one file
Args:
time_dt (dt.datetime): time of observation
list_obscfg (list of dict)
- obs_errors (False or str): Key to obsdict, which field contains the error values
+ obs_errors (list of float, False): contains observation errors, one for each observation
+ if False: use zero error
archive_obs_coords (bool, str): False or str (filepath where `obs_seq.in` will be saved)
"""
print('creating obs_seq.in:')
@@ -411,9 +412,12 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors='error_assimilate'
n_obs_3d_thistype = len(coords)
# define obs error
- obserr_std = np.zeros(n_obs_3d_thistype)
- if obs_errors:
- obserr_std += obscfg[obs_errors]
+ if obs_errors == False:
+ obs_errors = np.zeros(n_obs_3d_thistype)
+ assert len(obs_errors) == n_obs_3d_thistype, 'len(obs_errors) == n_obs_3d_thistype'
+ obserr_std = obs_errors #np.zeros(n_obs_3d_thistype)
+ #if obs_errors:
+ # obserr_std += obscfg[obs_errors]
sat_info = write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt)
--
GitLab