STY: Pep8 fixes

doc/source/conf.py

@@ -13,7 +13,7 @@
 # import os
 # import sys
 # sys.path.insert(0, os.path.abspath('.'))
-
+import sphinx_theme
 
 # -- Project information -----------------------------------------------------
 
@@ -46,7 +46,7 @@
     'gallery_dirs': 'source/auto_examples'
 }
 
-#Add any paths that contain templates here, relative to this directory.
+# Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
 
 # List of patterns, relative to source directory, that match files and
@@ -67,9 +67,9 @@
 # so a file named "default.css" will overwrite the builtin "default.css".
 html_static_path = ['_static']
 html_theme = 'neo_rtd_theme'
-import sphinx_theme
+
 html_theme_path = [sphinx_theme.get_html_theme_path()]
 
 # Numpy autodoc attributes
 numpydoc_show_class_members = True
-autosummary_generate = True
+autosummary_generate = True

examples/plot_moments.py

@@ -25,4 +25,4 @@
 my_ds["doppler_velocity_max_peak"].plot(ax=ax[1], y='range')
 plt.show()
 
-test_file.close()
+test_file.close()

examples/plot_power_spectra.py

@@ -17,8 +17,8 @@
 # Plot the power spectra for a given time and height
 my_time = datetime(2017, 8, 4, 0, 40, 59)
 fig, ax = plt.subplots(1, 2, figsize=(15, 5))
-my_ds["power_spectra_normed_interp"].sel(time=my_time, range=350., method='nearest').plot(ax=ax[0])
-my_ds["power_spectra_normed_interp"].sel(time=my_time, range=950., method='nearest').plot(ax=ax[1])
+out_ds["power_spectral_density"].sel(time=my_time, range=350., method='nearest').plot(ax=ax[0])
+out_ds["power_spectral_density"].sel(time=my_time, range=950., method='nearest').plot(ax=ax[1])
 plt.show()
 
-test_file.close()
+test_file.close()

examples/read_00_file.py

@@ -1,6 +1,11 @@
+"""
+
+Example on reading raw Halo Photonics data
+------------------------------------------
+"""
 import highiq
 
 ds = highiq.io.read_00_data('sgpdlacfC1.00.20220107.141001.raw.aet_Stare_107_20220107_14.raw',
                             'sgpdlprofcalC1.home_point')
 my_spectra = highiq.calc.get_psd(ds)
-print(ds)
+print(ds)

highiq/calc/moments.py

@@ -5,7 +5,6 @@
 
 def _gpu_calc_power(psd, dV, block_size=200, normed=True):
     shp = psd.shape
-
     power = np.zeros((shp[0], shp[1]))
     if len(shp) == 3:
         gpu_array = tf.constant(psd, dtype=tf.float32)
@@ -25,7 +24,6 @@ def _gpu_calc_power(psd, dV, block_size=200, normed=True):
 
 def _gpu_calc_velocity(psd, power, vel_bins, dV):
     shp = psd.shape
-    times = shp[0]
     gpu_array = tf.constant(psd, dtype=tf.float32)
     power_array = tf.constant(power, dtype=tf.float32)
     vel_bins_tiled = np.tile(vel_bins, (shp[0], shp[1], 1))
@@ -35,7 +33,6 @@ def _gpu_calc_velocity(psd, power, vel_bins, dV):
 
 
 def _gpu_calc_velocity_dumb(psd, vel_bins):
-    shp = psd.shape
     dV = np.diff(vel_bins)[0]
     vel_min = vel_bins.min()
     gpu_array = tf.constant(psd)
@@ -55,9 +52,8 @@ def _gpu_calc_spectral_width(psd, power, vel_bins, velocity, dV):
 
     velocity_array = tf.transpose(np.tile(velocity, (shp[2], 1, 1)), [1, 2, 0])
     vel_bins_tiled = np.tile(vel_bins, (times, shp[1], 1))
-    gpu_array = tf.math.sqrt(1 / power_array *
-                        tf.math.reduce_sum(
-                            (vel_bins_tiled - velocity_array)**2 * gpu_array, axis=2))
+    gpu_array = tf.math.sqrt(1 / power_array * tf.math.reduce_sum(
+                             (vel_bins_tiled - velocity_array)**2 * gpu_array, axis=2))
     specwidth = gpu_array.numpy()
     return specwidth
 
