From e4b57bc7b491223646ab1832206877f76c36106f Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Mon, 12 Jun 2023 11:53:30 -0500 Subject: [PATCH] STY: Pep8 fixes --- doc/source/conf.py | 8 ++-- examples/plot_moments.py | 2 +- examples/plot_power_spectra.py | 6 +-- examples/read_00_file.py | 7 +++- highiq/calc/moments.py | 24 +++-------- highiq/calc/spectra.py | 76 ++++++++++++++-------------------- highiq/io/__init__.py | 2 +- highiq/io/arm_data.py | 27 ++++-------- highiq/testing/__init__.py | 2 +- tests/test_highiq.py | 2 +- 10 files changed, 60 insertions(+), 96 deletions(-) diff --git a/doc/source/conf.py b/doc/source/conf.py index 653b7d3..d06f076 100644 --- a/doc/source/conf.py +++ b/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 \ No newline at end of file +autosummary_generate = True diff --git a/examples/plot_moments.py b/examples/plot_moments.py index 110662f..c8476a7 100644 --- a/examples/plot_moments.py +++ b/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() \ No newline at end of file +test_file.close() diff --git a/examples/plot_power_spectra.py b/examples/plot_power_spectra.py index 3ede4cd..7ffba71 100644 --- a/examples/plot_power_spectra.py +++ b/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() \ No newline at end of file +test_file.close() diff --git a/examples/read_00_file.py b/examples/read_00_file.py index 6ab95e8..7825ac3 100644 --- a/examples/read_00_file.py +++ b/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) \ No newline at end of file +print(ds) diff --git a/highiq/calc/moments.py b/highiq/calc/moments.py index b46a24a..fc63c2a 100644 --- a/highiq/calc/moments.py +++ b/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'])) diff --git a/highiq/calc/spectra.py b/highiq/calc/spectra.py index 89d24db..32d488a 100644 --- a/highiq/calc/spectra.py +++ b/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,15 +110,13 @@ 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!") @@ -134,7 +124,7 @@ def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, t 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,9 +133,8 @@ 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')) @@ -153,8 +142,6 @@ def get_psd(spectra, gate_resolution=60., wavelength=None, fs=None, nfft=1024, t 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 diff --git a/highiq/io/__init__.py b/highiq/io/__init__.py index 3b19376..586f6b8 100644 --- a/highiq/io/__init__.py +++ b/highiq/io/__init__.py @@ -13,4 +13,4 @@ load_arm_netcdf """ -from .arm_data import load_arm_netcdf, read_00_data \ No newline at end of file +from .arm_data import load_arm_netcdf, read_00_data diff --git a/highiq/io/arm_data.py b/highiq/io/arm_data.py index da71b7d..587d41c 100644 --- a/highiq/io/arm_data.py +++ b/highiq/io/arm_data.py @@ -9,15 +9,10 @@ import numpy as np global_attrs = {'number_of_lags': { - 'cor.M1': 7, - 'hou.M1': 7, - 'oli.M1': 7, - 'sgp.C1': 20}, - 'number_of_samples_per_lag': { - 'cor.M1': 3200, - 'hou.M1': 3200, - 'oli.M1': 1920, - 'sgp.C1': 4000}} + 'cor.M1': 7, 'hou.M1': 7, 'oli.M1': 7, 'sgp.C1': 20}, + 'number_of_samples_per_lag': + {'cor.M1': 3200, 'hou.M1': 3200, 'oli.M1': 1920, 'sgp.C1': 4000}} + def load_arm_netcdf(arm_file, **kwargs): """ @@ -41,7 +36,7 @@ def load_arm_netcdf(arm_file, **kwargs): if 'time' not in ds['acf'].dims: ds['acf'] = ds['acf'].expand_dims('time') if 'time' not in ds['acf_bkg'].dims: - ds['acf_bkg'] = ds['acf_bkg'].expand_dims('time') + ds['acf_bkg'] = ds['acf_bkg'].expand_dims(dim={'time': len(ds.time)}, axis=0) return ds @@ -77,8 +72,7 @@ def read_00_data(file_name, home_point, site='sgp.C1', **kwargs): num_rows_to_skip += 1 hfile = pd.read_csv(home_point, skiprows=num_rows_to_skip, names=home_point_columns) - background_vals = global_attrs['number_of_lags'][site]\ - * global_attrs['number_of_samples_per_lag'][site] * 2 + background_vals = global_attrs['number_of_lags'][site] * global_attrs['number_of_samples_per_lag'][site] * 2 nlags = global_attrs['number_of_lags'][site] num_samples = global_attrs['number_of_samples_per_lag'][site] background_bytes = 8 * background_vals @@ -112,7 +106,6 @@ def read_00_data(file_name, home_point, site='sgp.C1', **kwargs): target_altitude = float(hfile['target_altitude'][i]) lidar_latitude = float(hfile['lidar_latitude'][i]) lidar_longitude = float(hfile['lidar_longitude'][i]) - lidar_altitude = float(hfile['lidar_altitude'][i]) lidar_home_point = float(hfile['lidar_home_point'][i]) if site[:3] == "oli": @@ -183,8 +176,7 @@ def read_00_data(file_name, home_point, site='sgp.C1', **kwargs): acf = xr.DataArray(acf, dims=('time', 'nsamples', 'nlags', 'complex')) acf.attrs["long_name"] = "Autocorrelation function" lidar_ds = xr.Dataset({'acf': acf, 'acf_bkg': acf_bkg, - 'aziumth': azimuth, 'time': timestamp, - 'elevation': elevation}) + 'aziumth': azimuth, 'time': timestamp, 'elevation': elevation}) lidar_ds = lidar_ds.sortby('time') lidar_ds.attrs["dlon"] = lidar_longitude lidar_ds.attrs["dlat"] = lidar_latitude @@ -195,8 +187,3 @@ def read_00_data(file_name, home_point, site='sgp.C1', **kwargs): lidar_ds.attrs["sample_rate"] = "50 MHz" lidar_ds.attrs["wavelength"] = "1548 nm" return lidar_ds - - - - - diff --git a/highiq/testing/__init__.py b/highiq/testing/__init__.py index 10c8bff..927455d 100644 --- a/highiq/testing/__init__.py +++ b/highiq/testing/__init__.py @@ -9,4 +9,4 @@ import os DATA_PATH = os.path.join(os.path.dirname(__file__), 'data') -TEST_FILE = os.path.join(DATA_PATH, 'testsig.a1.20180801.000003.nc') \ No newline at end of file +TEST_FILE = os.path.join(DATA_PATH, 'testsig.a1.20180801.000003.nc') diff --git a/tests/test_highiq.py b/tests/test_highiq.py index b03cfcc..9be3395 100644 --- a/tests/test_highiq.py +++ b/tests/test_highiq.py @@ -17,7 +17,7 @@ def test_spectra(): psd = my_spectra['power_spectral_density'].sel(range=400, method='nearest') vel_bins = my_spectra['vel_bins'] dV = vel_bins[1] - vel_bins[0] - np.testing.assert_almost_equal(psd.values.sum()*dV.values, 50.25995985957032) + np.testing.assert_almost_equal(psd.values.sum() * dV.values, 50.25995985957032) my_ds.close() my_spectra.close()