diff --git a/dartwrf/obs/create_obsseq_out.py b/dartwrf/obs/create_obsseq_out.py
index 69e95e66fdbe6b02340e9b3f597e4559b1fdde3c..d01e1e869720778a37ac9f3b6ae250682b4bc2c8 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")
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 e108e70491dcf26489db1958fb24fbf4e865368c..19ee4d5e2620aa749c6102e033ddb2c25ff67519 100755
--- a/dartwrf/obs/obsseq.py
+++ b/dartwrf/obs/obsseq.py
@@ -57,9 +57,7 @@ def rad_to_degrees(rad):
class ObsRecord(pd.DataFrame):
- """Content of obsseq.ObsSeq instances
- provides additional methods for pd.DataFrame
- created inside an ObsSeq instance
+ """Basically a pd.DataFrame with additional methods
"""
@property
def _constructor(self):
@@ -68,18 +66,31 @@ class ObsRecord(pd.DataFrame):
return ObsRecord
def get_prior_Hx(self):
- """Return prior Hx array (n_obs, n_ens)"""
+ """Retrieve H(x_prior) for all ensemble members
+
+ Returns:
+ np.array (n_obs, n_ens)
+ """
return self._get_model_Hx('prior')
def get_posterior_Hx(self):
- """Return posterior Hx array (n_obs, n_ens)"""
+ """Retrieve H(x_posterior) for all ensemble members
+
+ Returns:
+ np.array (n_obs, n_ens)
+ """
return self._get_model_Hx('posterior')
def get_truth_Hx(self):
+ """Retrieve H(x_truth)
+
+ Returns:
+ np.array (n_obs,)
+ """
return self['truth'].values
def _get_model_Hx(self, what):
- """Return ensemble member info
+ """Retrieve a subset of the obs-sequence table, e.g. H(x_prior)
Args:
self (pd.DataFrame): usually self.self
@@ -106,6 +117,11 @@ class ObsRecord(pd.DataFrame):
return Hx.values
def get_lon_lat(self):
+ """Retrieve longitude and latitude of observations
+
+ Returns:
+ pd.DataFrame (n_obs, 2)
+ """
lats = np.empty(len(self), np.float32)
lons = lats.copy()
@@ -120,13 +136,13 @@ class ObsRecord(pd.DataFrame):
return pd.DataFrame(index=self.index, data=dict(lat=lats, lon=lons))
- def get_from_cartesian_grid(self, i, j, k):
+ def _get_from_cartesian_grid(self, i, j, k):
"""Get the observation using cartesian grid indices (ix, iy, iz)
"""
# find indices of observations within pandas.DataFrame
return self.iloc[self.i_obs_grid[i, j, k].ravel()]
- def determine_nlayers(self):
+ def _determine_nlayers(self):
nlayers = 1 # first guess
from dartwrf.exp_config import exp
@@ -203,7 +219,7 @@ class ObsRecord(pd.DataFrame):
print("window (#obs in x/y)=", win_obs)
# superob in case of multiple layers, only implemented for single obstype
- nlayers = self.determine_nlayers()
+ nlayers = self._determine_nlayers()
# indices of observations (starting from 0)
i_obs_grid = self.index.values
@@ -237,7 +253,7 @@ class ObsRecord(pd.DataFrame):
print("obs indices box=", i_obs_grid[i : i + win_obs, j : j + win_obs, k])
# find observations within superob window
- obs_box = self.get_from_cartesian_grid(slice(i, i + win_obs),
+ obs_box = self._get_from_cartesian_grid(slice(i, i + win_obs),
slice(j, j + win_obs),
k)
@@ -245,10 +261,10 @@ class ObsRecord(pd.DataFrame):
# save boundary of box to list, for plotting later
eps = dx_obs_lat_deg/2 # for plotting
eps2 = eps*0.8 # for plotting
- lat1, lon1 = self.get_from_cartesian_grid(i, j, k).get_lon_lat().values[0]
- lat2, lon2 = self.get_from_cartesian_grid(i+win_obs-1, j, k).get_lon_lat().values[0]
- lat3, lon3 = self.get_from_cartesian_grid(i, j+win_obs-1, k).get_lon_lat().values[0]
- lat4, lon4 = self.get_from_cartesian_grid(i+win_obs-1, j+win_obs-1, k).get_lon_lat().values[0]
+ lat1, lon1 = self._get_from_cartesian_grid(i, j, k).get_lon_lat().values[0]
+ lat2, lon2 = self._get_from_cartesian_grid(i+win_obs-1, j, k).get_lon_lat().values[0]
+ lat3, lon3 = self._get_from_cartesian_grid(i, j+win_obs-1, k).get_lon_lat().values[0]
+ lat4, lon4 = self._get_from_cartesian_grid(i+win_obs-1, j+win_obs-1, k).get_lon_lat().values[0]
boxes.append(([lat1-eps2, lat2+eps2, lat3-eps2, lat4+eps2],
[lon1-eps, lon2-eps, lon3+eps, lon4+eps]))
@@ -317,16 +333,17 @@ class ObsSeq(object):
def __init__(self, filepath):
self.ascii = open(filepath, "r").readlines()
- self.get_preamble_content()
- self.read_preamble()
+ self._get_preamble_content()
+ self._read_preamble()
self.df = ObsRecord(self.to_pandas())
def __str__(self):
return self.df.__str__()
- def get_preamble_content(self):
+ def _get_preamble_content(self):
"""Split the obs_seq.out file into two parts
+
1) First lines of obs_seq.out file until the first observation message
2) Observation contents
"""
@@ -339,7 +356,7 @@ class ObsSeq(object):
self.preamble = self.ascii[:i]
self.content = self.ascii[i:]
- def read_preamble(self):
+ def _read_preamble(self):
"""Defines
self.obstypes (tuple(kind_nr, kind_descriptor)) for each obs type
"""
@@ -394,7 +411,7 @@ class ObsSeq(object):
self.keys_for_values = keys
- def obs_to_dict(self):
+ def _obs_to_dict(self):
"""Convert an obs_seq.out file to a dictionary"""
obs_begin_str = "OBS "
@@ -495,7 +512,7 @@ class ObsSeq(object):
def to_pandas(self):
"""Create pd.DataFrame with rows=observations
"""
- obs_dict_list = self.obs_to_dict()
+ obs_dict_list = self._obs_to_dict()
# convert to pandas.DataFrame
# each observation is one line
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]