diff --git a/act/corrections/doppler_lidar.py b/act/corrections/doppler_lidar.py index 0f9660765b..c9347a9a5d 100644 --- a/act/corrections/doppler_lidar.py +++ b/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 diff --git a/act/plotting/TimeSeriesDisplay.py b/act/plotting/TimeSeriesDisplay.py index 867ef5310a..e6c9c68ff7 100644 --- a/act/plotting/TimeSeriesDisplay.py +++ b/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: diff --git a/act/retrievals/doppler_lidar.py b/act/retrievals/doppler_lidar.py index 252b24ede4..66013bcba3 100644 --- a/act/retrievals/doppler_lidar.py +++ b/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) diff --git a/act/tests/__init__.py b/act/tests/__init__.py index 3a91b595ac..7969395cdd 100644 --- a/act/tests/__init__.py +++ b/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) diff --git a/act/tests/test_correct.py b/act/tests/test_correct.py index 70799ac5bf..b22bd425fe 100644 --- a/act/tests/test_correct.py +++ b/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