Merge pull request #149 from kenkehoe/dl_update

Doppler Lidar Correction updates along with corrections to wind barb plotting.

act/corrections/doppler_lidar.py

@@ -1,11 +1,13 @@
 import numpy as np
 
 
-def correct_dl(obj, var_name='attenuated_backscatter', fill_value=1e-7):
+def correct_dl(obj, var_name='attenuated_backscatter', range_normalize=True,
+               fill_value=1e-7):
     """
     This procedure corrects doppler lidar data by filling all zero and
     negative values of backscatter with fill_value and then converting
-    the backscatter data into logarithmic space.
+    the backscatter data into logarithmic space. Will also range normalize
+    the values by multipling by the range^2 with default units.
 
     Parameters
     ----------
@@ -14,7 +16,9 @@ def correct_dl(obj, var_name='attenuated_backscatter', fill_value=1e-7):
         in linear space.
     var_name: str
         The variable name of data in the object.
-    fill_value: float
+    range_normalize : bool
+        Option to range normalize the data.
+    fill_value : float
         The fill_value to use. The fill_value is entered in linear space.
 
     Returns
@@ -25,19 +29,23 @@ def correct_dl(obj, var_name='attenuated_backscatter', fill_value=1e-7):
     """
     backscat = obj[var_name].values
 
-    # This will get the name of the coordinate dimension so it's not assumed
-    # via position or name
-    height_name = list(set(obj[var_name].dims) - set(['time']))[0]
-    height = obj[height_name].values
-    height = height ** 2
+    if range_normalize:
+        # This will get the name of the coordinate dimension so it's not assumed
+        # via position or name
+        height_name = list(set(obj[var_name].dims) - set(['time']))[0]
+        height = obj[height_name].values
+
+        height = height / np.max(height)
+
+        backscat = backscat * height ** 2
 
-    # Doing this trick with height to change the array shape so it
-    # will broadcast correclty against backscat
-    backscat = backscat * height[None, :]
     backscat[backscat <= 0] = fill_value
     obj[var_name].values = np.log10(backscat)
 
     # Updating the units to correctly indicate the values are log values
-    obj[var_name].attrs['units'] = 'log(' + obj[var_name].attrs['units'] + ')'
+    if range_normalize:
+        obj[var_name].attrs['units'] = 'log( Range normalized ' + obj[var_name].attrs['units'] + ')'
+    else:
+        obj[var_name].attrs['units'] = 'log(' + obj[var_name].attrs['units'] + ')'
 
     return obj

act/plotting/TimeSeriesDisplay.py

@@ -651,25 +651,25 @@ def plot_barbs_from_u_v(self, u_field, v_field, pres_field=None,
 
         else:
             if 'cmap' in kwargs.keys():
-                map_color = np.sqrt(np.power(u[::barb_step_y, ::barb_step_x], 2) +
-                                    np.power(v[::barb_step_y, ::barb_step_x], 2))
+                map_color = np.sqrt(np.power(u[::barb_step_x, ::barb_step_y], 2) +
+                                    np.power(v[::barb_step_x, ::barb_step_y], 2))
                 map_color[np.isnan(map_color)] = 0
 
                 ax = self.axes[subplot_index].barbs(
-                    xdata[::barb_step_y, ::barb_step_x],
-                    ydata[::barb_step_y, ::barb_step_x],
-                    u[::barb_step_y, ::barb_step_x],
-                    v[::barb_step_y, ::barb_step_x], map_color,
+                    xdata[::barb_step_x, ::barb_step_y],
+                    ydata[::barb_step_x, ::barb_step_y],
+                    u[::barb_step_x, ::barb_step_y],
+                    v[::barb_step_x, ::barb_step_y], map_color,
                     **kwargs)
                 plt.colorbar(ax, label='Wind Speed (' +
                              self._arm[dsname][u_field].attrs['units'] + ')')
 
             else:
                 ax = self.axes[subplot_index].barbs(
-                    xdata[::barb_step_y, ::barb_step_x],
-                    ydata[::barb_step_y, ::barb_step_x],
-                    u[::barb_step_y, ::barb_step_x],
-                    v[::barb_step_y, ::barb_step_x],
+                    xdata[::barb_step_x, ::barb_step_y],
+                    ydata[::barb_step_x, ::barb_step_y],
+                    u[::barb_step_x, ::barb_step_y],
+                    v[::barb_step_x, ::barb_step_y],
                     **kwargs)
 
         if day_night_background is True:

act/retrievals/doppler_lidar.py

@@ -8,7 +8,7 @@ def compute_winds_from_ppi(obj, elevation_name='elevation', azimuth_name='azimut
                            snr_name='signal_to_noise_ratio',
                            intensity_name=None,
                            snr_threshold=0.008, remove_all_missing=False,
-                           condition_limit=1.0e4):
+                           condition_limit=1.0e4, return_obj=None):
     """
     This function will convert a Doppler Lidar PPI scan into vertical distribution
     of horizontal wind direction and speed.
@@ -39,6 +39,10 @@ def compute_winds_from_ppi(obj, elevation_name='elevation', azimuth_name='azimut
     condition_limit: float
         Upper limit used with Normalized data to check if data should be converted
         from scan signal to noise ration to wind speeds and directions.
+    return_obj : None or  Xarray Dataset Object
+        If set to a Xarray Dataset Object the calculated winds object will
+        be concatinated onto this object. This is to allow looping over
+        this function for many scans and returning a single object.
 
     Returns
     -------
@@ -61,8 +65,6 @@ def compute_winds_from_ppi(obj, elevation_name='elevation', azimuth_name='azimut
 
     """
 
-    return_obj = None
-
     azimuth = obj[azimuth_name].values
     azimuth_rounded = np.round(azimuth).astype(int)
 

act/tests/__init__.py

@@ -22,4 +22,4 @@
                            EXAMPLE_CEIL_WILDCARD, EXAMPLE_ANL_CSV,
                            EXAMPLE_MPL_1SAMPLE, EXAMPLE_IRT25m20s,
                            EXAMPLE_MET_CONTOUR, EXAMPLE_NAV,
-                           EXAMPLE_AOSMET)
+                           EXAMPLE_AOSMET, EXAMPLE_DLPPI)

act/tests/test_correct.py

@@ -40,3 +40,24 @@ def test_correct_wind():
 
     assert round(obj['wind_speed_corrected'].values[800]) == 5.0
     assert round(obj['wind_direction_corrected'].values[800]) == 92.0
+
+
+def test_correct_dl():
+    # Test the DL correction script on a PPI dataset eventhough it will
+    # mostlikely be used on FPT scans. Doing this to save space with only
+    # one datafile in the repo.
+    files = act.tests.sample_files.EXAMPLE_DLPPI
+    obj = act.io.armfiles.read_netcdf(files)
+
+    new_obj = act.corrections.doppler_lidar.correct_dl(obj, fill_value=np.nan)
+    data = new_obj['attenuated_backscatter'].data
+    data[np.isnan(data)] = 0.
+    data = data * 100.
+    data = data.astype(np.int64)
+    assert np.sum(data) == -18633551
+
+    new_obj = act.corrections.doppler_lidar.correct_dl(obj, range_normalize=False)
+    data = new_obj['attenuated_backscatter'].data
+    data[np.isnan(data)] = 0.
+    data = data.astype(np.int64)
+    assert np.sum(data) == -224000