@@ -73,16 +69,14 @@ def _gpu_calc_skewness(psd, power, vel_bins, velocity, spec_width, dV):
     velocity_array = tf.transpose(np.tile(velocity, (shp[2], 1, 1)), [1, 2, 0])
     vel_bins_tiled = np.tile(vel_bins, (times, shp[1], 1))
     gpu_array = 1 / power_array * tf.math.reduce_sum(
-            (vel_bins_tiled - velocity_array)**3 * gpu_array, axis=2)
+        (vel_bins_tiled - velocity_array)**3 * gpu_array, axis=2)
     skewness = gpu_array.numpy()
     return skewness
 
 
 def _gpu_calc_kurtosis(psd, power, vel_bins, velocity, spec_width, dV):
     shp = psd.shape
-    times = shp[0]
     kurtosis = np.zeros((shp[0], shp[1]))
-
     gpu_array = tf.constant(psd.values, dtype=tf.float32)
     power_array = tf.constant(power, dtype=tf.float32)
     spec_width_array = tf.constant(spec_width, dtype=tf.float32)
@@ -121,8 +115,7 @@ def get_lidar_moments(spectra, snr_thresh=0, block_size_ratio=1.0, which_moments
         The database with the Doppler lidar moments.
     """
     if 'power_spectral_density' not in spectra.variables.keys():
-        raise ValueError("You must calculate the power spectra " + 
-                         "before calculating moments!")
+        raise ValueError("You must calculate the power spectra before calculating moments!")
     if which_moments is None:
         which_moments = ['snr', 'doppler_velocity', 'spectral_width',
                          'skewness', 'kurtosis']
@@ -137,13 +130,8 @@ def get_lidar_moments(spectra, snr_thresh=0, block_size_ratio=1.0, which_moments
     linear_psd_0filled = linear_psd.fillna(0).values
     power = _gpu_calc_power(
         linear_psd_0filled, dV)
-    nyquist_range = spectra['vel_bins'].values.max() - spectra['vel_bins'].values.min()
-    if ('doppler_velocity' in which_moments or 'spectral_width' in which_moments or
-        'skewness' in which_moments or 'kurtosis' in which_moments):
-        velocity = _gpu_calc_velocity(
-            linear_psd_0filled, power,
-            spectra['vel_bins'].values, dV)
-
+    velocity = _gpu_calc_velocity(linear_psd_0filled, power, spectra['vel_bins'].values, dV)
+
     if 'snr' in which_moments:
         power_with_noise = dV * spectra['power_spectral_density'].sum(dim='vel_bins')
         spectra['snr'] = power_with_noise / (dV * len(spectra['vel_bins']))

highiq/calc/spectra.py

@@ -9,13 +9,12 @@
 
 def _fast_expand(complex_array, factor, num_per_block=200):
     shp = complex_array.shape
-    times = shp[0]
     expanded_array = np.zeros((shp[0], shp[1], shp[2] * factor))
     weights = np.tile(np.arange(0, factor) / factor, (shp[0], shp[1], 1))
-    for l in range(shp[2]):
-        gpu_array = tf.constant(complex_array[:, :, l], dtype=tf.float32)
-        if l < shp[2] - 1:
-            gpu_array2 = tf.constant(complex_array[:, :, l + 1], dtype=tf.float32)
+    for i in range(shp[2]):
+        gpu_array = tf.constant(complex_array[:, :, i], dtype=tf.float32)
+        if i < shp[2] - 1:
+            gpu_array2 = tf.constant(complex_array[:, :, i + 1], dtype=tf.float32)
             diff_array = gpu_array2 - gpu_array
         else:
             diff_array = tf.zeros((shp[0], shp[1]))
@@ -25,9 +24,10 @@ def _fast_expand(complex_array, factor, num_per_block=200):
         diff_array = tf.transpose(
             np.tile(diff_array, (factor, 1, 1)), [1, 2, 0])
         temp_array = diff_array * weights + rep_array
-        expanded_array[:, :, factor * l:factor * (l + 1)] = temp_array.numpy()
+        expanded_array[:, :, factor * i:factor * (i + 1)] = temp_array.numpy()
     return expanded_array
 
+
 def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, time_window=None,
             acf_name='acf', acf_bkg_name='acf_bkg', block_size_ratio=1.0):
     """
@@ -68,45 +68,37 @@ def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, t
 
     Q_ = UnitRegistry().Quantity
     if fs is None:
-        if not "sample_rate" in spectra.attrs:
-            raise KeyError("If sample frequency is not specified, then ACT Dataset must contain a sample_rate " +
-                           "attribute!")
+        if "sample_rate" not in spectra.attrs:
+            raise KeyError("If sample frequency is not specified, then ACT Dataset must contain a sample_rate attribute!")
 
         fs_pint = Q_(spectra.attrs["sample_rate"])
         fs_pint = fs_pint.to("Hz")
         fs = fs_pint.magnitude
 
-
     if wavelength is None:
-        if not "wavelength" in spectra.attrs:
+        if "wavelength" not in spectra.attrs:
             raise KeyError("If wavelength is not specified, then the dataset must contain a wavelength attribute!")
         fs_pint = Q_(spectra.attrs["wavelength"])
         fs_pint = fs_pint.to("m")
         wavelength = fs_pint.magnitude
-    
+
     if time_window is not None:
         spectra = spectra.resample(time='%ds' % int(time_window)).mean()
     else:
         spectra[acf_bkg_name] = xr.DataArray(np.ones(spectra[acf_name].shape),
-                                            dims=spectra[acf_name].dims) * spectra[acf_bkg_name]
+                                             dims=spectra[acf_name].dims) * spectra[acf_bkg_name]
 
     num_gates = int(gate_resolution / 3)
