docs

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")

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)
+        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]

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

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]