ENH: Updating the sobel CBH retrieval method to be more robust for ot… (

ARM-DOE#646)

* ENH: Updating the sobel CBH retrieval method to be more robust for other instruments

* ENH: Updating test for new errors

* BUG: Bug fix for the pbl tests

* ENH: Updating for Max's comments

* ENH: Update

act/corrections/mpl.py

@@ -92,7 +92,7 @@ def correct_mpl(
         )
 
     # 1 - Remove negative height data
-    ds = ds.where(ds[height_var_name] > 0.0, drop=True)
+    ds = ds.where(ds[height_var_name].load() > 0.0, drop=True)
     height = ds[height_var_name].values
 
     # Get indices for calculating background
@@ -109,16 +109,16 @@ def correct_mpl(
 
     # Run through co and cross pol data for corrections
     co_bg = dummy[co_pol_var_name]
-    co_bg = co_bg.where(co_bg > -9998.0)
+    co_bg = co_bg.where(co_bg.load() > -9998.0)
     co_bg = co_bg.mean(dim='dim_0').values
 
     x_bg = dummy[cross_pol_var_name]
-    x_bg = x_bg.where(x_bg > -9998.0)
+    x_bg = x_bg.where(x_bg.load() > -9998.0)
     x_bg = x_bg.mean(dim='dim_0').values
 
     # Seems to be the fastest way of removing background signal at the moment
-    co_data = ds[co_pol_var_name].where(ds[co_pol_var_name] > 0).values
-    x_data = ds[cross_pol_var_name].where(ds[cross_pol_var_name] > 0).values
+    co_data = ds[co_pol_var_name].where(ds[co_pol_var_name].load() > 0).values
+    x_data = ds[cross_pol_var_name].where(ds[cross_pol_var_name].load() > 0).values
     for i in range(len(ds['time'].values)):
         co_data[i, :] = co_data[i, :] - co_bg[i]
         x_data[i, :] = x_data[i, :] - x_bg[i]

act/retrievals/cbh.py

@@ -15,6 +15,8 @@ def generic_sobel_cbh(
     var_thresh=None,
     fill_na=None,
     return_thresh=False,
+    filter_type='uniform',
+    edge_thresh=5.,
 ):
     """
     Function for calculating cloud base height from lidar/radar data
@@ -23,6 +25,9 @@ def generic_sobel_cbh(
     that there have been similar methods employed to detect boundary
     layer heights.
 
+    NOTE: The returned variable now appends the field name of the
+    data used to generate the CBH as part of the variable name.  cbh_sobel_[varname]
+
     Parameters
     ----------
     ds : ACT xarray.Dataset
@@ -34,7 +39,14 @@ def generic_sobel_cbh(
     var_thresh : float
         Thresholding for variable if needed.
     fill_na : float
-        What to fill nans with in DataArray if any.
+        Value to fill nans with in DataArray if any.
+    filter_type : string
+        Currently the only option is for uniform filtering.
+        uniform: Apply uniform filtering after the sobel filter?  Applies a standard area of 3x3 filtering
+        None: Excludes the filtering
+    edge_thresh : float
+        Threshold value for finding the edge after the sobel filtering.
+        If the signal is not strong, this may need to be lowered
 
     Returns
     -------
@@ -79,42 +91,46 @@ def generic_sobel_cbh(
         fill_na = var_thresh
 
     # Pull data into Standalone DataArray
-    data = ds[variable]
+    da = ds[variable]
 
     # Apply thresholds if set
     if var_thresh is not None:
-        data = data.where(data.values > var_thresh)
+        da = da.where(da.values > var_thresh)
 
     # Fill with fill_na values
-    data = data.fillna(fill_na)
+    da = da.fillna(fill_na)
 
     # If return_thresh is True, replace variable data with
     # thresholded data
     if return_thresh is True:
-        ds[variable].values = data.values
+        ds[variable].values = da.values
 
     # Apply Sobel filter to data and smooth the results
-    data = data.values
+    data = da.values.tolist()
     edge = ndimage.sobel(data)
-    edge = ndimage.uniform_filter(edge, size=3, mode='nearest')
+    if filter_type == 'uniform':
+        edge = ndimage.uniform_filter(edge, size=3, mode='nearest')
 
     # Create Data Array
-    edge_ds = xr.DataArray(edge, dims=ds[variable].dims)
+    edge_da = xr.DataArray(edge, dims=ds[variable].dims)
 
     # Filter some of the resulting edge data to get defined edges
-    edge_ds = edge_ds.where(edge_ds > 5.0)
-    edge_ds = edge_ds.fillna(fill_na)
+    edge_da = edge_da.where(edge_da > edge_thresh)
+    edge_da = edge_da.fillna(fill_na)
 
     # Do a diff along the height dimension to define edge
-    diff = edge_ds.diff(dim=1)
+    diff = edge_da.diff(dim=1).values
 
     # Get height variable to use for cbh
     height = ds[height_dim].values
 
     # Run through times and find the height
     cbh = []
     for i in range(np.shape(diff)[0]):
-        index = np.where(diff[i, :] > 5.0)[0]
+        try:
+            index = np.where(diff[i, :] > edge_thresh)[0]
+        except ValueError():
+            index = []
         if len(np.shape(height)) > 1:
             ht = height[i, :]
         else:
@@ -126,11 +142,12 @@ def generic_sobel_cbh(
             cbh.append(np.nan)
 
     # Create DataArray to add to the dataset
+    var_name = 'cbh_sobel_' + variable
     da = xr.DataArray(cbh, dims=['time'], coords=[ds['time'].values])
-    ds['cbh_sobel'] = da
-    ds['cbh_sobel'].attrs['long_name'] = ' '.join(
+    ds[var_name] = da
+    ds[var_name].attrs['long_name'] = ' '.join(
         ['CBH calculated from', variable, 'using sobel filter']
     )
-    ds['cbh_sobel'].attrs['units'] = ds[height_dim].attrs['units']
+    ds[var_name].attrs['units'] = ds[height_dim].attrs['units']
 
     return ds

act/tests/test_retrievals.py

@@ -63,9 +63,14 @@ def test_generic_sobel_cbh():
 
     ceil = ceil.resample(time='1min').nearest()
     ceil = act.retrievals.cbh.generic_sobel_cbh(
-        ceil, variable='backscatter', height_dim='range', var_thresh=1000.0, fill_na=0
+        ceil,
+        variable='backscatter',
+        height_dim='range',
+        var_thresh=1000.0,
+        fill_na=0.,
+        edge_thresh=5
     )
-    cbh = ceil['cbh_sobel'].values
+    cbh = ceil['cbh_sobel_backscatter'].values
     assert cbh[500] == 615.0
     assert cbh[1000] == 555.0
     ceil.close()
@@ -207,11 +212,11 @@ def test_calculate_pbl_liu_liang():
     assert ds['pblht_regime_liu_liang'].values == 'SBL'
 
     with np.testing.assert_raises(ValueError):
-        ds2 = ds.where(ds['alt'] < 1000.0, drop=True)
+        ds2 = ds.where(ds['alt'].load() < 1000.0, drop=True)
         ds2 = act.retrievals.sonde.calculate_pbl_liu_liang(ds2, smooth_height=15)
 
     with np.testing.assert_raises(ValueError):
-        ds2 = ds.where(ds['pres'] < 200.0, drop=True)
+        ds2 = ds.where(ds['pres'].load() < 200.0, drop=True)
         ds2 = act.retrievals.sonde.calculate_pbl_liu_liang(ds2, smooth_height=15)
 
     with np.testing.assert_raises(ValueError):

examples/retrievals/plot_cbh_sobel.py

@@ -0,0 +1,40 @@
+"""
+Cloud Base Height Retrievals
+----------------------------
+
+This example shows how to calculate the cloud base heights
+using the sobel edge detection method.  This can be used
+for vertical radar and lidar data.
+
+Author: Adam Theisen
+
+"""
+
+import glob
+from matplotlib import pyplot as plt
+import act
+import numpy as np
+
+# Read Ceilometer data for an example
+file = sorted(glob.glob(act.tests.sample_files.EXAMPLE_CEIL1))
+ds = act.io.armfiles.read_netcdf(file)
+
+ds = act.retrievals.cbh.generic_sobel_cbh(ds, variable='backscatter', height_dim='range',
+                                          var_thresh=1000.0, fill_na=0.)
+
+# Plot the cloud base height data
+display = act.plotting.TimeSeriesDisplay(ds, subplot_shape=(1, 2), figsize=(16, 6))
+display.plot('backscatter', subplot_index=(0, 0))
+title = 'SGP Ceilometer with Lidar-Calculated CBH Overplotted'
+display.plot('first_cbh', subplot_index=(0, 0), color='k', set_title=title)
+
+display.plot('backscatter', subplot_index=(0, 1))
+title = 'SGP Ceilometer with CBH Overplotted'
+display.plot('cbh_sobel_backscatter', color='k', subplot_index=(0, 1), set_title=title)
+
+diff = ds['first_cbh'].values - ds['cbh_sobel_backscatter'].values
+
+print("Average difference between ceilomter and sobel heights ", np.nanmean(diff))
+
+ds.close()
+plt.show()