-    complex_coeff_in = spectra[acf_name].isel(complex=0).values +\
-                       spectra[acf_name].isel(complex=1).values * 1j
-    
+    complex_coeff_in = spectra[acf_name].isel(complex=0).values + \
+        spectra[acf_name].isel(complex=1).values * 1j
+
     complex_coeff = np.zeros(
-            (complex_coeff_in.shape[0], int(complex_coeff_in.shape[1] / num_gates),
-                complex_coeff_in.shape[2]), dtype=np.complex128)
+        (complex_coeff_in.shape[0], int(complex_coeff_in.shape[1] / num_gates),
+            complex_coeff_in.shape[2]), dtype=np.complex128)
     for i in range(complex_coeff.shape[1]):
-        complex_coeff[:, i, :] = np.mean(complex_coeff_in[
-                :, (num_gates * i):(num_gates * i+1), :], axis=1)
-
-
+        complex_coeff[:, i, :] = np.mean(complex_coeff_in[:, (num_gates * i):(num_gates * i + 1), :], axis=1)
     del complex_coeff_in
-
-    ntimes = complex_coeff.shape[0]
     freq = np.fft.fftfreq(nfft) * fs
-
     spectra.attrs['nyquist_velocity'] = "%f m s-1" % (wavelength / (4 * 1 / fs))
     spectra['freq_bins'] = xr.DataArray(freq, dims=['freq'])
     spectra['freq_bins'].attrs["long_name"] = "Doppler spectra bins in frequency units"
@@ -118,23 +110,21 @@ def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, t
     spectra['vel_bins'] = xr.DataArray(vel_bins, dims=('vel_bins'), attrs=attrs_dict)
     spectra['freq_bins'] = spectra['freq_bins'][inds_sorted]
 
-    complex_coeff_bkg_in = (spectra[acf_bkg_name].isel(complex=0).values +
-                        spectra[acf_bkg_name].isel(complex=1).values * 1j)
+    complex_coeff_bkg_in = (spectra[acf_bkg_name].isel(complex=0).values + spectra[acf_bkg_name].isel(complex=1).values * 1j)
 
     complex_coeff_bkg = np.zeros(
-            (complex_coeff_bkg_in.shape[0], int(complex_coeff_bkg_in.shape[1] / num_gates),
-                complex_coeff_bkg_in.shape[2]), dtype=np.complex128)
+        (complex_coeff_bkg_in.shape[0], int(complex_coeff_bkg_in.shape[1] / num_gates),
+            complex_coeff_bkg_in.shape[2]), dtype=np.complex128)
     for i in range(complex_coeff.shape[1]):
