From 7df30c01181e06d00379d2f056dd8b847263e24f Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Tue, 18 Oct 2022 11:10:28 +0200
Subject: [PATCH] other
---
dartwrf/obs/__init__.py | 0
dartwrf/obsseq.py | 21 ++++--
dartwrf/old_functions.py | 145 +++++++++++++++++++++++++++++++++++++++
dartwrf/run_obs_diag.py | 6 +-
templates/input.nml | 2 +-
5 files changed, 164 insertions(+), 10 deletions(-)
create mode 100644 dartwrf/obs/__init__.py
create mode 100644 dartwrf/old_functions.py
diff --git a/dartwrf/obs/__init__.py b/dartwrf/obs/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/dartwrf/obsseq.py b/dartwrf/obsseq.py
index cd71c54..e3fa8ad 100755
--- a/dartwrf/obsseq.py
+++ b/dartwrf/obsseq.py
@@ -77,7 +77,13 @@ class ObsRecord(pd.DataFrame):
def get_posterior_Hx(self):
"""Return posterior Hx array (n_obs, n_ens)"""
- return self._get_model_Hx('posterior')
+ try:
+ return self._get_model_Hx('posterior')
+ except Exception as e:
+ # this is useful if we evaluate a 'posterior state'
+ # then the variable is called 'prior' while it really is a posterior
+ warnings.warn(str(e)+' returning prior from this file instead!')
+ return self._get_model_Hx('prior')
def get_truth_Hx(self):
return self['truth'].values
@@ -691,12 +697,15 @@ class ObsSeq(object):
if __name__ == "__main__":
# for testing purposes
- obs = ObsSeq(cluster.scriptsdir + "/../tests/obs_seq.orig.out")
- print(type(obs))
+ # f = cluster.scriptsdir + "/../tests/obs_seq.orig.out"
+ f = "/home/fs71386/lkugler/data/sim_archive/exp_v1.21_P3_wbub7_VIS_obs2-10_loc20/obs_seq_out/2008-07-30_12:30_obs_seq.out-orig"
+ obs = ObsSeq(f)
# select a subset (lat-lon)
- obs.df = obs.df[:].superob(window_km=50)
- print(type(obs.df))
+ obs.df = obs.df.superob(window_km=10)
+ # print(type(obs.df))
+
+ obs.plot(f_out="./map_obs_superobs.png")
# write to obs_seq.out in DART format
- obs.to_dart(f=cluster.dartrundir + "/obs_seq.out")
+ # obs.to_dart(f=cluster.dartrundir + "/obs_seq.out")
diff --git a/dartwrf/old_functions.py b/dartwrf/old_functions.py
new file mode 100644
index 0000000..c8ffebe
--- /dev/null
+++ b/dartwrf/old_functions.py
@@ -0,0 +1,145 @@
+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)
+
+ print('recycle DART output to be used as input')
+ for iens in range(1, exp.n_ens+1):
+ dart_output = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
+ dart_input = cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01'
+
+ #print('check for non-monotonic vertical pressure')
+
+ # convert link to file in order to be able to update the content
+ if os.path.islink(dart_input):
+ l = os.readlink(dart_input)
+ os.remove(dart_input)
+ copy(l, dart_input)
+
+ # print('move DART output to input: '+dart_output+' -> '+dart_input)
+ # os.rename(dart_output, dart_input) # probably doesnt work
+
+ print('updating', updates, 'in', dart_input, 'from', dart_output)
+ os.system(cluster.ncks+' -A -v '+updates+' '+dart_output+' '+dart_input)
+
+
+def run_Hx(time, obscfg):
+ """
+ # assumes that prior ensemble is already linked to advance_temp<i>/wrfout_d01
+ Creates
+ obs_seq.final (file): observations on (non-averaged) grid
+ """
+ get_obsseq_out(time)
+ print("running H(x) : obs operator on ensemble prior")
+
+ # set input.nml to calculate prior Hx
+ set_DART_nml(just_prior_values=True)
+
+ # run filter to calculate prior Hx
+ shell("mpirun -np 12 ./filter &> log.filter.preassim")
+
+ osf = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.final")
+
+ if getattr(exp, "superob_km", False):
+ df = osf.df
+ print("superobbing to", exp.superob_km, "km")
+ osf.df = osf.df.superob(window_km=exp.superob_km)
+ osf.to_dart(cluster.dartrundir + "/obs_seq.final")
+ return osf
+
+
+
+def replace_errors_obsseqout(f, new_errors):
+ """Replaces existing observation errors in obs_seq.final files
+
+ 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:
+ first_obs = second_obs = None
+ for i, line in enumerate(obsseq):
+ if " OBS " in line:
+ if first_obs is None:
+ first_obs = i
+ else:
+ second_obs = i
+ break
+ if not second_obs:
+ raise RuntimeError('just one OBS in this file?! '+str(f))
+ lines_between = second_obs - first_obs
+ lines_obserr_after_obsnr = lines_between - 1 # obserr line is one before ' OBS ' line
+
+ # replace values in list obsseq
+ i_obs = 0
+ for i, line in enumerate(obsseq):
+ if " OBS " in line:
+ line_error_obs_i = i + lines_obserr_after_obsnr
+
+ previous_err_var = obsseq[line_error_obs_i]
+ new_err_obs_i = new_errors[i_obs] ** 2 # variance in obs_seq.out
+ if debug:
+ print(
+ line.strip(),
+ "previous err var ",
+ float(previous_err_var.strip()),
+ "new error",
+ new_err_obs_i,
+ )
+ obsseq[line_error_obs_i] = " " + str(new_err_obs_i) + " \n"
+
+ i_obs += 1 # next iteration
+
+ with open(f, "w") as file:
+ for line in obsseq:
+ file.write(line)
+ print("replaced obs errors in", f)
+
+
+
+def obs_operator_nature(time):
+ print("getting true values in obs space from nature run")
+ prepare_nature_dart(time)
+ run_perfect_model_obs()
+ true, _ = read_truth_obs_obsseq(cluster.dartrundir + "/obs_seq.out")
+ return true
+
+
+
+def read_truth_obs_obsseq(f):
+ """Reads observed and true values from obs_seq.out/final files."""
+ obsseq = open(f, "r").readlines()
+ true = []
+ obs = []
+ # read observations from obs_seq.out
+ for i, line in enumerate(obsseq):
+ if " OBS " in line:
+ observed = float(obsseq[i + 1].strip())
+ truth = float(obsseq[i + 2].strip())
+ true.append(truth)
+ obs.append(observed)
+ return true, obs
+
+
+def read_prior_obs(f_obs_prior):
+ """
+ docstring
+ """
+ obsseq = open(f_obs_prior, "r").readlines()
+ OBSs = []
+ # read observations from obs_seq.final
+ for i, line in enumerate(obsseq):
+ if " OBS " in line:
+ observed = float(obsseq[i + 1].strip())
+ truth = float(obsseq[i + 2].strip())
+ prior_ensmean = float(obsseq[i + 3].strip())
+ prior_enssd = float(obsseq[i + 4].strip())
+ prior_ens = []
+ for j in range(5, 5 + exp.n_ens):
+ prior_ens.append(float(obsseq[i + j].strip()))
+
+ OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens)))
+ return OBSs
\ No newline at end of file
diff --git a/dartwrf/run_obs_diag.py b/dartwrf/run_obs_diag.py
index d3d3feb..fb327b5 100644
--- a/dartwrf/run_obs_diag.py
+++ b/dartwrf/run_obs_diag.py
@@ -1,6 +1,6 @@
import os, sys, shutil, glob
from config.cfg import exp, cluster
-from utils import symlink, copy, sed_inplace, append_file
+from utils import symlink, copy, sed_inplace, append_file, shell
rundir_program = '/home/fs71386/lkugler/data/run_DART/'
@@ -38,7 +38,7 @@ def run_obsdiag(filepaths, f_out='./obsdiag.nc'):
print('------ running obs_diag program')
os.chdir(rundir_program)
symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
- os.system('./obs_diag >& obs_diag.log') # caution, this overwrites obs_seq_to_netcdf
+ shell(cluster.container, './obs_diag >& obs_diag.log') # caution, this overwrites obs_seq_to_netcdf
# move output to archive
#outdir = os.path.dirname(f_out) #'/'.join(folder_obs_seq_final.split('/')[:-1])
@@ -56,7 +56,7 @@ def run_obs_seq_to_netcdf(filepaths, f_out='./obs_epoch.nc'):
print('------ running obs_seq_to_netcdf program')
#shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', rundir_program+'/obs_seq_to_netcdf')
os.chdir(rundir_program)
- os.system('./obs_seq_to_netcdf >& obs_seq_to_netcdf.log') # caution, overwrites its own binary?!
+ shell(cluster.container, './obs_seq_to_netcdf >& obs_seq_to_netcdf.log') # caution, overwrites its own binary?!
shutil.move(rundir_program+'/obs_epoch_001.nc', f_out)
print(f_out, 'saved.')
diff --git a/templates/input.nml b/templates/input.nml
index 97968eb..0e9d7a8 100644
--- a/templates/input.nml
+++ b/templates/input.nml
@@ -103,7 +103,7 @@
# localization radius, where the influence of an obs decreases
# to ~half at 300 km, and ~0 at the edges of the area.
&assim_tools_nml
- filter_kind = 1,
+ filter_kind = <filter_kind>,
cutoff = <cov_loc_radian>,
sort_obs_inc = .false.,
spread_restoration = .false.,
--
GitLab