From 200df500f754e13cb90d277b272a09f866a26522 Mon Sep 17 00:00:00 2001 From: Adam Theisen Date: Tue, 23 Aug 2022 15:50:06 -0500 Subject: [PATCH] Adding Reader for the NOAA FMCW data (#518) * ENH: Updating example to try and get it to run * ADD: Adding a reader for the NOAA FMCW radar --- act/io/__init__.py | 7 +- act/io/noaapsl.py | 164 ++++++++++++++++++++++++++++++ act/tests/test_io.py | 11 ++ examples/plot_noaa_fmcw_moment.py | 29 ++++++ 4 files changed, 210 insertions(+), 1 deletion(-) create mode 100644 examples/plot_noaa_fmcw_moment.py diff --git a/act/io/__init__.py b/act/io/__init__.py index 7cc2eec1f2..87c44cab1d 100644 --- a/act/io/__init__.py +++ b/act/io/__init__.py @@ -14,5 +14,10 @@ read_gml_ozone, read_gml_radiation, ) -from .noaapsl import read_psl_wind_profiler, read_psl_wind_profiler_temperature, read_psl_parsivel +from .noaapsl import ( + read_psl_wind_profiler, + read_psl_wind_profiler_temperature, + read_psl_parsivel, + read_psl_radar_fmcw_moment +) from .pysp2 import read_hk_file, read_sp2, read_sp2_dat diff --git a/act/io/noaapsl.py b/act/io/noaapsl.py index 8868ff3226..4f03abbb73 100644 --- a/act/io/noaapsl.py +++ b/act/io/noaapsl.py @@ -9,6 +9,7 @@ import numpy as np import pandas as pd import xarray as xr +import datetime as dt def read_psl_wind_profiler(filename, transpose=True): @@ -459,3 +460,166 @@ def read_psl_parsivel(files): obj[v].attrs = attrs[v] return obj + + +def read_psl_radar_fmcw_moment(files): + """ + Returns `xarray.Dataset` with stored data and metadata from + NOAA PSL FMCW Radar files. See References section for details. + + Parameters + ---------- + files : str or list + Name of file(s) to read. Currently does not support reading URLs but files can + be downloaded easily using the act.discovery.download_noaa_psl_data function. + + Return + ------ + obj : Xarray.dataset + Standard Xarray dataset with the data for the parsivel + + References + ---------- + Johnston, Paul E., James R. Jordan, Allen B. White, David A. Carter, David M. Costa, and Thomas E. Ayers. + "The NOAA FM-CW snow-level radar." Journal of Atmospheric and Oceanic Technology 34, no. 2 (2017): 249-267. + + """ + + # Set the initial dictionary to convert to xarray dataset + data = { + 'site': {'dims': ['file'], 'data': [], 'attrs': {'long_name': 'NOAA site code'}}, + 'lat': {'dims': ['file'], 'data': [], 'attrs': {'long_name': 'North Latitude', 'units': 'degree_N'}}, + 'lon': {'dims': ['file'], 'data': [], 'attrs': {'long_name': 'East Longitude', 'units': 'degree_E'}}, + 'alt': {'dims': ['file'], 'data': [], 'attrs': {'long_name': 'Altitude above mean sea level', 'units': 'm'}}, + 'freq': {'dims': ['file'], 'data': [], 'attrs': {'long_name': 'Operating Frequency; Ignore for FMCW', 'units': 'Hz'}}, + 'azimuth': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Azimuth angle', 'units': 'deg'}}, + 'elevation': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Elevation angle', 'units': 'deg'}}, + 'beam_direction_code': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Beam direction code', 'units': ''}}, + 'year': {'dims': ['time'], 'data': [], 'attrs': {'long_name': '2-digit year', 'units': ''}}, + 'day_of_year': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Day of the year', 'units': ''}}, + 'hour': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Hour of the day', 'units': ''}}, + 'minute': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Minutes', 'units': ''}}, + 'second': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Seconds', 'units': ''}}, + 'interpulse_period': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Interpulse Period', 'units': 'ms'}}, + 'pulse_width': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Pulse width', 'units': 'ns'}}, + 'first_range_gate': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Range to first range gate', 'units': 'm'}}, + 'range_between_gates': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Distance between range gates', 'units': 'm'}}, + 'n_gates': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Number of range gates', 'units': 'count'}}, + 'n_coherent_integration': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Number of cohrent integration', 'units': 'count'}}, + 'n_averaged_spectra': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Number of average spectra', 'units': 'count'}}, + 'n_points_spectrum': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Number of points in spectra', 'units': 'count'}}, + 'n_code_bits': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Number of code bits', 'units': 'count'}}, + 'radial_velocity': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'Radial velocity', 'units': 'm/s'}}, + 'snr': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'Signal-to-noise ratio - not range corrected', 'units': 'dB'}}, + 'signal_power': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'Signal Power - not range corrected', 'units': 'dB'}}, + 'spectral_width': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'Spectral width', 'units': 'm/s'}}, + 'noise_amplitude': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'noise_amplitude', 'units': 'dB'}}, + 'qc_variable': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'QC Value - not used', 'units': ''}}, + 'time': {'dims': ['time'], 'data': [], 'attrs': {'long_name': 'Datetime', 'units': ''}}, + 'range': {'dims': ['range'], 'data': [], 'attrs': {'long_name': 'Range', 'units': 'm'}}, + 'reflectivity_uncalibrated': {'dims': ['time', 'range'], 'data': [], 'attrs': {'long_name': 'Range', 'units': 'dB'}}, + } + + # Separate out the names as they will be accessed in different parts of the code + h1_names = ['site', 'lat', 'lon', 'alt', 'freq'] + h2_names = ['azimuth', 'elevation', 'beam_direction_code', 'year', 'day_of_year', 'hour', 'minute', + 'second'] + h3_names = ['interpulse_period', 'pulse_width', 'first_range_gate', 'range_between_gates', + 'n_gates', 'n_coherent_integration', 'n_averaged_spectra', 'n_points_spectrum', + 'n_code_bits'] + names = {'radial_velocity': 0, 'snr': 1, 'signal_power': 2, 'spectral_width': 3, + 'noise_amplitude': 4, 'qc_variable': 5} + + # Run through each file and read the data in + for f in files: + # Read in the first line of the file which has site, lat, lon, etc... + df = str(pd.read_table(f, nrows=0).columns[0]).split(' ') + ctr = 0 + for d in df: + if len(d) > 0: + if d == 'lat' or d == 'lon': + data[h1_names[ctr]]['data'].append(float(d) / 100.) + else: + data[h1_names[ctr]]['data'].append(d) + ctr += 1 + + # Set counts and errors + error = False + error_ct = 0 + ct = 0 + # Loop through while there's no errors i.e. eof + while error is False: + try: + # Read in the initial headers to get information used to parse data + if ct == 0: + df = str(pd.read_table(f, nrows=0, skiprows=[0]).columns[0]).split(' ') + ctr = 0 + for d in df: + if len(d) > 0: + data[h2_names[ctr]]['data'].append(d) + ctr += 1 + # Read in third row of header information + df = str(pd.read_table(f, nrows=0, skiprows=[0, 1]).columns[0]).split(' ') + ctr = 0 + for d in df: + if len(d) > 0: + data[h3_names[ctr]]['data'].append(d) + ctr += 1 + # Read in the data based on number of gates + df = pd.read_csv(f, skiprows=[0, 1, 2], nrows=int(data['n_gates']['data'][-1]) - 1, delim_whitespace=True, + names=list(names.keys())) + index2 = 0 + else: + # Set indices for parsing data, reading 2 headers and then the columns of data + index1 = ct * int(data['n_gates']['data'][-1]) + index2 = index1 + int(data['n_gates']['data'][-1]) + 2 * ct + 4 + df = str(pd.read_table(f, nrows=0, skiprows=list(range(index2 - 1))).columns[0]).split(' ') + ctr = 0 + for d in df: + if len(d) > 0: + data[h2_names[ctr]]['data'].append(d) + ctr += 1 + df = str(pd.read_table(f, nrows=0, skiprows=list(range(index2))).columns[0]).split(' ') + ctr = 0 + for d in df: + if len(d) > 0: + data[h3_names[ctr]]['data'].append(d) + ctr += 1 + df = pd.read_csv(f, skiprows=list(range(index2 + 1)), nrows=int(data['n_gates']['data'][-1]) - 1, delim_whitespace=True, + names=list(names.keys())) + + # Add data from the columns to the dictionary + for n in names: + data[n]['data'].append(df[n].to_list()) + + # Calculate the range based on number of gates, range to first gate and range between gates + if len(data['range']['data']) == 0: + ranges = float(data['first_range_gate']['data'][-1]) + np.array(range(int(data['n_gates']['data'][-1]) - 1)) * \ + float(data['range_between_gates']['data'][-1]) + data['range']['data'] = ranges + + # Calculate a time + time = dt.datetime( + int('20' + data['year']['data'][-1]), 1, 1, int(data['hour']['data'][-1]), + int(data['minute']['data'][-1]), int(data['second']['data'][-1])) + \ + dt.timedelta(days=int(data['day_of_year']['data'][-1])) + data['time']['data'].append(time) + + # Range correct the snr which converts it essentially to an uncalibrated reflectivity + snr_rc = data['snr']['data'][-1] - 20. * np.log10(1. / (ranges / 1000.) ** 2) + data['reflectivity_uncalibrated']['data'].append(snr_rc) + except Exception as e: + # Handle errors, if end of file then continue on, if something else + # try the next block of data but if it errors another time in this file move on + if isinstance(e, pd.errors.EmptyDataError) or error_ct > 1: + error = True + else: + print(e) + pass + error_ct += 1 + ct += 1 + + # Convert dictionary to Dataset + obj = xr.Dataset().from_dict(data) + + return obj diff --git a/act/tests/test_io.py b/act/tests/test_io.py index 67fe6234af..560742cd1e 100644 --- a/act/tests/test_io.py +++ b/act/tests/test_io.py @@ -526,3 +526,14 @@ def test_read_psl_parsivel(): obj = act.io.noaapsl.read_psl_parsivel('https://downloads.psl.noaa.gov/psd2/data/realtime/DisdrometerParsivel/Stats/ctd/2022/002/ctd2200201_stats.txt') assert 'number_density_drops' in obj + + +def test_read_psl_fmcw_moment(): + result = act.discovery.download_noaa_psl_data( + site='kps', instrument='Radar FMCW Moment', startdate='20220815', hour='06' + ) + obj = act.io.noaapsl.read_psl_radar_fmcw_moment([result[-1]]) + assert 'range' in obj + np.testing.assert_almost_equal(obj['reflectivity_uncalibrated'].mean(), 2.37, decimal=2) + assert obj['range'].max() == 10040. + assert len(obj['time'].values) == 115 diff --git a/examples/plot_noaa_fmcw_moment.py b/examples/plot_noaa_fmcw_moment.py new file mode 100644 index 0000000000..26c897bb60 --- /dev/null +++ b/examples/plot_noaa_fmcw_moment.py @@ -0,0 +1,29 @@ +""" +NOAA FMCW Moment Data +--------------------- + +ARM and NOAA have campaigns going on in the Crested Butte, CO region +and as part of that campaign NOAA has FMCW radars deployed that could +benefit the broader ARM and NOAA communities. This example shows how +easy it is to download and read in NOAA PSL data. + +""" + +import act +import os +import matplotlib.pyplot as plt + +# Use the ACT downloader to download a file from the +# Kettle Ponds site on 8/15/2022 at 2300 UTC +result = act.discovery.download_noaa_psl_data( + site='kps', instrument='Radar FMCW Moment', startdate='20220815', hour='23' +) + +# Read in the .raw file. Spectra data are also downloaded +obj = act.io.noaapsl.read_psl_radar_fmcw_moment([result[-1]]) + +# As Part of the reading in, the script calculates ranges and +# corrects the SNR for range +display = act.plotting.TimeSeriesDisplay(obj) +display.plot('reflectivity_uncalibrated', cmap='jet', vmin=-20, vmax=40) +plt.show()