-        complex_coeff_bkg[:, i, :] = np.mean(complex_coeff_bkg_in[:,
-                (num_gates * i):(num_gates * i+1), :], axis=1)
+        complex_coeff_bkg[:, i, :] = np.mean(complex_coeff_bkg_in[:, (num_gates * i):(num_gates * i + 1), :], axis=1)
     num_lags = complex_coeff_bkg_in.shape[2]
     if nfft < num_lags:
         raise RuntimeError("Number of points in FFT < number of lags in sample!")
     pad_after = int((nfft - num_lags))
     pad_before = 0
     pad_lengths = tf.constant([[0, 0], [0, 0], [pad_before, pad_after]])
     frames = tf.pad(complex_coeff, pad_lengths, mode='CONSTANT', constant_values=0)
-   
+
     power = np.abs(tf.signal.fft(frames).numpy())
     attrs_dict = {'long_name': 'Range', 'units': 'm'}
     spectra['range'] = xr.DataArray(
@@ -143,18 +133,15 @@ def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, t
     spectra['power'] = xr.DataArray(
         power[:, :, inds_sorted], dims=(('time', 'range', 'vel_bins')))
     frames = tf.pad(complex_coeff_bkg, pad_lengths, mode='CONSTANT', constant_values=0)
-
     power = np.abs(tf.signal.fft(frames).numpy())
-   
+
     spectra['power_bkg'] = xr.DataArray(
         power[:, :, inds_sorted], dims=('time', 'range', 'vel_bins'))
 
     # Subtract background noise
     spectra['power_spectral_density'] = spectra['power'] / spectra['power_bkg']
     spectra['power_spectral_density'].attrs["long_name"] = "Power spectral density"
     spectra['power_spectral_density'].attrs["units"] = ''
-
-
     spectra['range'].attrs['long_name'] = "Range"
     spectra['range'].attrs['units'] = 'm'
     spectra['vel_bins'].attrs['long_name'] = "Doppler velocity"
@@ -186,40 +173,37 @@ def calc_num_peaks(my_spectra, **kwargs):
     shp = my_array.shape
     num_peaks = np.zeros((shp[0], shp[1]))
     peak_loc = np.nan * np.ones((shp[0], shp[1], 5))
-    
+
     vel_bins = my_spectra['vel_bins'].values
 
-    if not 'prominence' in kwargs.keys():
+    if 'prominence' not in kwargs.keys():
         prominence = 0.01
     else:
         prominence = kwargs.pop('prominence')
 
-    if not 'width' in kwargs.keys():
+    if 'width' not in kwargs.keys():
         width = 8
     else:
         width = kwargs.pop('width')
 
-    if not 'height' in kwargs.keys():
+    if 'height' not in kwargs.keys():
         height = 1.5
     else:
         height = kwargs.pop('height')
-        
+
     for i in range(shp[0]):
         for j in range(shp[1]):
-            max_spec = spectra[i, j].max()
-            peaks = find_peaks(my_array[i, j], prominence=prominence, width=width, **kwargs)[0]
+            peaks = find_peaks(my_array[i, j], prominence=prominence, width=width, height=height, **kwargs)[0]
             num_peaks[i, j] = len(peaks)
             for k in range(len(peaks)):
                 if k > 4:
                     continue
                 peak_loc[i, j, k] = vel_bins[peaks[k]]
 
-
     my_spectra['npeaks'] = xr.DataArray(num_peaks, dims=('time', 'range'))
     my_spectra['npeaks'].attrs['long_name'] = "Number of peaks in Doppler spectra"
     my_spectra['npeaks'].attrs['units'] = "1"
-
-    my_spectra['peak_velocities'] =  xr.DataArray(peak_loc, dims=('time', 'range', 'peak_no'))
+    my_spectra['peak_velocities'] = xr.DataArray(peak_loc, dims=('time', 'range', 'peak_no'))
     my_spectra['peak_velocities'].attrs['long_name'] = "Dopper velocity peaks"
     my_spectra['peak_velocities'].attrs['units'] = 'm s-1'
     return my_spectra

highiq/io/__init__.py

@@ -13,4 +13,4 @@
     load_arm_netcdf
 
 """
-from .arm_data import load_arm_netcdf, read_00_data
+from .arm_data import load_arm_netcdf, read_00_data