diff --git a/config/exp_template.py b/config/exp_template.py
index 37b4a1534a7d2a4430c95f64b1dafe5ab91eaf63..0c97c5191a86ba237bb2befb2986bad77c826597 100755
--- a/config/exp_template.py
+++ b/config/exp_template.py
@@ -29,6 +29,9 @@ exp.dart_nml = {'&assim_tools_nml':
output_sd='.true.',
stages_to_write='output',
),
+ '&quality_control_nml':
+ dict(outlier_threshold='-1',
+ ),
'&location_nml':
dict(horiz_dist_only='.true.',
),
diff --git a/dartwrf/assim_synth_obs.py b/dartwrf/assim_synth_obs.py
index a23f804e773979f543176f5ba17bc17105711163..5b335443d88785eda7dbf610e841cab89bdad4cd 100755
--- a/dartwrf/assim_synth_obs.py
+++ b/dartwrf/assim_synth_obs.py
@@ -285,8 +285,11 @@ def qc_obs(time, oso):
print('saved', f_out_dart)
-def evaluate(assim_time, obs_seq_out=False,
- output_format="%Y-%m-%d_%H:%M_obs_seq.final-eval_posterior_allobs"):
+def evaluate(assim_time,
+ obs_seq_out=False,
+ prior=True,
+ posterior=False,
+ output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate"):
"""Calculates either prior or posterior obs space values.
Note: Depends on a prepared input_list.txt, which defines the ensemble (prior or posterior).
@@ -301,6 +304,8 @@ def evaluate(assim_time, obs_seq_out=False,
Returns:
obsseq.ObsSeq
"""
+ if prior == posterior:
+ raise ValueError('either prior or posterior must be True, the other must be False')
os.makedirs(cluster.dart_rundir, exist_ok=True) # create directory to run DART in
os.chdir(cluster.dart_rundir)
@@ -313,6 +318,11 @@ def evaluate(assim_time, obs_seq_out=False,
print("prepare nature")
prepare_nature_dart(assim_time) # link WRF files to DART directory
+ if posterior:
+ write_list_of_inputfiles_posterior(time)
+ if prior:
+ write_list_of_inputfiles_prior()
+
if obs_seq_out:
copy(obs_seq_out, cluster.dart_rundir+'/obs_seq.out')
else:
@@ -433,7 +443,7 @@ def link_DART_binaries_and_RTTOV_files():
print('prepared DART & RTTOV links in', cluster.dart_rundir)
except Exception as e:
- if any([hasattr(obscfg, 'sat_channel') for obscfg in exp.observations]):
+ if any(['sat_channel' in obscfg for obscfg in exp.observations]):
raise e
else:
pass # we dont need RTTOV anyway
@@ -489,20 +499,21 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
# is any observation error parametrized?
if any(error_is_parametrized) or do_QC:
- print(" (optional) evaluate prior for all observations (incl rejected)")
- evaluate(time, output_format="%Y-%m-%d_%H:%M_obs_seq.final-eval_prior_allobs")
+ print(" (optional) evaluate prior for all observations (incl rejected) ")
+ evaluate(time, prior=True,
+ output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate_prior")
print(" assign observation-errors for assimilation ")
set_obserr_assimilate_in_obsseqout(oso, outfile=cluster.dart_rundir + "/obs_seq.out")
if do_QC:
- print(" 2.3) reject observations? ")
+ print(" reject observations? ")
qc_obs(time, oso)
if prior_inflation_type == '2':
prepare_inflation_2(time, prior_init_time)
- print(" 3) run filter ")
+ print(" run filter ")
dart_nml.write_namelist()
filter(nproc=nproc)
archive_filteroutput(time)
@@ -510,12 +521,16 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
if prior_inflation_type == '2':
archive_inflation_2(time)
+ print(" evaluate posterior in observation-space")
if do_QC:
- print(" 4) evaluate posterior observations for all observations (incl rejected)")
- write_list_of_inputfiles_posterior(time)
-
- evaluate(time, obs_seq_out=cluster.archivedir+'/obs_seq_out/'+time.strftime('%Y-%m-%d_%H:%M_obs_seq.out-beforeQC'),
- output_format="%Y-%m-%d_%H:%M_obs_seq.final-eval_posterior_allobs")
+ f_oso = '%Y-%m-%d_%H:%M_obs_seq.out-beforeQC' # includes all observations (including rejected ones in qc_obs())
+ else:
+ f_oso = '%Y-%m-%d_%H:%M_obs_seq.out'
+ # evaluate() separately after ./filter is crucial when assimilating cloud variables
+ # as the obs_seq.final returned by ./filter was produced using un-clamped cloud variables
+ evaluate(time, obs_seq_out=cluster.archivedir+'/obs_seq_out/'+time.strftime(f_oso),
+ prior=False, posterior=True,
+ output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate_posterior")
if __name__ == "__main__":
diff --git a/dartwrf/dart_nml.py b/dartwrf/dart_nml.py
index 84a2412ea0f70675fe4fac77edc522152fc67c84..85716b640511a9f48c2f04aa9d33de66e9830379 100644
--- a/dartwrf/dart_nml.py
+++ b/dartwrf/dart_nml.py
@@ -1,5 +1,5 @@
+import warnings
from dartwrf.utils import append_file
-
from dartwrf.exp_config import exp
from dartwrf.server_config import cluster
@@ -191,36 +191,36 @@ def _get_list_of_localizations():
loc_horiz_rad = to_radian_horizontal(loc_horiz_km)
l_loc_horiz_rad.append(loc_horiz_rad)
- # compute vertical localization
-
- # do we have vertical localization?
- if not hasattr(obscfg, "loc_vert_km") and not hasattr(obscfg, "loc_vert_scaleheight"):
- l_loc_vert_km.append(-1)
- l_loc_vert_scaleheight.append(-1)
- # if not add dummy value
-
- # choose either localization by height or by scale height
- if hasattr(obscfg, "loc_vert_km") and hasattr(obscfg, "loc_vert_scaleheight"):
- raise ValueError("Observation config contains both loc_vert_km and loc_vert_scaleheight. Please choose one.")
-
- elif hasattr(obscfg, "loc_vert_km"): # localization by height
+ try: # localization by height
loc_vert_km = obscfg["loc_vert_km"]
vert_norm_hgt = to_vertical_normalization(loc_vert_km, loc_horiz_km)
l_loc_vert_km.append(vert_norm_hgt)
- elif hasattr(obscfg, "loc_vert_scaleheight"): # localization by scale height
- loc_vert_scaleheight = obscfg["loc_vert_scaleheight"]
+ # set the other (unused) list to a dummy value
+ l_loc_vert_scaleheight = [-1,]
- # no conversion necessary, take the values as defined in obscfg
- l_loc_vert_scaleheight.append(loc_vert_scaleheight)
+ except AttributeError: # localization by scale height
+ try:
+ loc_vert_scaleheight = obscfg["loc_vert_scaleheight"]
- # set the other (unused) list to a dummy value
- if len(l_loc_vert_km) > 0:
- l_loc_vert_scaleheight = [-1,]
- else:
- l_loc_vert_km = [-1,]
-
+ # no conversion necessary, take the values as defined in obscfg
+ l_loc_vert_scaleheight.append(loc_vert_scaleheight)
+
+ # set the other (unused) list to a dummy value
+ l_loc_vert_km = [-1,]
+
+ except AttributeError:
+
+ # do we have vertical localization?
+ # check parameter horiz_dist_only == true
+ if exp.dart_nml['&location_nml']['horiz_dist_only'][0] == '.true.':
+ # no vertical localization
+ l_loc_vert_km.append(-1)
+ l_loc_vert_scaleheight.append(-1)
+ else:
+ raise ValueError('Neither `loc_vert_km` nor `loc_vert_scaleheight` defined in obscfg.')
+
return l_obstypes, l_loc_horiz_rad, l_loc_vert_km, l_loc_vert_scaleheight
@@ -280,15 +280,6 @@ def write_namelist(just_prior_values=False):
list_obstypes, list_loc_horiz_rad, list_loc_vert_km, list_loc_vert_scaleheight = _get_list_of_localizations()
nml = read_namelist(cluster.dart_srcdir + "/input.nml")
-
- # dont compute posterior, just evaluate prior
- if just_prior_values:
- nml['&filter_nml']['compute_posterior'] = [['.false.']]
- nml['&filter_nml']['output_members'] = [['.false.']]
- nml['&filter_nml']['output_mean'] = [['.false.']]
- nml['&filter_nml']['output_sd'] = [['.false.']]
- nml['&obs_kind_nml']['assimilate_these_obs_types'] = [[]]
- nml['&obs_kind_nml']['evaluate_these_obs_types'] = [list_obstypes]
if len(list_obstypes) > 0:
# make sure that observations defined in `exp.observations` are assimilated
@@ -305,6 +296,15 @@ def write_namelist(just_prior_values=False):
nml['&location_nml']['special_vert_normalization_levels'] = [[-1,]]
nml['&location_nml']['special_vert_normalization_pressures'] = [[-1,]]
+ # dont compute posterior, just evaluate prior
+ if just_prior_values:
+ nml['&filter_nml']['compute_posterior'] = [['.false.']]
+ nml['&filter_nml']['output_members'] = [['.false.']]
+ nml['&filter_nml']['output_mean'] = [['.false.']]
+ nml['&filter_nml']['output_sd'] = [['.false.']]
+ nml['&obs_kind_nml']['assimilate_these_obs_types'] = [[]]
+ nml['&obs_kind_nml']['evaluate_these_obs_types'] = [list_obstypes]
+
# overwrite namelist parameters as defined in the experiment configuration
for section, sdata in exp.dart_nml.items():
@@ -330,8 +330,8 @@ def write_namelist(just_prior_values=False):
# fail if horiz_dist_only == false but observations contain a satellite channel
if nml['&location_nml']['horiz_dist_only'][0] == '.false.':
for obscfg in exp.observations:
- if hasattr(obscfg, "sat_channel"):
- raise ValueError("Selected vertical localization, but observations contain satellite obs -> Not possible.")
+ if 'sat_channel' in obscfg:
+ warnings.warn("Selected vertical localization, but observations contain satellite obs -> Bug in DART.")
# write to file
write_namelist_from_dict(nml, cluster.dart_rundir + "/input.nml")
diff --git a/dartwrf/evaluate_obs_space.py b/dartwrf/evaluate_obs_space.py
index 0436d926e3e075bbafb91b70a4558052c6e05f86..27269103e251934b50eaf87005c7878f2e049a53 100755
--- a/dartwrf/evaluate_obs_space.py
+++ b/dartwrf/evaluate_obs_space.py
@@ -1,49 +1,59 @@
import os, sys, shutil, warnings
-import time as time_module
import datetime as dt
-import numpy as np
from dartwrf.exp_config import exp
from dartwrf.server_config import cluster
from dartwrf import assim_synth_obs as aso
-def get_previous_obsseq_file(time):
- oso_input = cluster.archivedir+'/obs_seq_out' + time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out-beforeQC")
+def evaluate_one_time(init, valid):
+ """Evaluate the ensemble forecast in observation space at a given time, apart from the analysis time.
+
+ Args:
+ init (datetime): initialization time of the forecast
+ valid (datetime): valid time of the forecast
+
+ Returns:
+ None
+ """
+
+ # # prepare nature and prior ensemble
+ aso.prepare_nature_dart(valid)
+ aso.prepare_prior_ensemble(valid, prior_init_time=init, prior_valid_time=valid, prior_path_exp=cluster.archivedir)
+
- if not os.path.isfile(oso_input): # fallback
- oso_input = cluster.archivedir+'/obs_seq_out' + time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out")
+ if use_other_obsseq: # use a different obsseq file
+ f_obs_seq_out = use_other_obsseq
+ shutil.copy(f_obs_seq_out, cluster.dart_rundir+'/obs_seq.out')
+
+ else: # from same exp
+ # use the last assimilation obsseq file for the observation locations (note: observed values are not valid)
+ f_obs_seq_out = cluster.archivedir+valid.strftime("/obs_seq_out/%Y-%m-%d_%H:%M_obs_seq.out")
- return oso_input
+ from dartwrf.obs import obsseq
+ oso = obsseq.ObsSeq(f_obs_seq_out)
+ oso.df['observations'] = -9999
+ oso.to_dart(cluster.dart_rundir+'/obs_seq.out')
+ aso.evaluate(valid, output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate")
if __name__ == "__main__":
"""Evaluate the ensemble forecast in observation space at a given time, apart from the analysis time.
Note: Observations are not assimilated. This is only for evaluation purposes.
- """
- init = dt.datetime.strptime(sys.argv[1], "%Y-%m-%d_%H:%M")
- time = dt.datetime.strptime(sys.argv[2], "%Y-%m-%d_%H:%M")
+ Usage: python3 evaluate_obs_space.py init1,valid1 init2,valid2 ...
+ """
+ args = sys.argv[1:]
+ init_valid_tuples = [a.split(',') for a in args]
use_other_obsseq = False
# we need an existing run_DART folder
aso.prepare_run_DART_folder()
- # # prepare nature and prior ensemble
- aso.prepare_nature_dart(time)
- aso.prepare_prior_ensemble(time, prior_init_time=init, prior_valid_time=time, prior_path_exp=cluster.archivedir)
-
- # tell DART to use the prior as input
- aso.write_list_of_inputfiles_prior()
-
- if use_other_obsseq: # use a different obsseq file
- oso_input = use_other_obsseq
- else: # from same exp
-
- # use the last assimilation obsseq file for the observation locations (note: observed values are not valid)
- oso_input = get_previous_obsseq_file(init)
- shutil.copy(oso_input, cluster.dart_rundir+'/obs_seq.out')
+ for (init, valid) in init_valid_tuples:
- aso.evaluate(time, output_format="%Y-%m-%d_%H:%M_obs_seq.final-evaluate")
\ No newline at end of file
+ init = dt.datetime.strptime(init, "%Y-%m-%d_%H:%M")
+ valid = dt.datetime.strptime(valid, "%Y-%m-%d_%H:%M")
+ evaluate_one_time(init, valid)
\ No newline at end of file
diff --git a/dartwrf/obs/create_obsseq_in.py b/dartwrf/obs/create_obsseq_in.py
index ebcbed58a36b22a5351e3603308739dfa9618f53..ca990517a3bd7eb52e5e9b7c28240a509af508c3 100755
--- a/dartwrf/obs/create_obsseq_in.py
+++ b/dartwrf/obs/create_obsseq_in.py
@@ -42,7 +42,7 @@ def round_to_day(dtobj):
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)
@@ -61,14 +61,7 @@ def get_dart_date(time_dt):
return dart_date_day, secs_thatday
-def write_tuple_to_pickle(fpath_out, tuple):
- import pickle
- os.makedirs(os.path.dirname(fpath_out), exist_ok=True)
- with open(fpath_out, 'wb') as f:
- pickle.dump(tuple, f)
- print(fpath_out, 'saved.')
-
-def write_file(msg, output_path='./'):
+def _write_file(msg, output_path='./'):
try:
os.remove(output_path)
except OSError:
@@ -79,7 +72,7 @@ def write_file(msg, output_path='./'):
print(output_path, 'saved.')
-def append_hgt_to_coords(coords, heights):
+def _append_hgt_to_coords(coords, heights):
coords2 = []
for i in range(len(coords)):
for hgt_m in heights:
@@ -108,7 +101,7 @@ obs_kind_definitions
first: 1 last: """+n_obs_str
-def determine_vert_coords(sat_channel, kind, obscfg):
+def _determine_vert_coords(sat_channel, kind, obscfg):
if not sat_channel:
if 'SYNOP' in kind:
vert_coord_sys = "-1" # meters AGL
@@ -124,6 +117,15 @@ def determine_vert_coords(sat_channel, kind, obscfg):
def write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt):
"""Writes metadata str for an observation inside obs_seq.out
+
+ Args:
+ sat_channel (int or False): False if not a satellite observation
+ lat0 (float): latitude of point on earth
+ lon0 (float): longitude of point on earth
+ time_dt (dt.datetime): time of observation
+
+ Returns:
+ str
"""
if sat_channel:
sun_az = str(get_azimuth(lat0, lon0, time_dt))
@@ -146,6 +148,9 @@ def write_section(obs, last=False):
Args:
obs (object)
last (bool): True if this is the last observation in the obs_seq file
+
+ Returns:
+ str
"""
lon_rad = str(degr_to_rad(obs['lon']))
lat_rad = str(degr_to_rad(obs['lat']))
@@ -188,7 +193,7 @@ def create_obs_seq_in(time_dt, list_obscfg,
os.makedirs(os.path.dirname(output_path), exist_ok=True)
print('creating obs_seq.in:')
- time_dt = add_timezone_UTC(time_dt)
+ time_dt = _add_timezone_UTC(time_dt)
dart_date_day, secs_thatday = get_dart_date(time_dt)
txt = ''
@@ -216,8 +221,8 @@ def create_obs_seq_in(time_dt, list_obscfg,
sat_channel = obscfg.get('sat_channel', False)
# add observation locations in the vertical at different levels
- vert_coord_sys, vert_coords = determine_vert_coords(sat_channel, kind, obscfg)
- coords = append_hgt_to_coords(coords, vert_coords)
+ vert_coord_sys, vert_coords = _determine_vert_coords(sat_channel, kind, obscfg)
+ coords = _append_hgt_to_coords(coords, vert_coords)
n_obs_3d_thistype = len(coords)
# user defined generation error
@@ -252,57 +257,19 @@ def create_obs_seq_in(time_dt, list_obscfg,
pretxt = preamble(n_obs_total, list_kinds)
alltxt = pretxt + txt
- write_file(alltxt, output_path=output_path)
+ _write_file(alltxt, output_path=output_path)
if __name__ == '__main__':
# for testing
time_dt = dt.datetime(2008, 7, 30, 13, 0)
-
- vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
- kind='MSG_4_SEVIRI_BDRF', sat_channel=1,
- n_obs=961, obs_locations='square_array_evenly_on_grid',
- error_generate=0.03, error_assimilate=0.06,
- cov_loc_radius_km=32)
-
- # wv62 = dict(plotname='Brightness temperature WV 6.2µm', plotunits='[K]',
- # kind='MSG_4_SEVIRI_TB', sat_channel=5,
- # n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
- # error_generate=1., error_assimilate=False,
- # cov_loc_radius_km=20)
-
- # wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
- # kind='MSG_4_SEVIRI_TB', sat_channel=6,
- # n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
- # error_generate=1., error_assimilate=False,
- # cov_loc_radius_km=20)
-
- # ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
- # kind='MSG_4_SEVIRI_TB', sat_channel=9,
- # n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
- # error_generate=5., error_assimilate=10.,
- # cov_loc_radius_km=32)
-
- radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
- kind='RADAR_REFLECTIVITY',
- n_obs=1, obs_locations=[(45,0),],
- error_generate=2.5, error_assimilate=5.,
- heights=np.arange(1000, 15001, 1000),
- cov_loc_radius_km=20, cov_loc_vert_km=4)
-
- # t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
- # kind='SYNOP_TEMPERATURE',
- # n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
- # error_generate=0.1, error_assimilate=1.,
- # cov_loc_radius_km=20, cov_loc_vert_km=3)
-
- # psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
- # kind='SYNOP_SURFACE_PRESSURE',
- # n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
- # error_generate=50., error_assimilate=100.,
- # cov_loc_radius_km=32, cov_loc_vert_km=5)
-
+ radar = dict(var_name='Radar reflectivity', unit='[dBz]',
+ kind='RADAR_REFLECTIVITY',
+ n_obs=5776, obs_locations='square_array_evenly_on_grid',
+ error_generate=2.5, error_assimilate=2.5,
+ heights=range(2000, 14001, 2000),
+ loc_horiz_km=20, loc_vert_km=2.5)
create_obs_seq_in(time_dt, [radar],
output_path=utils.userhome+'/run_DART/obs_seq.in')
diff --git a/dartwrf/obs/create_obsseq_out.py b/dartwrf/obs/create_obsseq_out.py
index 69e95e66fdbe6b02340e9b3f597e4559b1fdde3c..af9eaa503b05dc6803879784d1db802fc5970ff9 100644
--- a/dartwrf/obs/create_obsseq_out.py
+++ b/dartwrf/obs/create_obsseq_out.py
@@ -63,14 +63,16 @@ def generate_obsseq_out(time):
print('saved', f_oso)
return oso
-def run_perfect_model_obs(nproc=12, verbose=True):
- """Run the perfect_model_obs program to generate observations
+def run_perfect_model_obs(nproc=12):
+ """Run the ./perfect_model_obs program
+
+ Args:
+ nproc (int): number of processors to use
Returns:
- None
+ None, creates obs_seq.out
"""
- if verbose:
- print("generating observations - running ./perfect_model_obs")
+ print("running ./perfect_model_obs")
os.chdir(cluster.dart_rundir)
try_remove(cluster.dart_rundir + "/obs_seq.out")
@@ -90,6 +92,9 @@ if __name__ == '__main__':
Args:
times (str): comma-separated list of times of the observations
+
+ Returns:
+ None, creates obs_seq.out in cluster.archivedir
"""
import argparse
import datetime as dt
diff --git a/dartwrf/obs/error_models.py b/dartwrf/obs/error_models.py
index 08c3236fbf24c8bb85fdd76d379adff6f964d51d..9ec0eae7ee641a267205feba6cb070aa9b369f91 100644
--- a/dartwrf/obs/error_models.py
+++ b/dartwrf/obs/error_models.py
@@ -3,15 +3,17 @@ from scipy.interpolate import interp1d
def calc_obserr_WV(channel, Hx_nature, Hx_prior):
"""Calculate parametrized error (for assimilation)
+
Args:
- channel (str): satellite channel
- Hx_nature (np.array): dim (observations)
- Hx_prior (np.array): dim (ensemble_members, observations)
+ channel (str): satellite channel
+ Hx_nature (np.array): H(x_nature) with dimension (observations)
+ Hx_prior (np.array): H(x_prior) with dimension (ensemble_members, observations)
Returns
- np.array Observation error std-deviation with dim (observations)
+ np.array Observation error std-deviation with dimension (observations)
"""
- assert channel in ['WV62', 'WV73']
+ if channel not in ['WV62', 'WV73']:
+ raise NotImplementedError("channel not implemented: " + channel)
debug = False
n_obs = len(Hx_nature)
@@ -21,13 +23,13 @@ def calc_obserr_WV(channel, Hx_nature, Hx_prior):
bt_x_ens = Hx_prior[:, iobs]
# compute Cloud impact for every pair (ensemble, truth)
- CIs = [cloudimpact(channel, bt_x, bt_y) for bt_x in bt_x_ens]
+ CIs = [_cloudimpact(channel, bt_x, bt_y) for bt_x in bt_x_ens]
mean_CI = np.mean(CIs)
if channel == 'WV62':
- oe_model = OE_model_harnisch_WV62(mean_CI)
+ oe_model = _OE_model_harnisch_WV62(mean_CI)
elif channel == 'WV73':
- oe_model = OE_model_harnisch_WV73(mean_CI)
+ oe_model = _OE_model_harnisch_WV73(mean_CI)
if debug:
print("BT_nature=", bt_y, "=> mean_CI=", mean_CI, "=> OE_assim=", oe_model)
@@ -36,7 +38,7 @@ def calc_obserr_WV(channel, Hx_nature, Hx_prior):
return OEs
-def cloudimpact(channel, bt_mod, bt_obs):
+def _cloudimpact(channel, bt_mod, bt_obs):
"""
follows Harnisch 2016, Figure 3
"""
@@ -53,7 +55,7 @@ def cloudimpact(channel, bt_mod, bt_obs):
return ci
-def OE_model_harnisch_WV62(ci):
+def _OE_model_harnisch_WV62(ci):
if ci >= 0 and ci < 7.5:
# Kelvin, fit of Fig 7a, Harnisch 2016
x_ci = [0, 2.5, 4.5, 5.5, 7.5] # average cloud impact [K]
@@ -64,7 +66,7 @@ def OE_model_harnisch_WV62(ci):
return 6.5
-def OE_model_harnisch_WV73(ci):
+def _OE_model_harnisch_WV73(ci):
if ci >= 0 and ci < 16:
# Kelvin, fit of Fig 7b, Harnisch 2016
x_ci = [0, 5, 10.5, 13, 16] # average cloud impact [K]
diff --git a/dartwrf/obs/obsseq.py b/dartwrf/obs/obsseq.py
index cbfa62213a43d63d74464593938079eff07c760e..49f3bd80498f2a3888014166b294af31cb187d44 100755
--- a/dartwrf/obs/obsseq.py
+++ b/dartwrf/obs/obsseq.py
@@ -36,7 +36,6 @@ import os, warnings
import numpy as np
import pandas as pd
-
def _plot_box(m, lat, lon, label="", **kwargs):
""""Draw bounding box
diff --git a/dartwrf/obs/obsseq_2dim.py b/dartwrf/obs/obsseq_2dim.py
index 7239906206672579edf2320578784032fc0bec3a..40afa4a7c51999d5a518d73fc3405d6d16cf2be0 100755
--- a/dartwrf/obs/obsseq_2dim.py
+++ b/dartwrf/obs/obsseq_2dim.py
@@ -14,18 +14,23 @@ Example:
"""
from copy import copy
-import os, sys, shutil, warnings
-import time as time_module
+import os, sys, warnings
import datetime as dt
import numpy as np
from dartwrf.server_config import cluster
from dartwrf import utils
-from dartwrf import assim_synth_obs as aso
from dartwrf.obs import obsseq
def _get_n_obs_per_layer(oso):
- """Get number of observations per layer"""
+ """Determine number of observations per layer from obsseq.ObsSeq object
+
+ Args:
+ oso (obsseq.ObsSeq): obsseq object
+
+ Returns:
+ int
+ """
height_all = np.array([a[2] for a in oso.df.loc3d])
height_first = height_all[0]
diff --git a/dartwrf/workflows.py b/dartwrf/workflows.py
index f7914d2a985b0c9bac4e7df7c585c0c95e1df6c5..9122a18278e6bb0dc79788d35066b0a3b58a497f 100644
--- a/dartwrf/workflows.py
+++ b/dartwrf/workflows.py
@@ -117,7 +117,8 @@ class WorkFlows(object):
_copy_dartwrf_to_archive() # includes config files
# we set the path from where python should import dartwrf modules
- self.cluster.python = 'export PYTHONPATH='+self.cluster.scripts_rundir+'; '+self.cluster.python
+ # pythonpath is the folder above the dartwrf folder
+ self.cluster.python = 'export PYTHONPATH='+self.cluster.scripts_rundir+'/../; '+self.cluster.python
# Set paths and backup scripts
self.cluster.log_dir = self.cluster.archivedir+'/logs/'
@@ -322,7 +323,26 @@ class WorkFlows(object):
id = self.cluster.run_job("obsseq_netcdf", cfg_update={"time": "10", "mail-type": "FAIL,END"},
depends_on=[depends_on])
return id
+
+ def evaluate_plus1(self, list_assim_times, depends_on=None):
+ list_of_tuples = [(init, (init+dt.timedelta(minutes=1))) for init in list_assim_times]
+ arg = ' '.join([ttuple[0].strftime('%Y-%m-%d_%H:%M,')+ttuple[1].strftime('%Y-%m-%d_%H:%M') for ttuple in list_of_tuples])
+
+ cmd = self.cluster.python+' '+self.cluster.scripts_rundir+'/evaluate_obs_space.py '+arg
+ id = self.cluster.run_job(cmd, 'eval+1'+self.exp.expname, cfg_update={"ntasks": "12", "mem": "50G", "ntasks-per-node": "12", "ntasks-per-core": "2",
+ "time": "15", "mail-type": "FAIL"},
+ depends_on=[depends_on])
+ return id
+
+ def evaluate_plus0(self, list_assim_times, depends_on=None):
+ list_of_tuples = [(init, init) for init in list_assim_times]
+ arg = ' '.join([ttuple[0].strftime('%Y-%m-%d_%H:%M,')+ttuple[1].strftime('%Y-%m-%d_%H:%M') for ttuple in list_of_tuples])
+ cmd = self.cluster.python+' '+self.cluster.scripts_rundir+'/evaluate_obs_space.py '+arg
+ id = self.cluster.run_job(cmd, 'eval+0'+self.exp.expname, cfg_update={"ntasks": "12", "mem": "50G", "ntasks-per-node": "12", "ntasks-per-core": "2",
+ "time": "15", "mail-type": "FAIL"},
+ depends_on=[depends_on])
+ return id
def verify_sat(self, depends_on=None):
cmd = self.cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+self.exp.expname+' has_node sat verif1d FSS BS'
diff --git a/docs/source/dartwrf.obs.rst b/docs/source/dartwrf.obs.rst
index 79a1e0e57472de8c7562aacda2720bb20ffb10f3..6106c161cb073f643cf7a5b703934c589d9f841e 100644
--- a/docs/source/dartwrf.obs.rst
+++ b/docs/source/dartwrf.obs.rst
@@ -12,6 +12,22 @@ dartwrf.obs.calculate\_obs\_locations module
:undoc-members:
:show-inheritance:
+dartwrf.obs.create\_obsseq\_in module
+-------------------------------------
+
+.. automodule:: dartwrf.obs.create_obsseq_in
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+dartwrf.obs.create\_obsseq\_out module
+--------------------------------------
+
+.. automodule:: dartwrf.obs.create_obsseq_out
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
dartwrf.obs.error\_models module
--------------------------------
@@ -20,6 +36,46 @@ dartwrf.obs.error\_models module
:undoc-members:
:show-inheritance:
+dartwrf.obs.obskind module
+--------------------------
+
+.. automodule:: dartwrf.obs.obskind
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+dartwrf.obs.obsseq module
+-------------------------
+
+.. automodule:: dartwrf.obs.obsseq
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+dartwrf.obs.obsseq\_2dim module
+-------------------------------
+
+.. automodule:: dartwrf.obs.obsseq_2dim
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+dartwrf.obs.obsseq\_to\_netcdf module
+-------------------------------------
+
+.. automodule:: dartwrf.obs.obsseq_to_netcdf
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+dartwrf.obs.run\_obs\_diag module
+---------------------------------
+
+.. automodule:: dartwrf.obs.run_obs_diag
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
Module contents
---------------
diff --git a/docs/source/dartwrf.rst b/docs/source/dartwrf.rst
index d8bab9e076cbad9ecc730a1b560fc0105e690691..9d0143c46d34ab052303f4ca8160869ca0e72ec1 100644
--- a/docs/source/dartwrf.rst
+++ b/docs/source/dartwrf.rst
@@ -28,14 +28,6 @@ dartwrf.create\_obs\_upfront module
:undoc-members:
:show-inheritance:
-dartwrf.create\_obsseq module
------------------------------
-
-.. automodule:: dartwrf.create_obsseq
- :members:
- :undoc-members:
- :show-inheritance:
-
dartwrf.create\_wbubble\_wrfinput module
----------------------------------------
@@ -84,38 +76,6 @@ dartwrf.exp\_config module
:undoc-members:
:show-inheritance:
-dartwrf.obskind module
-----------------------
-
-.. automodule:: dartwrf.obskind
- :members:
- :undoc-members:
- :show-inheritance:
-
-dartwrf.obsseq module
----------------------
-
-.. automodule:: dartwrf.obsseq
- :members:
- :undoc-members:
- :show-inheritance:
-
-dartwrf.obsseq\_2dim module
----------------------------
-
-.. automodule:: dartwrf.obsseq_2dim
- :members:
- :undoc-members:
- :show-inheritance:
-
-dartwrf.obsseq\_to\_netcdf module
----------------------------------
-
-.. automodule:: dartwrf.obsseq_to_netcdf
- :members:
- :undoc-members:
- :show-inheritance:
-
dartwrf.prep\_IC\_prior module
------------------------------
@@ -140,14 +100,6 @@ dartwrf.prepare\_wrfrundir module
:undoc-members:
:show-inheritance:
-dartwrf.run\_obs\_diag module
------------------------------
-
-.. automodule:: dartwrf.run_obs_diag
- :members:
- :undoc-members:
- :show-inheritance:
-
dartwrf.server\_config module
-----------------------------
diff --git a/evaluate_plus1.py b/evaluate_plus1.py
new file mode 100755
index 0000000000000000000000000000000000000000..b128c3b1dc44866356b1f7d0d4a814784f306034
--- /dev/null
+++ b/evaluate_plus1.py
@@ -0,0 +1,23 @@
+#!/usr/bin/python3
+"""
+Evaluate the forecast in observation space one minute after the analysis time.
+"""
+import datetime as dt
+from dartwrf.workflows import WorkFlows
+
+w = WorkFlows(exp_config='exp_nonlin.py', server_config='jet.py')
+
+
+assim_times = [dt.datetime(2008,7,30,12),
+ dt.datetime(2008,7,30,12,30),
+ dt.datetime(2008,7,30,13),
+ dt.datetime(2008,7,30,13,30),
+ dt.datetime(2008,7,30,14),]
+
+# generate observations at +1 minute after the assimilation time
+obs_times = [each+dt.timedelta(minutes=1) for each in assim_times]
+w.generate_obsseq_out(obs_times)
+
+
+# evaluate the forecast at +1 minute after the assimilation time
+w.evaluate_plus1(assim_times)
\ No newline at end of file
diff --git a/tests/test_inputnml.py b/tests/test_inputnml.py
index 6117ef08c0bafd7cd837036c3e1c513be4c3f868..78bf73170855de7af78fab3093dcf3e95e890104 100644
--- a/tests/test_inputnml.py
+++ b/tests/test_inputnml.py
@@ -53,5 +53,13 @@ def test_input_nml():
os.remove(test_output)
+def test_get_list_of_localizations():
+
+ output = dart_nml._get_list_of_localizations()
+ assert (['SYNOP_TEMPERATURE'], [0.0015698587127158557], [1274000.0], [-1]) == output
+
+
if __name__ == '__main__':
test_input_nml()
+
+ test_get_list_of_localizations()
\ No newline at end of file