From c9d88ec7ff03adff0f77fa3c1beca5a6ceecde17 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 5 Jan 2024 15:08:57 -0600 Subject: [PATCH 01/48] ENH: Add first cut at nexrad reader --- xradar/io/backends/__init__.py | 2 + xradar/io/backends/common.py | 211 +++++ xradar/io/backends/nexrad_archive.py | 287 +++++++ xradar/io/backends/nexrad_common.py | 247 ++++++ xradar/io/backends/nexrad_interpolate.pyx | 110 +++ xradar/io/backends/nexrad_level2.py | 960 ++++++++++++++++++++++ 6 files changed, 1817 insertions(+) create mode 100644 xradar/io/backends/nexrad_archive.py create mode 100644 xradar/io/backends/nexrad_common.py create mode 100644 xradar/io/backends/nexrad_interpolate.pyx create mode 100644 xradar/io/backends/nexrad_level2.py diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 4edf3d6..509b11f 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -24,5 +24,7 @@ from .iris import * # noqa from .odim import * # noqa from .rainbow import * # noqa +from .nexrad_level2 import * # noqa +from .nexrad_archive import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/common.py b/xradar/io/backends/common.py index d954bce..b338498 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -12,11 +12,17 @@ """ +import bz2 +import gzip import io +import itertools import struct from collections import OrderedDict +from collections.abc import MutableMapping +import fsspec import h5netcdf +import netCDF4 import numpy as np import xarray as xr from datatree import DataTree @@ -229,3 +235,208 @@ def _unpack_dictionary(buffer, dictionary, rawdata=False): UINT1 = {"fmt": "B", "dtype": "unit8"} UINT2 = {"fmt": "H", "dtype": "uint16"} UINT4 = {"fmt": "I", "dtype": "unint32"} + + +def prepare_for_read(filename, storage_options={"anon": True}): + """ + Return a file like object read for reading. + + Open a file for reading in binary mode with transparent decompression of + Gzip and BZip2 files. The resulting file-like object should be closed. + + Parameters + ---------- + filename : str or file-like object + Filename or file-like object which will be opened. File-like objects + will not be examined for compressed data. + + storage_options : dict, optional + Parameters passed to the backend file-system such as Google Cloud Storage, + Amazon Web Service S3. + + Returns + ------- + file_like : file-like object + File like object from which data can be read. + + """ + # if a file-like object was provided, return + if hasattr(filename, "read"): # file-like object + return filename + + # look for compressed data by examining the first few bytes + fh = fsspec.open(filename, mode="rb", compression="infer", **storage_options).open() + magic = fh.read(3) + fh.close() + + # If the data is still compressed, use gunzip/bz2 to uncompress the data + if magic.startswith(b"\x1f\x8b"): + return gzip.GzipFile(filename, "rb") + + if magic.startswith(b"BZh"): + return bz2.BZ2File(filename, "rb") + + return fsspec.open( + filename, mode="rb", compression="infer", **storage_options + ).open() + + +def stringarray_to_chararray(arr, numchars=None): + """ + Convert an string array to a character array with one extra dimension. + + Parameters + ---------- + arr : array + Array with numpy dtype 'SN', where N is the number of characters + in the string. + + numchars : int + Number of characters used to represent the string. If numchar > N + the results will be padded on the right with blanks. The default, + None will use N. + + Returns + ------- + chararr : array + Array with dtype 'S1' and shape = arr.shape + (numchars, ). + + """ + carr = netCDF4.stringtochar(arr) + if numchars is None: + return carr + + arr_numchars = carr.shape[-1] + if numchars <= arr_numchars: + raise ValueError("numchars must be >= %i" % (arr_numchars)) + chararr = np.zeros(arr.shape + (numchars,), dtype="S1") + chararr[..., :arr_numchars] = carr[:] + return chararr + + +def _test_arguments(dic): + """Issue a warning if receive non-empty argument dict.""" + if dic: + import warnings + + warnings.warn("Unexpected arguments: %s" % dic.keys()) + + +def make_time_unit_str(dtobj): + """Return a time unit string from a datetime object.""" + return "seconds since " + dtobj.strftime("%Y-%m-%dT%H:%M:%SZ") + + +class LazyLoadDict(MutableMapping): + """ + A dictionary-like class supporting lazy loading of specified keys. + + Keys which are lazy loaded are specified using the set_lazy method. + The callable object which produces the specified key is provided as the + second argument to this method. This object gets called when the value + of the key is loaded. After this initial call the results is cached + in the traditional dictionary which is used for supplemental access to + this key. + + Testing for keys in this dictionary using the "key in d" syntax will + result in the loading of a lazy key, use "key in d.keys()" to prevent + this evaluation. + + The comparison methods, __cmp__, __ge__, __gt__, __le__, __lt__, __ne__, + nor the view methods, viewitems, viewkeys, viewvalues, are implemented. + Neither is the the fromkeys method. + + This is from Py-ART. + + Parameters + ---------- + dic : dict + Dictionary containing key, value pairs which will be stored and + evaluated traditionally. This dictionary referenced not copied into + the LazyLoadDictionary and hence changed to this dictionary may change + the original. If this behavior is not desired copy dic in the + initalization. + + Examples + -------- + >>> d = LazyLoadDict({'key1': 'value1', 'key2': 'value2'}) + >>> d.keys() + ['key2', 'key1'] + >>> lazy_func = lambda : 999 + >>> d.set_lazy('lazykey1', lazy_func) + >>> d.keys() + ['key2', 'key1', 'lazykey1'] + >>> d['lazykey1'] + 999 + + """ + + def __init__(self, dic): + """initalize.""" + self._dic = dic + self._lazyload = {} + + # abstract methods + def __setitem__(self, key, value): + """Set a key which will not be stored and evaluated traditionally.""" + self._dic[key] = value + if key in self._lazyload: + del self._lazyload[key] + + def __getitem__(self, key): + """Get the value of a key, evaluating a lazy key if needed.""" + if key in self._lazyload: + value = self._lazyload[key]() + self._dic[key] = value + del self._lazyload[key] + return self._dic[key] + + def __delitem__(self, key): + """Remove a lazy or traditional key from the dictionary.""" + if key in self._lazyload: + del self._lazyload[key] + else: + del self._dic[key] + + def __iter__(self): + """Iterate over all lazy and traditional keys.""" + return itertools.chain(self._dic.copy(), self._lazyload.copy()) + + def __len__(self): + """Return the number of traditional and lazy keys.""" + return len(self._dic) + len(self._lazyload) + + # additional class to mimic dict behavior + def __str__(self): + """Return a string representation of the object.""" + if len(self._dic) == 0 or len(self._lazyload) == 0: + seperator = "" + else: + seperator = ", " + lazy_reprs = [(repr(k), repr(v)) for k, v in self._lazyload.items()] + lazy_strs = ["{}: LazyLoad({})".format(*r) for r in lazy_reprs] + lazy_str = ", ".join(lazy_strs) + "}" + return str(self._dic)[:-1] + seperator + lazy_str + + def has_key(self, key): + """True if dictionary has key, else False.""" + return key in self + + def copy(self): + """ + Return a copy of the dictionary. + + Lazy keys are not evaluated in the original or copied dictionary. + """ + dic = self.__class__(self._dic.copy()) + # load all lazy keys into the copy + for key, value_callable in self._lazyload.items(): + dic.set_lazy(key, value_callable) + return dic + + # lazy dictionary specific methods + def set_lazy(self, key, value_callable): + """Set a lazy key to load from a callable object.""" + if key in self._dic: + del self._dic[key] + self._lazyload[key] = value_callable diff --git a/xradar/io/backends/nexrad_archive.py b/xradar/io/backends/nexrad_archive.py new file mode 100644 index 0000000..d65fc3d --- /dev/null +++ b/xradar/io/backends/nexrad_archive.py @@ -0,0 +1,287 @@ +import warnings + +import numpy as np + +from ...model import ( + get_altitude_attrs, + get_azimuth_attrs, + get_elevation_attrs, + get_latitude_attrs, + get_longitude_attrs, + get_moment_attrs, + get_range_attrs, + get_time_attrs, +) +from .common import LazyLoadDict, prepare_for_read +from .nexrad_common import get_nexrad_location +from .nexrad_interpolate import _fast_interpolate_scan_2, _fast_interpolate_scan_4 +from .nexrad_level2 import NEXRADLevel2File + +nexrad_mapping = { + "REF": "DBZH", + "VEL": "VRADH", + "SW": "WRADH", + "ZDR": "ZDR", + "PHI": "PHIDP", + "RHO": "RHOHV", + "CFP": "CCORH", +} + + +class _NEXRADLevel2StagedField: + """ + A class to facilitate on demand loading of field data from a Level 2 file. + """ + + def __init__(self, nfile, moment, max_ngates, scans): + """initialize.""" + self.nfile = nfile + self.moment = moment + self.max_ngates = max_ngates + self.scans = scans + + def __call__(self): + """Return the array containing the field data.""" + return self.nfile.get_data(self.moment, self.max_ngates, scans=self.scans) + + +def _find_range_params(scan_info): + """Return range parameters, first_gate, gate_spacing, last_gate.""" + min_first_gate = 999999 + min_gate_spacing = 999999 + max_last_gate = 0 + for scan_params in scan_info: + ngates = scan_params["ngates"][0] + for i, moment in enumerate(scan_params["moments"]): + first_gate = scan_params["first_gate"][i] + gate_spacing = scan_params["gate_spacing"][i] + last_gate = first_gate + gate_spacing * (ngates - 0.5) + + min_first_gate = min(min_first_gate, first_gate) + min_gate_spacing = min(min_gate_spacing, gate_spacing) + max_last_gate = max(max_last_gate, last_gate) + return min_first_gate, min_gate_spacing, max_last_gate + + +def _find_scans_to_interp(scan_info, first_gate, gate_spacing): + """Return a dict indicating what moments/scans need interpolation.""" + moments = {m for scan in scan_info for m in scan["moments"]} + interpolate = {moment: [] for moment in moments} + for scan_num, scan in enumerate(scan_info): + for moment in moments: + if moment not in scan["moments"]: + continue + index = scan["moments"].index(moment) + first = scan["first_gate"][index] + spacing = scan["gate_spacing"][index] + if first != first_gate or spacing != gate_spacing: + interpolate[moment].append(scan_num) + # for proper interpolation the gate spacing of the scan to be + # interpolated should be 1/4th the spacing of the radar + if spacing == gate_spacing * 4: + interpolate["multiplier"] = "4" + elif spacing == gate_spacing * 2: + interpolate["multiplier"] = "2" + else: + raise ValueError("Gate spacing is neither 1/4 or 1/2") + # assert first_gate + 1.5 * gate_spacing == first + # remove moments with no scans needing interpolation + interpolate = {k: v for k, v in interpolate.items() if len(v) != 0} + return interpolate + + +def _interpolate_scan(mdata, start, end, moment_ngates, multiplier, linear_interp=True): + """Interpolate a single NEXRAD moment scan from 1000 m to 250 m.""" + fill_value = -9999 + data = mdata.filled(fill_value) + scratch_ray = np.empty((data.shape[1],), dtype=data.dtype) + if multiplier == "4": + _fast_interpolate_scan_4( + data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp + ) + else: + _fast_interpolate_scan_2( + data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp + ) + mdata[:] = np.ma.array(data, mask=(data == fill_value)) + + +def open_dataset( + filename, + field_names=None, + additional_metadata=None, + file_field_names=False, + exclude_fields=None, + include_fields=None, + delay_field_loading=False, + station=None, + scans=None, + linear_interp=True, + storage_options={"anon": True}, + **kwargs, +): + # Load the data file in using NEXRADLevel2File Class + nfile = NEXRADLevel2File(prepare_for_read(filename)) + + # Access the scan info and load in the available moments + scan_info = nfile.scan_info(scans) + available_moments = {m for scan in scan_info for m in scan["moments"]} + first_gate, gate_spacing, last_gate = _find_range_params(scan_info) + + # Interpolate to 360 degrees where neccessary + interpolate = _find_scans_to_interp(scan_info, first_gate, gate_spacing) + + # Deal with time + time_start, _time = nfile.get_times(scans) + time = get_time_attrs(time_start) + time["data"] = _time + + # range + _range = get_range_attrs() + first_gate, gate_spacing, last_gate = _find_range_params(scan_info) + _range["data"] = np.arange(first_gate, last_gate, gate_spacing, "float32") + _range["meters_to_center_of_first_gate"] = float(first_gate) + _range["meters_between_gates"] = float(gate_spacing) + + # metadata + metadata = { + "Conventions": "CF/Radial instrument_parameters", + "version": "1.3", + "title": "", + "institution": "", + "references": "", + "source": "", + "history": "", + "comment": "", + "intrument_name": "", + "original_container": "NEXRAD Level II", + } + + vcp_pattern = nfile.get_vcp_pattern() + if vcp_pattern is not None: + metadata["vcp_pattern"] = vcp_pattern + if "icao" in nfile.volume_header.keys(): + metadata["instrument_name"] = nfile.volume_header["icao"].decode() + + # scan_type + + # latitude, longitude, altitude + latitude = get_latitude_attrs() + longitude = get_longitude_attrs() + altitude = get_altitude_attrs() + + if nfile._msg_type == "1" and station is not None: + lat, lon, alt = get_nexrad_location(station) + elif ( + "icao" in nfile.volume_header.keys() + and nfile.volume_header["icao"].decode()[0] == "T" + ): + lat, lon, alt = get_nexrad_location(nfile.volume_header["icao"].decode()) + else: + lat, lon, alt = nfile.location() + latitude["data"] = np.array([lat], dtype="float64") + longitude["data"] = np.array([lon], dtype="float64") + altitude["data"] = np.array([alt], dtype="float64") + + # Sweep information + sweep_number = { + "units": "count", + "standard_name": "sweep_number", + "long_name": "Sweep number", + } + + sweep_mode = { + "units": "unitless", + "standard_name": "sweep_mode", + "long_name": "Sweep mode", + "comment": 'Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi"', + } + sweep_start_ray_index = { + "long_name": "Index of first ray in sweep, 0-based", + "units": "count", + } + sweep_end_ray_index = { + "long_name": "Index of last ray in sweep, 0-based", + "units": "count", + } + + if scans is None: + nsweeps = int(nfile.nscans) + else: + nsweeps = len(scans) + sweep_number["data"] = np.arange(nsweeps, dtype="int32") + sweep_mode["data"] = np.array(nsweeps * ["azimuth_surveillance"], dtype="S") + + rays_per_scan = [s["nrays"] for s in scan_info] + sweep_end_ray_index["data"] = np.cumsum(rays_per_scan, dtype="int32") - 1 + + rays_per_scan.insert(0, 0) + sweep_start_ray_index["data"] = np.cumsum(rays_per_scan[:-1], dtype="int32") + + # azimuth, elevation, fixed_angle + azimuth = get_azimuth_attrs() + elevation = get_elevation_attrs() + fixed_angle = { + "long_name": "Target angle for sweep", + "units": "degrees", + "standard_name": "target_fixed_angle", + } + + azimuth["data"] = nfile.get_azimuth_angles(scans) + elevation["data"] = nfile.get_elevation_angles(scans).astype("float32") + fixed_agl = [] + for i in nfile.get_target_angles(scans): + if i > 180: + i = i - 360.0 + warnings.warn( + "Fixed_angle(s) greater than 180 degrees present." + " Assuming angle to be negative so subtrating 360", + UserWarning, + ) + else: + i = i + fixed_agl.append(i) + fixed_angles = np.array(fixed_agl, dtype="float32") + fixed_angle["data"] = fixed_angles + + # fields + max_ngates = len(_range["data"]) + available_moments = {m for scan in scan_info for m in scan["moments"]} + interpolate = _find_scans_to_interp( + scan_info, + first_gate, + gate_spacing, + ) + + fields = {} + for moment in available_moments: + dic = get_moment_attrs(nexrad_mapping[moment]) + dic["_FillValue"] = -9999 + if delay_field_loading and moment not in interpolate: + dic = LazyLoadDict(dic) + data_call = _NEXRADLevel2StagedField(nfile, moment, max_ngates, scans) + dic.set_lazy("data", data_call) + else: + mdata = nfile.get_data(moment, max_ngates, scans=scans) + if moment in interpolate: + interp_scans = interpolate[moment] + warnings.warn( + "Gate spacing is not constant, interpolating data in " + + f"scans {interp_scans} for moment {moment}.", + UserWarning, + ) + for scan in interp_scans: + idx = scan_info[scan]["moments"].index(moment) + moment_ngates = scan_info[scan]["ngates"][idx] + start = sweep_start_ray_index["data"][scan] + end = sweep_end_ray_index["data"][scan] + if interpolate["multiplier"] == "4": + multiplier = "4" + else: + multiplier = "2" + _interpolate_scan( + mdata, start, end, moment_ngates, multiplier, linear_interp + ) + dic["data"] = mdata + fields[nexrad_mapping[moment]] = dic + return fields diff --git a/xradar/io/backends/nexrad_common.py b/xradar/io/backends/nexrad_common.py new file mode 100644 index 0000000..e34c277 --- /dev/null +++ b/xradar/io/backends/nexrad_common.py @@ -0,0 +1,247 @@ +""" +Data and functions common to all types of NEXRAD files. + +""" +# The functions in this module are intended to be used in other +# nexrad related modules. The functions are not and should not be +# exported into the pyart.io namespace. + + +def get_nexrad_location(station): + """ + Return the latitude, longitude and altitude of a NEXRAD station. + + Parameters + ---------- + station : str + Four letter NEXRAD station ICAO name. + + Returns + ------- + lat, lon, alt : float + Latitude (in degrees), longitude (in degrees), and altitude + (in meters above mean sea level) of the NEXRAD station. + + """ + loc = NEXRAD_LOCATIONS[station.upper()] + + # Convert from feet to meters for elevation units + loc["elev"] = loc["elev"] * 0.3048 + + return loc["lat"], loc["lon"], loc["elev"] + + +# Locations of NEXRAD locations was retrieved from NOAA's +# Historical Observing Metadata Repository (HOMR) on +# 2014-Mar-27. http://www.ncdc.noaa.gov/homr/ +# The data below was extracted with: +# cut -c 10-14,107-115,117-127,128-133 + +NEXRAD_LOCATIONS = { + "KABR": {"lat": 45.45583, "lon": -98.41306, "elev": 1302}, + "KABX": {"lat": 35.14972, "lon": -106.82333, "elev": 5870}, + "KAKQ": {"lat": 36.98389, "lon": -77.0075, "elev": 112}, + "KAMA": {"lat": 35.23333, "lon": -101.70889, "elev": 3587}, + "KAMX": {"lat": 25.61056, "lon": -80.41306, "elev": 14}, + "KAPX": {"lat": 44.90722, "lon": -84.71972, "elev": 1464}, + "KARX": {"lat": 43.82278, "lon": -91.19111, "elev": 1276}, + "KATX": {"lat": 48.19472, "lon": -122.49444, "elev": 494}, + "KBBX": {"lat": 39.49611, "lon": -121.63167, "elev": 173}, + "KBGM": {"lat": 42.19972, "lon": -75.985, "elev": 1606}, + "KBHX": {"lat": 40.49833, "lon": -124.29194, "elev": 2402}, + "KBIS": {"lat": 46.77083, "lon": -100.76028, "elev": 1658}, + "KBLX": {"lat": 45.85389, "lon": -108.60611, "elev": 3598}, + "KBMX": {"lat": 33.17194, "lon": -86.76972, "elev": 645}, + "KBOX": {"lat": 41.95583, "lon": -71.1375, "elev": 118}, + "KBRO": {"lat": 25.91556, "lon": -97.41861, "elev": 23}, + "KBUF": {"lat": 42.94861, "lon": -78.73694, "elev": 693}, + "KBYX": {"lat": 24.59694, "lon": -81.70333, "elev": 8}, + "KCAE": {"lat": 33.94861, "lon": -81.11861, "elev": 231}, + "KCBW": {"lat": 46.03917, "lon": -67.80694, "elev": 746}, + "KCBX": {"lat": 43.49083, "lon": -116.23444, "elev": 3061}, + "KCCX": {"lat": 40.92306, "lon": -78.00389, "elev": 2405}, + "KCLE": {"lat": 41.41306, "lon": -81.86, "elev": 763}, + "KCLX": {"lat": 32.65556, "lon": -81.04222, "elev": 97}, + "KCRI": {"lat": 35.2383, "lon": -97.4602, "elev": 1201}, + "KCRP": {"lat": 27.78389, "lon": -97.51083, "elev": 45}, + "KCXX": {"lat": 44.51111, "lon": -73.16639, "elev": 317}, + "KCYS": {"lat": 41.15194, "lon": -104.80611, "elev": 6128}, + "KDAX": {"lat": 38.50111, "lon": -121.67667, "elev": 30}, + "KDDC": {"lat": 37.76083, "lon": -99.96833, "elev": 2590}, + "KDFX": {"lat": 29.2725, "lon": -100.28028, "elev": 1131}, + "KDGX": {"lat": 32.28, "lon": -89.98444, "elev": -99999}, + "KDIX": {"lat": 39.94694, "lon": -74.41111, "elev": 149}, + "KDLH": {"lat": 46.83694, "lon": -92.20972, "elev": 1428}, + "KDMX": {"lat": 41.73111, "lon": -93.72278, "elev": 981}, + "KDOX": {"lat": 38.82556, "lon": -75.44, "elev": 50}, + "KDTX": {"lat": 42.69972, "lon": -83.47167, "elev": 1072}, + "KDVN": {"lat": 41.61167, "lon": -90.58083, "elev": 754}, + "KDYX": {"lat": 32.53833, "lon": -99.25417, "elev": 1517}, + "KEAX": {"lat": 38.81028, "lon": -94.26417, "elev": 995}, + "KEMX": {"lat": 31.89361, "lon": -110.63028, "elev": 5202}, + "KENX": {"lat": 42.58639, "lon": -74.06444, "elev": 1826}, + "KEOX": {"lat": 31.46028, "lon": -85.45944, "elev": 434}, + "KEPZ": {"lat": 31.87306, "lon": -106.6975, "elev": 4104}, + "KESX": {"lat": 35.70111, "lon": -114.89139, "elev": 4867}, + "KEVX": {"lat": 30.56417, "lon": -85.92139, "elev": 140}, + "KEWX": {"lat": 29.70361, "lon": -98.02806, "elev": 633}, + "KEYX": {"lat": 35.09778, "lon": -117.56, "elev": 2757}, + "KFCX": {"lat": 37.02417, "lon": -80.27417, "elev": 2868}, + "KFDR": {"lat": 34.36222, "lon": -98.97611, "elev": 1267}, + "KFDX": {"lat": 34.63528, "lon": -103.62944, "elev": 4650}, + "KFFC": {"lat": 33.36333, "lon": -84.56583, "elev": 858}, + "KFSD": {"lat": 43.58778, "lon": -96.72889, "elev": 1430}, + "KFSX": {"lat": 34.57444, "lon": -111.19833, "elev": -99999}, + "KFTG": {"lat": 39.78667, "lon": -104.54528, "elev": 5497}, + "KFWS": {"lat": 32.57278, "lon": -97.30278, "elev": 683}, + "KGGW": {"lat": 48.20639, "lon": -106.62417, "elev": 2276}, + "KGJX": {"lat": 39.06222, "lon": -108.21306, "elev": 9992}, + "KGLD": {"lat": 39.36694, "lon": -101.7, "elev": 3651}, + "KGRB": {"lat": 44.49833, "lon": -88.11111, "elev": 682}, + "KGRK": {"lat": 30.72167, "lon": -97.38278, "elev": 538}, + "KGRR": {"lat": 42.89389, "lon": -85.54472, "elev": 778}, + "KGSP": {"lat": 34.88306, "lon": -82.22028, "elev": 940}, + "KGWX": {"lat": 33.89667, "lon": -88.32889, "elev": 476}, + "KGYX": {"lat": 43.89139, "lon": -70.25694, "elev": 409}, + "KHDX": {"lat": 33.07639, "lon": -106.12222, "elev": 4222}, + "KHGX": {"lat": 29.47194, "lon": -95.07889, "elev": 18}, + "KHNX": {"lat": 36.31417, "lon": -119.63111, "elev": 243}, + "KHPX": {"lat": 36.73667, "lon": -87.285, "elev": 576}, + "KHTX": {"lat": 34.93056, "lon": -86.08361, "elev": 1760}, + "KICT": {"lat": 37.65444, "lon": -97.4425, "elev": 1335}, + "KICX": {"lat": 37.59083, "lon": -112.86222, "elev": 10600}, + "KILN": {"lat": 39.42028, "lon": -83.82167, "elev": 1056}, + "KILX": {"lat": 40.15056, "lon": -89.33667, "elev": 582}, + "KIND": {"lat": 39.7075, "lon": -86.28028, "elev": 790}, + "KINX": {"lat": 36.175, "lon": -95.56444, "elev": 668}, + "KIWA": {"lat": 33.28917, "lon": -111.66917, "elev": 1353}, + "KIWX": {"lat": 41.40861, "lon": -85.7, "elev": 960}, + "KJAX": {"lat": 30.48444, "lon": -81.70194, "elev": 33}, + "KJGX": {"lat": 32.675, "lon": -83.35111, "elev": 521}, + "KJKL": {"lat": 37.59083, "lon": -83.31306, "elev": 1364}, + "KLBB": {"lat": 33.65417, "lon": -101.81361, "elev": 3259}, + "KLCH": {"lat": 30.125, "lon": -93.21583, "elev": 13}, + "KLGX": {"lat": 47.1158, "lon": -124.1069, "elev": 252}, + "KLIX": {"lat": 30.33667, "lon": -89.82528, "elev": 24}, + "KLNX": {"lat": 41.95778, "lon": -100.57583, "elev": 2970}, + "KLOT": {"lat": 41.60444, "lon": -88.08472, "elev": 663}, + "KLRX": {"lat": 40.73972, "lon": -116.80278, "elev": 6744}, + "KLSX": {"lat": 38.69889, "lon": -90.68278, "elev": 608}, + "KLTX": {"lat": 33.98917, "lon": -78.42917, "elev": 64}, + "KLVX": {"lat": 37.97528, "lon": -85.94389, "elev": 719}, + "KLWX": {"lat": 38.97628, "lon": -77.48751, "elev": -99999}, + "KLZK": {"lat": 34.83639, "lon": -92.26194, "elev": 568}, + "KMAF": {"lat": 31.94333, "lon": -102.18889, "elev": 2868}, + "KMAX": {"lat": 42.08111, "lon": -122.71611, "elev": 7513}, + "KMBX": {"lat": 48.3925, "lon": -100.86444, "elev": 1493}, + "KMHX": {"lat": 34.77583, "lon": -76.87639, "elev": 31}, + "KMKX": {"lat": 42.96778, "lon": -88.55056, "elev": 958}, + "KMLB": {"lat": 28.11306, "lon": -80.65444, "elev": 99}, + "KMOB": {"lat": 30.67944, "lon": -88.23972, "elev": 208}, + "KMPX": {"lat": 44.84889, "lon": -93.56528, "elev": 946}, + "KMQT": {"lat": 46.53111, "lon": -87.54833, "elev": 1411}, + "KMRX": {"lat": 36.16833, "lon": -83.40194, "elev": 1337}, + "KMSX": {"lat": 47.04111, "lon": -113.98611, "elev": 7855}, + "KMTX": {"lat": 41.26278, "lon": -112.44694, "elev": 6460}, + "KMUX": {"lat": 37.15528, "lon": -121.8975, "elev": 3469}, + "KMVX": {"lat": 47.52806, "lon": -97.325, "elev": 986}, + "KMXX": {"lat": 32.53667, "lon": -85.78972, "elev": 400}, + "KNKX": {"lat": 32.91889, "lon": -117.04194, "elev": 955}, + "KNQA": {"lat": 35.34472, "lon": -89.87333, "elev": 282}, + "KOAX": {"lat": 41.32028, "lon": -96.36639, "elev": 1148}, + "KOHX": {"lat": 36.24722, "lon": -86.5625, "elev": 579}, + "KOKX": {"lat": 40.86556, "lon": -72.86444, "elev": 85}, + "KOTX": {"lat": 47.68056, "lon": -117.62583, "elev": 2384}, + "KPAH": {"lat": 37.06833, "lon": -88.77194, "elev": 392}, + "KPBZ": {"lat": 40.53167, "lon": -80.21833, "elev": 1185}, + "KPDT": {"lat": 45.69056, "lon": -118.85278, "elev": 1515}, + "KPOE": {"lat": 31.15528, "lon": -92.97583, "elev": 408}, + "KPUX": {"lat": 38.45944, "lon": -104.18139, "elev": 5249}, + "KRAX": {"lat": 35.66528, "lon": -78.49, "elev": 348}, + "KRGX": {"lat": 39.75417, "lon": -119.46111, "elev": 8299}, + "KRIW": {"lat": 43.06611, "lon": -108.47667, "elev": 5568}, + "KRLX": {"lat": 38.31194, "lon": -81.72389, "elev": 1080}, + "KRTX": {"lat": 45.715, "lon": -122.96417, "elev": -99999}, + "KSFX": {"lat": 43.10583, "lon": -112.68528, "elev": 4474}, + "KSGF": {"lat": 37.23528, "lon": -93.40028, "elev": 1278}, + "KSHV": {"lat": 32.45056, "lon": -93.84111, "elev": 273}, + "KSJT": {"lat": 31.37111, "lon": -100.49222, "elev": 1890}, + "KSOX": {"lat": 33.81778, "lon": -117.635, "elev": 3027}, + "KSRX": {"lat": 35.29056, "lon": -94.36167, "elev": -99999}, + "KTBW": {"lat": 27.70528, "lon": -82.40194, "elev": 41}, + "KTFX": {"lat": 47.45972, "lon": -111.38444, "elev": 3714}, + "KTLH": {"lat": 30.3975, "lon": -84.32889, "elev": 63}, + "KTLX": {"lat": 35.33306, "lon": -97.2775, "elev": 1213}, + "KTWX": {"lat": 38.99694, "lon": -96.2325, "elev": 1367}, + "KTYX": {"lat": 43.75583, "lon": -75.68, "elev": 1846}, + "KUDX": {"lat": 44.125, "lon": -102.82944, "elev": 3016}, + "KUEX": {"lat": 40.32083, "lon": -98.44167, "elev": 1976}, + "KVAX": {"lat": 30.89, "lon": -83.00194, "elev": 178}, + "KVBX": {"lat": 34.83806, "lon": -120.39583, "elev": 1233}, + "KVNX": {"lat": 36.74083, "lon": -98.1275, "elev": 1210}, + "KVTX": {"lat": 34.41167, "lon": -119.17861, "elev": 2726}, + "KVWX": {"lat": 38.2600, "lon": -87.7247, "elev": -99999}, + "KYUX": {"lat": 32.49528, "lon": -114.65583, "elev": 174}, + "LPLA": {"lat": 38.73028, "lon": -27.32167, "elev": 3334}, + "PABC": {"lat": 60.79278, "lon": -161.87417, "elev": 162}, + "PACG": {"lat": 56.85278, "lon": -135.52917, "elev": 270}, + "PAEC": {"lat": 64.51139, "lon": -165.295, "elev": 54}, + "PAHG": {"lat": 60.725914, "lon": -151.35146, "elev": 243}, + "PAIH": {"lat": 59.46194, "lon": -146.30111, "elev": 67}, + "PAKC": {"lat": 58.67944, "lon": -156.62944, "elev": 63}, + "PAPD": {"lat": 65.03556, "lon": -147.49917, "elev": 2593}, + "PGUA": {"lat": 13.45444, "lon": 144.80833, "elev": 264}, + "PHKI": {"lat": 21.89417, "lon": -159.55222, "elev": 179}, + "PHKM": {"lat": 20.12556, "lon": -155.77778, "elev": 3812}, + "PHMO": {"lat": 21.13278, "lon": -157.18, "elev": 1363}, + "PHWA": {"lat": 19.095, "lon": -155.56889, "elev": 1370}, + "RKJK": {"lat": 35.92417, "lon": 126.62222, "elev": 78}, + "RKSG": {"lat": 36.95972, "lon": 127.01833, "elev": 52}, + "RODN": {"lat": 26.30194, "lon": 127.90972, "elev": 218}, + "TJUA": {"lat": 18.1175, "lon": -66.07861, "elev": 2794}, + "TJFK": {"lat": 40.5668, "lon": -73.8874, "elev": 112}, + "TADW": {"lat": 38.6704, "lon": -76.8446, "elev": 346}, + "TATL": {"lat": 33.6433, "lon": -84.2524, "elev": 1075}, + "TBNA": {"lat": 35.9767, "lon": -86.6618, "elev": 817}, + "TBOS": {"lat": 42.1515, "lon": -70.9302, "elev": 264}, + "TBWI": {"lat": 39.0870, "lon": -76.6276, "elev": 297}, + "TCLT": {"lat": 35.3269, "lon": -80.8772, "elev": 871}, + "TCMH": {"lat": 39.9878, "lon": -82.71, "elev": 1148}, + "TCVG": {"lat": 38.8799, "lon": -84.5737, "elev": 1053}, + "TDAL": {"lat": 32.9076, "lon": -96.9568, "elev": 622}, + "TDAY": {"lat": 39.9875, "lon": -84.1102, "elev": 1019}, + "TDCA": {"lat": 38.7474, "lon": -76.9509, "elev": 345}, + "TDEN": {"lat": 39.7256, "lon": -104.5431, "elev": 5701}, + "TDFW": {"lat": 33.0396, "lon": -96.8974, "elev": 585}, + "TDTW": {"lat": 42.0710, "lon": -83.4704, "elev": 772}, + "TEWR": {"lat": 40.5880, "lon": -74.2503, "elev": 136}, + "TFLL": {"lat": 26.1263, "lon": -80.3478, "elev": 120}, + "THOU": {"lat": 29.5328, "lon": -95.2444, "elev": 117}, + "TIAD": {"lat": 39.0675, "lon": -77.5012, "elev": 473}, + "TIAH": {"lat": 30.0297, "lon": -95.5708, "elev": 253}, + "TICH": {"lat": 37.4069, "lon": -97.4764, "elev": 1351}, + "TIDS": {"lat": 39.5978, "lon": -86.4085, "elev": 847}, + "TLAS": {"lat": 36.1292, "lon": -115.0147, "elev": 2058}, + "TLVE": {"lat": 41.2805, "lon": -81.9659, "elev": 931}, + "TMCI": {"lat": 39.4488, "lon": -94.7396, "elev": 1090}, + "TMCO": {"lat": 28.2584, "lon": -81.3133, "elev": 169}, + "TMDW": {"lat": 41.69, "lon": -87.8034, "elev": 763}, + "TMEM": {"lat": 34.8867, "lon": -90.0007, "elev": 483}, + "TMIA": {"lat": 25.7555, "lon": -80.4932, "elev": 125}, + "TMKE": {"lat": 42.7619, "lon": -87.9994, "elev": 933}, + "TMSP": {"lat": 44.8197, "lon": -92.9392, "elev": 1121}, + "TMSY": {"lat": 29.9385, "lon": -90.3811, "elev": 99}, + "TOKC": {"lat": 35.2474, "lon": -97.5395, "elev": 1308}, + "TORD": {"lat": 41.7712, "lon": -87.8363, "elev": 744}, + "TPBI": {"lat": 26.6572, "lon": -80.2586, "elev": 133}, + "TPHL": {"lat": 39.9084, "lon": -75.0426, "elev": 153}, + "TPHX": {"lat": 33.3678, "lon": -112.1580, "elev": 1089}, + "TPIT": {"lat": 40.4641, "lon": -80.4697, "elev": 1386}, + "TRDU": {"lat": 35.9898, "lon": -78.6787, "elev": 515}, + "TSDF": {"lat": 38.0109, "lon": -85.5995, "elev": 731}, + "TSJU": {"lat": 18.4313, "lon": -66.1722, "elev": 157}, + "TSLC": {"lat": 40.9341, "lon": -111.9214, "elev": 4295}, + "TSTL": {"lat": 38.7668, "lon": -90.4698, "elev": 647}, + "TTPA": {"lat": 27.8196, "lon": -82.5179, "elev": 93}, + "TTUL": {"lat": 36.0236, "lon": -95.8175, "elev": 823}, +} diff --git a/xradar/io/backends/nexrad_interpolate.pyx b/xradar/io/backends/nexrad_interpolate.pyx new file mode 100644 index 0000000..2930d29 --- /dev/null +++ b/xradar/io/backends/nexrad_interpolate.pyx @@ -0,0 +1,110 @@ +""" +Interpolation of NEXRAD moments from 1000 meter to 250 meter gate spacing. + +""" + +def _fast_interpolate_scan_4( + float[:, :] data, float[:] scratch_ray, float fill_value, + int start, int end, int moment_ngates, int linear_interp): + """ Interpolate a single NEXRAD moment scan from 1000 m to 250 m. """ + # This interpolation scheme is only valid for NEXRAD data where a 4:1 + # (1000 m : 250 m) interpolation is needed. + # + # The scheme here performs a linear interpolation between pairs of gates + # in a ray when the both of the gates are not masked (below threshold). + # When one of the gates is masked the interpolation changes to a nearest + # neighbor interpolation. Nearest neighbor is also performed at the end + # points until the new range bin would be centered beyond half of the range + # spacing of the original range. + # + # Nearest neighbor interpolation is performed when linear_interp is False, + # this is equivalent to repeating each gate four times in each ray. + # + # No transformation of the raw data is performed prior to interpolation, so + # reflectivity will be interpolated in dB units, velocity in m/s, etc, + # this may not be the best method for interpolation. + # + # This method was adapted from Radx and Py-ART + cdef int ray_num, i, interp_ngates + cdef float gate_val, next_val, delta + + interp_ngates = 4 * moment_ngates # number of gates interpolated + + for ray_num in range(start, end+1): + + # repeat each gate value 4 times + for i in range(moment_ngates): + gate_val = data[ray_num, i] + scratch_ray[i*4 + 0] = gate_val + scratch_ray[i*4 + 1] = gate_val + scratch_ray[i*4 + 2] = gate_val + scratch_ray[i*4 + 3] = gate_val + + if linear_interp: + # linear interpolate + for i in range(2, interp_ngates - 4, 4): + gate_val = scratch_ray[i] + next_val = scratch_ray[i+4] + if gate_val == fill_value or next_val == fill_value: + continue + delta = (next_val - gate_val) / 4. + scratch_ray[i+0] = gate_val + delta * 0.5 + scratch_ray[i+1] = gate_val + delta * 1.5 + scratch_ray[i+2] = gate_val + delta * 2.5 + scratch_ray[i+3] = gate_val + delta * 3.5 + + for i in range(interp_ngates): + data[ray_num, i] = scratch_ray[i] + + +def _fast_interpolate_scan_2( + float[:, :] data, float[:] scratch_ray, float fill_value, + int start, int end, int moment_ngates, int linear_interp): + """ Interpolate a single NEXRAD moment scan from 300 m to 150 m. """ + # This interpolation scheme is only valid for NEXRAD TWDR data where a 2:1 + # (300 m : 150 m) interpolation is needed. + # + # The scheme here performs a linear interpolation between pairs of gates + # in a ray when the both of the gates are not masked (below threshold). + # When one of the gates is masked the interpolation changes to a nearest + # neighbor interpolation. Nearest neighbor is also performed at the end + # points until the new range bin would be centered beyond half of the range + # spacing of the original range. + # + # Nearest neighbor interpolation is performed when linear_interp is False, + # this is equivalent to repeating each gate four times in each ray. + # + # No transformation of the raw data is performed prior to interpolation, so + # reflectivity will be interpolated in dB units, velocity in m/s, etc, + # this may not be the best method for interpolation. + # + # This method was adapted from Radx + cdef int ray_num, i, interp_ngates + cdef float gate_val, next_val, delta + + interp_ngates = 2 * moment_ngates - 1 # number of gates interpolated + + for ray_num in range(start, end+1): + + # repeat each gate value 4 times + for i in range(moment_ngates): + gate_val = data[ray_num, i] + if i == moment_ngates - 1: + scratch_ray[i*2 + 0] = gate_val + else: + scratch_ray[i*2 + 0] = gate_val + scratch_ray[i*2 + 1] = gate_val + + if linear_interp: + # linear interpolate + for i in range(1, interp_ngates - 2, 2): + gate_val = scratch_ray[i] + next_val = scratch_ray[i+2] + if gate_val == fill_value or next_val == fill_value: + continue + delta = (next_val - gate_val) / 2. + scratch_ray[i+0] = gate_val + delta * 0.5 + scratch_ray[i+1] = gate_val + delta * 1.5 + + for i in range(interp_ngates): + data[ray_num, i] = scratch_ray[i] diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py new file mode 100644 index 0000000..e44bc14 --- /dev/null +++ b/xradar/io/backends/nexrad_level2.py @@ -0,0 +1,960 @@ +""" +Functions for reading NEXRAD level 2 files. +""" + + +# This file has been adapted from is Py-ART, the Python ARM Radar Toolkit +# https://github.com/ARM-DOE/pyart + + +import bz2 +import struct +import warnings +from datetime import datetime, timedelta + +import numpy as np + + +class _NEXRADLevel2StagedField: + """ + A class to facilitate on demand loading of field data from a Level 2 file. + """ + + def __init__(self, nfile, moment, max_ngates, scans): + """initialize.""" + self.nfile = nfile + self.moment = moment + self.max_ngates = max_ngates + self.scans = scans + + def __call__(self): + """Return the array containing the field data.""" + return self.nfile.get_data(self.moment, self.max_ngates, scans=self.scans) + + +class NEXRADLevel2File: + """ + Class for accessing data in a NEXRAD (WSR-88D) Level II file. + + NEXRAD Level II files [1]_, also know as NEXRAD Archive Level II or + WSR-88D Archive level 2, are available from the NOAA National Climate Data + Center [2]_ as well as on the UCAR THREDDS Data Server [3]_. Files with + uncompressed messages and compressed messages are supported. This class + supports reading both "message 31" and "message 1" type files. + + Parameters + ---------- + filename : str + Filename of Archive II file to read. + + Attributes + ---------- + radial_records : list + Radial (1 or 31) messages in the file. + nscans : int + Number of scans in the file. + scan_msgs : list of arrays + Each element specifies the indices of the message in the + radial_records attribute which belong to a given scan. + volume_header : dict + Volume header. + vcp : dict + VCP information dictionary. + _records : list + A list of all records (message) in the file. + _fh : file-like + File like object from which data is read. + _msg_type : '31' or '1': + Type of radial messages in file. + + References + ---------- + .. [1] http://www.roc.noaa.gov/WSR88D/Level_II/Level2Info.aspx + .. [2] http://www.ncdc.noaa.gov/ + .. [3] http://thredds.ucar.edu/thredds/catalog.html + + """ + + def __init__(self, filename): + """initalize the object.""" + # read in the volume header and compression_record + if hasattr(filename, "read"): + fh = filename + else: + fh = open(filename, "rb") + size = _structure_size(VOLUME_HEADER) + self.volume_header = _unpack_structure(fh.read(size), VOLUME_HEADER) + compression_record = fh.read(COMPRESSION_RECORD_SIZE) + + # read the records in the file, decompressing as needed + compression_slice = slice(CONTROL_WORD_SIZE, CONTROL_WORD_SIZE + 2) + compression_or_ctm_info = compression_record[compression_slice] + if compression_or_ctm_info == b"BZ": + buf = _decompress_records(fh) + # The 12-byte compression record previously held the Channel Terminal + # Manager (CTM) information. Bytes 4 through 6 contain the size of the + # record (2432) as a big endian unsigned short, which is encoded as + # b'\t\x80' == struct.pack('>H', 2432). + # Newer files zero out this section. + elif compression_or_ctm_info in (b"\x00\x00", b"\t\x80"): + buf = fh.read() + else: + raise OSError("unknown compression record") + self._fh = fh + + # read the records from the buffer + self._records = [] + buf_length = len(buf) + pos = 0 + while pos < buf_length: + pos, dic = _get_record_from_buf(buf, pos) + self._records.append(dic) + + # pull out radial records (1 or 31) which contain the moment data. + self.radial_records = [r for r in self._records if r["header"]["type"] == 31] + self._msg_type = "31" + if len(self.radial_records) == 0: + self.radial_records = [r for r in self._records if r["header"]["type"] == 1] + self._msg_type = "1" + if len(self.radial_records) == 0: + raise ValueError("No MSG31 records found, cannot read file") + elev_nums = np.array( + [m["msg_header"]["elevation_number"] for m in self.radial_records] + ) + self.scan_msgs = [ + np.where(elev_nums == i + 1)[0] for i in range(elev_nums.max()) + ] + self.nscans = len(self.scan_msgs) + + # pull out the vcp record + msg_5 = [r for r in self._records if r["header"]["type"] == 5] + + if len(msg_5): + self.vcp = msg_5[0] + else: + # There is no VCP Data.. This is uber dodgy + warnings.warn( + "No MSG5 detected. Setting to meaningless data. " + "Rethink your life choices and be ready for errors." + "Specifically fixed angle data will be missing" + ) + + self.vcp = None + return + + def close(self): + """Close the file.""" + self._fh.close() + + def location(self): + """ + Find the location of the radar. + + Returns all zeros if location is not available. + + Returns + ------- + latitude : float + Latitude of the radar in degrees. + longitude : float + Longitude of the radar in degrees. + height : int + Height of radar and feedhorn in meters above mean sea level. + + """ + if self._msg_type == "31": + dic = self.radial_records[0]["VOL"] + height = dic["height"] + dic["feedhorn_height"] + return dic["lat"], dic["lon"], height + else: + return 0.0, 0.0, 0.0 + + def scan_info(self, scans=None): + """ + Return a list of dictionaries with scan information. + + Parameters + ---------- + scans : list ot None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + scan_info : list, optional + A list of the scan performed with a dictionary with keys + 'moments', 'ngates', 'nrays', 'first_gate' and 'gate_spacing' + for each scan. The 'moments', 'ngates', 'first_gate', and + 'gate_spacing' keys are lists of the NEXRAD moments and gate + information for that moment collected during the specific scan. + The 'nrays' key provides the number of radials collected in the + given scan. + + """ + info = [] + + if scans is None: + scans = range(self.nscans) + for scan in scans: + nrays = self.get_nrays(scan) + if nrays < 2: + self.nscans -= 1 + continue + msg31_number = self.scan_msgs[scan][0] + msg = self.radial_records[msg31_number] + nexrad_moments = ["REF", "VEL", "SW", "ZDR", "PHI", "RHO", "CFP"] + moments = [f for f in nexrad_moments if f in msg] + ngates = [msg[f]["ngates"] for f in moments] + gate_spacing = [msg[f]["gate_spacing"] for f in moments] + first_gate = [msg[f]["first_gate"] for f in moments] + info.append( + { + "nrays": nrays, + "ngates": ngates, + "gate_spacing": gate_spacing, + "first_gate": first_gate, + "moments": moments, + } + ) + return info + + def get_vcp_pattern(self): + """ + Return the numerical volume coverage pattern (VCP) or None if unknown. + """ + if self.vcp is None: + return None + else: + return self.vcp["msg5_header"]["pattern_number"] + + def get_nrays(self, scan): + """ + Return the number of rays in a given scan. + + Parameters + ---------- + scan : int + Scan of interest (0 based). + + Returns + ------- + nrays : int + Number of rays (radials) in the scan. + + """ + return len(self.scan_msgs[scan]) + + def get_range(self, scan_num, moment): + """ + Return an array of gate ranges for a given scan and moment. + + Parameters + ---------- + scan_num : int + Scan number (0 based). + moment : 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'RHO', or 'CFP' + Moment of interest. + + Returns + ------- + range : ndarray + Range in meters from the antenna to the center of gate (bin). + + """ + dic = self.radial_records[self.scan_msgs[scan_num][0]][moment] + ngates = dic["ngates"] + first_gate = dic["first_gate"] + gate_spacing = dic["gate_spacing"] + return np.arange(ngates) * gate_spacing + first_gate + + # helper functions for looping over scans + def _msg_nums(self, scans): + """Find the all message number for a list of scans.""" + return np.concatenate([self.scan_msgs[i] for i in scans]) + + def _radial_array(self, scans, key): + """ + Return an array of radial header elements for all rays in scans. + """ + msg_nums = self._msg_nums(scans) + temp = [self.radial_records[i]["msg_header"][key] for i in msg_nums] + return np.array(temp) + + def _radial_sub_array(self, scans, key): + """ + Return an array of RAD or msg_header elements for all rays in scans. + """ + msg_nums = self._msg_nums(scans) + if self._msg_type == "31": + tmp = [self.radial_records[i]["RAD"][key] for i in msg_nums] + else: + tmp = [self.radial_records[i]["msg_header"][key] for i in msg_nums] + return np.array(tmp) + + def get_times(self, scans=None): + """ + Retrieve the times at which the rays were collected. + + Parameters + ---------- + scans : list or None + Scans (0-based) to retrieve ray (radial) collection times from. + None (the default) will return the times for all scans in the + volume. + + Returns + ------- + time_start : Datetime + Initial time. + time : ndarray + Offset in seconds from the initial time at which the rays + in the requested scans were collected. + + """ + if scans is None: + scans = range(self.nscans) + days = self._radial_array(scans, "collect_date") + secs = self._radial_array(scans, "collect_ms") / 1000.0 + offset = timedelta(days=int(days[0]) - 1, seconds=int(secs[0])) + time_start = datetime(1970, 1, 1) + offset + time = secs - int(secs[0]) + (days - days[0]) * 86400 + return time_start, time + + def get_azimuth_angles(self, scans=None): + """ + Retrieve the azimuth angles of all rays in the requested scans. + + Parameters + ---------- + scans : list ot None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + angles : ndarray + Azimuth angles in degress for all rays in the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "1": + scale = 180 / (4096 * 8.0) + else: + scale = 1.0 + return self._radial_array(scans, "azimuth_angle") * scale + + def get_elevation_angles(self, scans=None): + """ + Retrieve the elevation angles of all rays in the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which ray (radial) azimuth angles will be + retrieved. None (the default) will return the angles for + all scans in the volume. + + Returns + ------- + angles : ndarray + Elevation angles in degress for all rays in the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "1": + scale = 180 / (4096 * 8.0) + else: + scale = 1.0 + return self._radial_array(scans, "elevation_angle") * scale + + def get_target_angles(self, scans=None): + """ + Retrieve the target elevation angle of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the target elevation angles will be + retrieved. None (the default) will return the angles for all + scans in the volume. + + Returns + ------- + angles : ndarray + Target elevation angles in degress for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + if self._msg_type == "31": + if self.vcp is not None: + cut_parameters = self.vcp["cut_parameters"] + else: + cut_parameters = [{"elevation_angle": 0.0}] * self.nscans + scale = 360.0 / 65536.0 + return np.array( + [cut_parameters[i]["elevation_angle"] * scale for i in scans], + dtype="float32", + ) + else: + scale = 180 / (4096 * 8.0) + msgs = [self.radial_records[self.scan_msgs[i][0]] for i in scans] + return np.round( + np.array( + [m["msg_header"]["elevation_angle"] * scale for m in msgs], + dtype="float32", + ), + 1, + ) + + def get_nyquist_vel(self, scans=None): + """ + Retrieve the Nyquist velocities of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the Nyquist velocities will be + retrieved. None (the default) will return the velocities for all + scans in the volume. + + Returns + ------- + velocities : ndarray + Nyquist velocities (in m/s) for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + return self._radial_sub_array(scans, "nyquist_vel") * 0.01 + + def get_unambigous_range(self, scans=None): + """ + Retrieve the unambiguous range of the requested scans. + + Parameters + ---------- + scans : list or None + Scans (0 based) for which the unambiguous range will be retrieved. + None (the default) will return the range for all scans in the + volume. + + Returns + ------- + unambiguous_range : ndarray + Unambiguous range (in meters) for the requested scans. + + """ + if scans is None: + scans = range(self.nscans) + # unambiguous range is stored in tenths of km, x100 for meters + return self._radial_sub_array(scans, "unambig_range") * 100.0 + + def get_data(self, moment, max_ngates, scans=None, raw_data=False): + """ + Retrieve moment data for a given set of scans. + + Masked points indicate that the data was not collected, below + threshold or is range folded. + + Parameters + ---------- + moment : 'REF', 'VEL', 'SW', 'ZDR', 'PHI', 'RHO', or 'CFP' + Moment for which to to retrieve data. + max_ngates : int + Maximum number of gates (bins) in any ray. + requested. + raw_data : bool + True to return the raw data, False to perform masking as well as + applying the appropiate scale and offset to the data. When + raw_data is True values of 1 in the data likely indicate that + the gate was not present in the sweep, in some cases in will + indicate range folded data. + scans : list or None. + Scans to retrieve data from (0 based). None (the default) will + get the data for all scans in the volume. + + Returns + ------- + data : ndarray + + """ + if scans is None: + scans = range(self.nscans) + + # determine the number of rays + msg_nums = self._msg_nums(scans) + nrays = len(msg_nums) + # extract the data + set_datatype = False + data = np.ones((nrays, max_ngates), ">B") + for i, msg_num in enumerate(msg_nums): + msg = self.radial_records[msg_num] + if moment not in msg.keys(): + continue + if not set_datatype: + data = data.astype(">" + _bits_to_code(msg, moment)) + set_datatype = True + + ngates = min(msg[moment]["ngates"], max_ngates, len(msg[moment]["data"])) + data[i, :ngates] = msg[moment]["data"][:ngates] + # return raw data if requested + if raw_data: + return data + + # mask, scan and offset, assume that the offset and scale + # are the same in all scans/gates + for scan in scans: # find a scan which contains the moment + msg_num = self.scan_msgs[scan][0] + msg = self.radial_records[msg_num] + if moment in msg.keys(): + offset = np.float32(msg[moment]["offset"]) + scale = np.float32(msg[moment]["scale"]) + mask = data <= 1 + scaled_data = (data - offset) / scale + return np.ma.array(scaled_data, mask=mask) + + # moment is not present in any scan, mask all values + return np.ma.masked_less_equal(data, 1) + + +def _bits_to_code(msg, moment): + """ + Convert number of bits to the proper code for unpacking. + Based on the code found in MetPy: + https://github.com/Unidata/MetPy/blob/40d5c12ab341a449c9398508bd41 + d010165f9eeb/src/metpy/io/_tools.py#L313-L321 + """ + if msg["header"]["type"] == 1: + word_size = msg[moment]["data"].dtype + if word_size == "uint16": + return "H" + elif word_size == "uint8": + return "B" + else: + warnings.warn(('Unsupported bit size: %s. Returning "B"', word_size)) + return "B" + + elif msg["header"]["type"] == 31: + word_size = msg[moment]["word_size"] + if word_size == 16: + return "H" + elif word_size == 8: + return "B" + else: + warnings.warn(('Unsupported bit size: %s. Returning "B"', word_size)) + return "B" + else: + raise TypeError("Unsupported msg type %s", msg["header"]["type"]) + + +def _decompress_records(file_handler): + """ + Decompressed the records from an BZ2 compressed Archive 2 file. + """ + file_handler.seek(0) + cbuf = file_handler.read() # read all data from the file + decompressor = bz2.BZ2Decompressor() + skip = _structure_size(VOLUME_HEADER) + CONTROL_WORD_SIZE + buf = decompressor.decompress(cbuf[skip:]) + while len(decompressor.unused_data): + cbuf = decompressor.unused_data + decompressor = bz2.BZ2Decompressor() + buf += decompressor.decompress(cbuf[CONTROL_WORD_SIZE:]) + + return buf[COMPRESSION_RECORD_SIZE:] + + +def _get_record_from_buf(buf, pos): + """Retrieve and unpack a NEXRAD record from a buffer.""" + dic = {"header": _unpack_from_buf(buf, pos, MSG_HEADER)} + msg_type = dic["header"]["type"] + + if msg_type == 31: + new_pos = _get_msg31_from_buf(buf, pos, dic) + elif msg_type == 5: + # Sometimes we encounter incomplete buffers + try: + new_pos = _get_msg5_from_buf(buf, pos, dic) + except struct.error: + warnings.warn( + "Encountered incomplete MSG5. File may be corrupt.", RuntimeWarning + ) + new_pos = pos + RECORD_SIZE + elif msg_type == 29: + new_pos = _get_msg29_from_buf(pos, dic) + warnings.warn("Message 29 encountered, not parsing.", RuntimeWarning) + elif msg_type == 1: + new_pos = _get_msg1_from_buf(buf, pos, dic) + else: # not message 31 or 1, no decoding performed + new_pos = pos + RECORD_SIZE + + return new_pos, dic + + +def _get_msg29_from_buf(pos, dic): + msg_size = dic["header"]["size"] + if msg_size == 65535: + msg_size = dic["header"]["segments"] << 16 | dic["header"]["seg_num"] + msg_header_size = _structure_size(MSG_HEADER) + new_pos = pos + msg_header_size + msg_size + return new_pos + + +def _get_msg31_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG31 record from a buffer.""" + msg_size = dic["header"]["size"] * 2 - 4 + msg_header_size = _structure_size(MSG_HEADER) + new_pos = pos + msg_header_size + msg_size + mbuf = buf[pos + msg_header_size : new_pos] + msg_31_header = _unpack_from_buf(mbuf, 0, MSG_31) + block_pointers = [ + v for k, v in msg_31_header.items() if k.startswith("block_pointer") and v > 0 + ] + for block_pointer in block_pointers: + block_name, block_dic = _get_msg31_data_block(mbuf, block_pointer) + dic[block_name] = block_dic + + dic["msg_header"] = msg_31_header + return new_pos + + +def _get_msg31_data_block(buf, ptr): + """Unpack a msg_31 data block into a dictionary.""" + block_name = buf[ptr + 1 : ptr + 4].decode("ascii").strip() + + if block_name == "VOL": + dic = _unpack_from_buf(buf, ptr, VOLUME_DATA_BLOCK) + elif block_name == "ELV": + dic = _unpack_from_buf(buf, ptr, ELEVATION_DATA_BLOCK) + elif block_name == "RAD": + dic = _unpack_from_buf(buf, ptr, RADIAL_DATA_BLOCK) + elif block_name in ["REF", "VEL", "SW", "ZDR", "PHI", "RHO", "CFP"]: + dic = _unpack_from_buf(buf, ptr, GENERIC_DATA_BLOCK) + ngates = dic["ngates"] + ptr2 = ptr + _structure_size(GENERIC_DATA_BLOCK) + if dic["word_size"] == 16: + data = np.frombuffer(buf[ptr2 : ptr2 + ngates * 2], ">u2") + elif dic["word_size"] == 8: + data = np.frombuffer(buf[ptr2 : ptr2 + ngates], ">u1") + else: + warnings.warn( + 'Unsupported bit size: %s. Returning array dtype "B"', dic["word_size"] + ) + dic["data"] = data + else: + dic = {} + return block_name, dic + + +def _get_msg1_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG1 record from a buffer.""" + msg_header_size = _structure_size(MSG_HEADER) + msg1_header = _unpack_from_buf(buf, pos + msg_header_size, MSG_1) + dic["msg_header"] = msg1_header + + sur_nbins = int(msg1_header["sur_nbins"]) + doppler_nbins = int(msg1_header["doppler_nbins"]) + + sur_step = int(msg1_header["sur_range_step"]) + doppler_step = int(msg1_header["doppler_range_step"]) + + sur_first = int(msg1_header["sur_range_first"]) + doppler_first = int(msg1_header["doppler_range_first"]) + if doppler_first > 2**15: + doppler_first = doppler_first - 2**16 + + if msg1_header["sur_pointer"]: + offset = pos + msg_header_size + msg1_header["sur_pointer"] + data = np.frombuffer(buf[offset : offset + sur_nbins], ">u1") + dic["REF"] = { + "ngates": sur_nbins, + "gate_spacing": sur_step, + "first_gate": sur_first, + "data": data, + "scale": 2.0, + "offset": 66.0, + } + if msg1_header["vel_pointer"]: + offset = pos + msg_header_size + msg1_header["vel_pointer"] + data = np.frombuffer(buf[offset : offset + doppler_nbins], ">u1") + dic["VEL"] = { + "ngates": doppler_nbins, + "gate_spacing": doppler_step, + "first_gate": doppler_first, + "data": data, + "scale": 2.0, + "offset": 129.0, + } + if msg1_header["doppler_resolution"] == 4: + # 1 m/s resolution velocity, offset remains 129. + dic["VEL"]["scale"] = 1.0 + if msg1_header["width_pointer"]: + offset = pos + msg_header_size + msg1_header["width_pointer"] + data = np.frombuffer(buf[offset : offset + doppler_nbins], ">u1") + dic["SW"] = { + "ngates": doppler_nbins, + "gate_spacing": doppler_step, + "first_gate": doppler_first, + "data": data, + "scale": 2.0, + "offset": 129.0, + } + return pos + RECORD_SIZE + + +def _get_msg5_from_buf(buf, pos, dic): + """Retrieve and unpack a MSG1 record from a buffer.""" + msg_header_size = _structure_size(MSG_HEADER) + msg5_header_size = _structure_size(MSG_5) + msg5_elev_size = _structure_size(MSG_5_ELEV) + + dic["msg5_header"] = _unpack_from_buf(buf, pos + msg_header_size, MSG_5) + dic["cut_parameters"] = [] + for i in range(dic["msg5_header"]["num_cuts"]): + pos2 = pos + msg_header_size + msg5_header_size + msg5_elev_size * i + dic["cut_parameters"].append(_unpack_from_buf(buf, pos2, MSG_5_ELEV)) + return pos + RECORD_SIZE + + +def _structure_size(structure): + """Find the size of a structure in bytes.""" + return struct.calcsize(">" + "".join([i[1] for i in structure])) + + +def _unpack_from_buf(buf, pos, structure): + """Unpack a structure from a buffer.""" + size = _structure_size(structure) + return _unpack_structure(buf[pos : pos + size], structure) + + +def _unpack_structure(string, structure): + """Unpack a structure from a string.""" + fmt = ">" + "".join([i[1] for i in structure]) # NEXRAD is big-endian + lst = struct.unpack(fmt, string) + return dict(zip([i[0] for i in structure], lst)) + + +# NEXRAD Level II file structures and sizes +# The deails on these structures are documented in: +# "Interface Control Document for the Achive II/User" RPG Build 12.0 +# Document Number 2620010E +# and +# "Interface Control Document for the RDA/RPG" Open Build 13.0 +# Document Number 2620002M +# Tables and page number refer to those in the second document unless +# otherwise noted. +RECORD_SIZE = 2432 +COMPRESSION_RECORD_SIZE = 12 +CONTROL_WORD_SIZE = 4 + +# format of structure elements +# section 3.2.1, page 3-2 +CODE1 = "B" +CODE2 = "H" +INT1 = "B" +INT2 = "H" +INT4 = "I" +REAL4 = "f" +REAL8 = "d" +SINT1 = "b" +SINT2 = "h" +SINT4 = "i" + +# Figure 1 in Interface Control Document for the Archive II/User +# page 7-2 +VOLUME_HEADER = ( + ("tape", "9s"), + ("extension", "3s"), + ("date", "I"), + ("time", "I"), + ("icao", "4s"), +) + +# Table II Message Header Data +# page 3-7 +MSG_HEADER = ( + ("size", INT2), # size of data, no including header + ("channels", INT1), + ("type", INT1), + ("seq_id", INT2), + ("date", INT2), + ("ms", INT4), + ("segments", INT2), + ("seg_num", INT2), +) + +# Table XVII Digital Radar Generic Format Blocks (Message Type 31) +# pages 3-87 to 3-89 +MSG_31 = ( + ("id", "4s"), # 0-3 + ("collect_ms", INT4), # 4-7 + ("collect_date", INT2), # 8-9 + ("azimuth_number", INT2), # 10-11 + ("azimuth_angle", REAL4), # 12-15 + ("compress_flag", CODE1), # 16 + ("spare_0", INT1), # 17 + ("radial_length", INT2), # 18-19 + ("azimuth_resolution", CODE1), # 20 + ("radial_spacing", CODE1), # 21 + ("elevation_number", INT1), # 22 + ("cut_sector", INT1), # 23 + ("elevation_angle", REAL4), # 24-27 + ("radial_blanking", CODE1), # 28 + ("azimuth_mode", SINT1), # 29 + ("block_count", INT2), # 30-31 + ("block_pointer_1", INT4), # 32-35 Volume Data Constant XVII-E + ("block_pointer_2", INT4), # 36-39 Elevation Data Constant XVII-F + ("block_pointer_3", INT4), # 40-43 Radial Data Constant XVII-H + ("block_pointer_4", INT4), # 44-47 Moment "REF" XVII-{B/I} + ("block_pointer_5", INT4), # 48-51 Moment "VEL" + ("block_pointer_6", INT4), # 52-55 Moment "SW" + ("block_pointer_7", INT4), # 56-59 Moment "ZDR" + ("block_pointer_8", INT4), # 60-63 Moment "PHI" + ("block_pointer_9", INT4), # 64-67 Moment "RHO" + ("block_pointer_10", INT4), # Moment "CFP" +) + + +# Table III Digital Radar Data (Message Type 1) +# pages 3-7 to +MSG_1 = ( + ("collect_ms", INT4), # 0-3 + ("collect_date", INT2), # 4-5 + ("unambig_range", SINT2), # 6-7 + ("azimuth_angle", CODE2), # 8-9 + ("azimuth_number", INT2), # 10-11 + ("radial_status", CODE2), # 12-13 + ("elevation_angle", INT2), # 14-15 + ("elevation_number", INT2), # 16-17 + ("sur_range_first", CODE2), # 18-19 + ("doppler_range_first", CODE2), # 20-21 + ("sur_range_step", CODE2), # 22-23 + ("doppler_range_step", CODE2), # 24-25 + ("sur_nbins", INT2), # 26-27 + ("doppler_nbins", INT2), # 28-29 + ("cut_sector_num", INT2), # 30-31 + ("calib_const", REAL4), # 32-35 + ("sur_pointer", INT2), # 36-37 + ("vel_pointer", INT2), # 38-39 + ("width_pointer", INT2), # 40-41 + ("doppler_resolution", CODE2), # 42-43 + ("vcp", INT2), # 44-45 + ("spare_1", "8s"), # 46-53 + ("spare_2", "2s"), # 54-55 + ("spare_3", "2s"), # 56-57 + ("spare_4", "2s"), # 58-59 + ("nyquist_vel", SINT2), # 60-61 + ("atmos_attenuation", SINT2), # 62-63 + ("threshold", SINT2), # 64-65 + ("spot_blank_status", INT2), # 66-67 + ("spare_5", "32s"), # 68-99 + # 100+ reflectivity, velocity and/or spectral width data, CODE1 +) + +# Table XI Volume Coverage Pattern Data (Message Type 5 & 7) +# pages 3-51 to 3-54 +MSG_5 = ( + ("msg_size", INT2), + ("pattern_type", CODE2), + ("pattern_number", INT2), + ("num_cuts", INT2), + ("clutter_map_group", INT2), + ("doppler_vel_res", CODE1), # 2: 0.5 degrees, 4: 1.0 degrees + ("pulse_width", CODE1), # 2: short, 4: long + ("spare", "10s"), # halfwords 7-11 (10 bytes, 5 halfwords) +) + +MSG_5_ELEV = ( + ("elevation_angle", CODE2), # scaled by 360/65536 for value in degrees. + ("channel_config", CODE1), + ("waveform_type", CODE1), + ("super_resolution", CODE1), + ("prf_number", INT1), + ("prf_pulse_count", INT2), + ("azimuth_rate", CODE2), + ("ref_thresh", SINT2), + ("vel_thresh", SINT2), + ("sw_thresh", SINT2), + ("zdr_thres", SINT2), + ("phi_thres", SINT2), + ("rho_thres", SINT2), + ("edge_angle_1", CODE2), + ("dop_prf_num_1", INT2), + ("dop_prf_pulse_count_1", INT2), + ("spare_1", "2s"), + ("edge_angle_2", CODE2), + ("dop_prf_num_2", INT2), + ("dop_prf_pulse_count_2", INT2), + ("spare_2", "2s"), + ("edge_angle_3", CODE2), + ("dop_prf_num_3", INT2), + ("dop_prf_pulse_count_3", INT2), + ("spare_3", "2s"), +) + +# Table XVII-B Data Block (Descriptor of Generic Data Moment Type) +# pages 3-90 and 3-91 +GENERIC_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), # VEL, REF, SW, RHO, PHI, ZDR + ("reserved", INT4), + ("ngates", INT2), + ("first_gate", SINT2), + ("gate_spacing", SINT2), + ("thresh", SINT2), + ("snr_thres", SINT2), + ("flags", CODE1), + ("word_size", INT1), + ("scale", REAL4), + ("offset", REAL4), + # then data +) + +# Table XVII-E Data Block (Volume Data Constant Type) +# page 3-92 +VOLUME_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("version_major", INT1), + ("version_minor", INT1), + ("lat", REAL4), + ("lon", REAL4), + ("height", SINT2), + ("feedhorn_height", INT2), + ("refl_calib", REAL4), + ("power_h", REAL4), + ("power_v", REAL4), + ("diff_refl_calib", REAL4), + ("init_phase", REAL4), + ("vcp", INT2), + ("spare", "2s"), +) + +# Table XVII-F Data Block (Elevation Data Constant Type) +# page 3-93 +ELEVATION_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("atmos", SINT2), + ("refl_calib", REAL4), +) + +# Table XVII-H Data Block (Radial Data Constant Type) +# pages 3-93 +RADIAL_DATA_BLOCK = ( + ("block_type", "1s"), + ("data_name", "3s"), + ("lrtup", INT2), + ("unambig_range", SINT2), + ("noise_h", REAL4), + ("noise_v", REAL4), + ("nyquist_vel", SINT2), + ("spare", "2s"), +) From f7e42d228f4e74f69e3eb35657a70189293c15ae Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 8 Jan 2024 09:49:07 -0600 Subject: [PATCH 02/48] FIX: Fix incorrect module registration --- setup.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 setup.py diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..db3bcee --- /dev/null +++ b/setup.py @@ -0,0 +1,27 @@ +import numpy +from Cython.Build import cythonize +from Cython.Compiler import Options +from setuptools import Extension, setup + +# These are optional +Options.docstrings = True +Options.annotate = False + +# Modules to be compiled and include_dirs when necessary +extensions = [ + Extension( + "xradar.io.backends.nexrad_interpolate", + sources=["xradar/io/backends/nexrad_interpolate.pyx"], + include_dirs=[numpy.get_include()], + ), +] + + +# This is the function that is executed +setup( + name="xradar", # Required + # external to be compiled + ext_modules=cythonize( + extensions, compiler_directives={"language_level": "3", "cpow": True} + ), +) From 34b52b02956c73d187816b5cc87dc5a84eeaba1f Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 8 Jan 2024 09:53:50 -0600 Subject: [PATCH 03/48] FIX: Fix the requirements for the package --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 788ca1f..9d24919 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,6 +45,8 @@ requires = [ "setuptools>=45", "wheel", "setuptools_scm[toml]>=7.0", + "cython", + "numpy" ] build-backend = "setuptools.build_meta" From 9330b173d069596481a4bf13c10034ff2cf3e35b Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 18 Jan 2024 15:58:52 -0600 Subject: [PATCH 04/48] ADD: Add new xarray dataset builder --- xradar/io/backends/nexrad_archive.py | 194 +++++++++++++++++++++++++-- 1 file changed, 182 insertions(+), 12 deletions(-) diff --git a/xradar/io/backends/nexrad_archive.py b/xradar/io/backends/nexrad_archive.py index d65fc3d..cb3c9c2 100644 --- a/xradar/io/backends/nexrad_archive.py +++ b/xradar/io/backends/nexrad_archive.py @@ -1,8 +1,20 @@ import warnings import numpy as np - -from ...model import ( +import xarray as xr +from xarray.backends.common import BackendEntrypoint +from xarray.core import indexing +from xarray.core.variable import Variable + +from xradar import util +from xradar.io.backends.common import LazyLoadDict, prepare_for_read +from xradar.io.backends.nexrad_common import get_nexrad_location +from xradar.io.backends.nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) +from xradar.io.backends.nexrad_level2 import NEXRADLevel2File +from xradar.model import ( get_altitude_attrs, get_azimuth_attrs, get_elevation_attrs, @@ -12,10 +24,6 @@ get_range_attrs, get_time_attrs, ) -from .common import LazyLoadDict, prepare_for_read -from .nexrad_common import get_nexrad_location -from .nexrad_interpolate import _fast_interpolate_scan_2, _fast_interpolate_scan_4 -from .nexrad_level2 import NEXRADLevel2File nexrad_mapping = { "REF": "DBZH", @@ -106,7 +114,7 @@ def _interpolate_scan(mdata, start, end, moment_ngates, multiplier, linear_inter mdata[:] = np.ma.array(data, mask=(data == fill_value)) -def open_dataset( +def format_nexrad_data( filename, field_names=None, additional_metadata=None, @@ -118,10 +126,14 @@ def open_dataset( scans=None, linear_interp=True, storage_options={"anon": True}, + first_dimension=None, + group=None, **kwargs, ): # Load the data file in using NEXRADLevel2File Class - nfile = NEXRADLevel2File(prepare_for_read(filename)) + nfile = NEXRADLevel2File( + prepare_for_read(filename, storage_options=storage_options) + ) # Access the scan info and load in the available moments scan_info = nfile.scan_info(scans) @@ -179,9 +191,9 @@ def open_dataset( lat, lon, alt = get_nexrad_location(nfile.volume_header["icao"].decode()) else: lat, lon, alt = nfile.location() - latitude["data"] = np.array([lat], dtype="float64") - longitude["data"] = np.array([lon], dtype="float64") - altitude["data"] = np.array([alt], dtype="float64") + latitude["data"] = np.array(lat, dtype="float64") + longitude["data"] = np.array(lon, dtype="float64") + altitude["data"] = np.array(alt, dtype="float64") # Sweep information sweep_number = { @@ -284,4 +296,162 @@ def open_dataset( ) dic["data"] = mdata fields[nexrad_mapping[moment]] = dic - return fields + instrument_parameters = {} + + return create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters, + first_dimension=first_dimension, + group=group, + ) + + +def create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters, + first_dimension=None, + group=None, +): + """ + Creates an xarray dataset given a selection of fields defined as a collection of dictionaries + """ + if first_dimension is None: + first_dimension = "azimuth" + + if group is not None: + encoding = {"group": f"sweep_{group}"} + else: + encoding = {} + group = 0 + + group_slice = slice( + sweep_start_ray_index["data"][group], sweep_end_ray_index["data"][group] + ) + + # Track the number of sweeps + nsweeps = len(sweep_number["data"]) + + # Handle the coordinates + azimuth["data"] = indexing.LazilyOuterIndexedArray(azimuth["data"][group_slice]) + elevation["data"] = indexing.LazilyOuterIndexedArray(elevation["data"][group_slice]) + time["data"] = indexing.LazilyOuterIndexedArray(time["data"][group_slice]) + sweep_mode["data"] = sweep_mode["data"][group] + sweep_number["data"] = sweep_number["data"][group] + fixed_angle["data"] = fixed_angle["data"][group] + + coords = { + "azimuth": dictionary_to_xarray_variable(azimuth, encoding, first_dimension), + "elevation": dictionary_to_xarray_variable( + elevation, encoding, first_dimension + ), + "time": dictionary_to_xarray_variable(time, encoding, first_dimension), + "range": dictionary_to_xarray_variable(_range, encoding, "range"), + "longitude": dictionary_to_xarray_variable(longitude), + "latitude": dictionary_to_xarray_variable(latitude), + "altitude": dictionary_to_xarray_variable(altitude), + "sweep_mode": dictionary_to_xarray_variable(sweep_mode), + "sweep_number": dictionary_to_xarray_variable(sweep_number), + "sweep_fixed_angle": dictionary_to_xarray_variable(fixed_angle), + } + data_vars = { + field: dictionary_to_xarray_variable( + data, dims=(first_dimension, "range"), subset=group_slice + ) + for (field, data) in fields.items() + } + + ds = xr.Dataset(data_vars=data_vars, coords=coords, attrs=metadata) + + ds.attrs["nsweeps"] = nsweeps + + return ds + + +def dictionary_to_xarray_variable(variable_dict, encoding=None, dims=None, subset=None): + attrs = variable_dict.copy() + data = variable_dict.copy() + attrs.pop("data") + if dims is None: + dims = () + if subset is not None: + data["data"] = data["data"][subset] + return Variable(dims, data["data"], encoding, attrs) + + +class NexradBackendEntrypoint(BackendEntrypoint): + """ + Xarray BackendEntrypoint for NEXRAD Level2 Data + """ + + description = "Open NEXRAD Level2 files in Xarray" + url = "tbd" + + def open_dataset( + self, + filename_or_obj, + *, + mask_and_scale=True, + decode_times=True, + concat_characters=True, + decode_coords=True, + drop_variables=None, + use_cftime=None, + decode_timedelta=None, + group=None, + first_dim="auto", + reindex_angle=False, + fix_second_angle=False, + site_coords=True, + optional=True, + ): + ds = format_nexrad_data(filename_or_obj, group=group) + + # reassign azimuth/elevation/time coordinates + ds = ds.assign_coords({"azimuth": ds.azimuth}) + ds = ds.assign_coords({"elevation": ds.elevation}) + ds = ds.assign_coords({"time": ds.time}) + + ds.encoding["engine"] = "NexradLevel2" + + # handle duplicates and reindex + if decode_coords and reindex_angle is not False: + ds = ds.pipe(util.remove_duplicate_rays) + ds = ds.pipe(util.reindex_angle, **reindex_angle) + ds = ds.pipe(util.ipol_time, **reindex_angle) + + dim0 = "elevation" if ds.sweep_mode.load() == "rhi" else "azimuth" + + # todo: could be optimized + if first_dim == "time": + ds = ds.swap_dims({dim0: "time"}) + ds = ds.sortby("time") + else: + ds = ds.sortby(dim0) + + return ds From fbd942cb07b5661ab9eda41e85e7d4778c3d0cb8 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 13:57:56 -0600 Subject: [PATCH 05/48] FIX: Fix last few steps --- pyproject.toml | 1 + xradar/io/backends/__init__.py | 3 +- xradar/io/backends/nexrad_archive.py | 457 ------------------------ xradar/io/backends/nexrad_level2.py | 516 ++++++++++++++++++++++++++- 4 files changed, 510 insertions(+), 467 deletions(-) delete mode 100644 xradar/io/backends/nexrad_archive.py diff --git a/pyproject.toml b/pyproject.toml index 9d24919..9ac1069 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ gamic = "xradar.io.backends:GamicBackendEntrypoint" iris = "xradar.io.backends:IrisBackendEntrypoint" odim = "xradar.io.backends:OdimBackendEntrypoint" rainbow = "xradar.io.backends:RainbowBackendEntrypoint" +nexradlevel2 = "xradar.io.backends:NexradLevel2BackendEntrypoint" [build-system] requires = [ diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 509b11f..7eded4a 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -15,6 +15,7 @@ .. automodule:: xradar.io.backends.furuno .. automodule:: xradar.io.backends.rainbow .. automodule:: xradar.io.backends.iris +.. automodule:: xradar.io.backends.nexrad_level2 """ @@ -25,6 +26,6 @@ from .odim import * # noqa from .rainbow import * # noqa from .nexrad_level2 import * # noqa -from .nexrad_archive import * # noqa +from .nexrad_interpolate import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/nexrad_archive.py b/xradar/io/backends/nexrad_archive.py deleted file mode 100644 index cb3c9c2..0000000 --- a/xradar/io/backends/nexrad_archive.py +++ /dev/null @@ -1,457 +0,0 @@ -import warnings - -import numpy as np -import xarray as xr -from xarray.backends.common import BackendEntrypoint -from xarray.core import indexing -from xarray.core.variable import Variable - -from xradar import util -from xradar.io.backends.common import LazyLoadDict, prepare_for_read -from xradar.io.backends.nexrad_common import get_nexrad_location -from xradar.io.backends.nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) -from xradar.io.backends.nexrad_level2 import NEXRADLevel2File -from xradar.model import ( - get_altitude_attrs, - get_azimuth_attrs, - get_elevation_attrs, - get_latitude_attrs, - get_longitude_attrs, - get_moment_attrs, - get_range_attrs, - get_time_attrs, -) - -nexrad_mapping = { - "REF": "DBZH", - "VEL": "VRADH", - "SW": "WRADH", - "ZDR": "ZDR", - "PHI": "PHIDP", - "RHO": "RHOHV", - "CFP": "CCORH", -} - - -class _NEXRADLevel2StagedField: - """ - A class to facilitate on demand loading of field data from a Level 2 file. - """ - - def __init__(self, nfile, moment, max_ngates, scans): - """initialize.""" - self.nfile = nfile - self.moment = moment - self.max_ngates = max_ngates - self.scans = scans - - def __call__(self): - """Return the array containing the field data.""" - return self.nfile.get_data(self.moment, self.max_ngates, scans=self.scans) - - -def _find_range_params(scan_info): - """Return range parameters, first_gate, gate_spacing, last_gate.""" - min_first_gate = 999999 - min_gate_spacing = 999999 - max_last_gate = 0 - for scan_params in scan_info: - ngates = scan_params["ngates"][0] - for i, moment in enumerate(scan_params["moments"]): - first_gate = scan_params["first_gate"][i] - gate_spacing = scan_params["gate_spacing"][i] - last_gate = first_gate + gate_spacing * (ngates - 0.5) - - min_first_gate = min(min_first_gate, first_gate) - min_gate_spacing = min(min_gate_spacing, gate_spacing) - max_last_gate = max(max_last_gate, last_gate) - return min_first_gate, min_gate_spacing, max_last_gate - - -def _find_scans_to_interp(scan_info, first_gate, gate_spacing): - """Return a dict indicating what moments/scans need interpolation.""" - moments = {m for scan in scan_info for m in scan["moments"]} - interpolate = {moment: [] for moment in moments} - for scan_num, scan in enumerate(scan_info): - for moment in moments: - if moment not in scan["moments"]: - continue - index = scan["moments"].index(moment) - first = scan["first_gate"][index] - spacing = scan["gate_spacing"][index] - if first != first_gate or spacing != gate_spacing: - interpolate[moment].append(scan_num) - # for proper interpolation the gate spacing of the scan to be - # interpolated should be 1/4th the spacing of the radar - if spacing == gate_spacing * 4: - interpolate["multiplier"] = "4" - elif spacing == gate_spacing * 2: - interpolate["multiplier"] = "2" - else: - raise ValueError("Gate spacing is neither 1/4 or 1/2") - # assert first_gate + 1.5 * gate_spacing == first - # remove moments with no scans needing interpolation - interpolate = {k: v for k, v in interpolate.items() if len(v) != 0} - return interpolate - - -def _interpolate_scan(mdata, start, end, moment_ngates, multiplier, linear_interp=True): - """Interpolate a single NEXRAD moment scan from 1000 m to 250 m.""" - fill_value = -9999 - data = mdata.filled(fill_value) - scratch_ray = np.empty((data.shape[1],), dtype=data.dtype) - if multiplier == "4": - _fast_interpolate_scan_4( - data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp - ) - else: - _fast_interpolate_scan_2( - data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp - ) - mdata[:] = np.ma.array(data, mask=(data == fill_value)) - - -def format_nexrad_data( - filename, - field_names=None, - additional_metadata=None, - file_field_names=False, - exclude_fields=None, - include_fields=None, - delay_field_loading=False, - station=None, - scans=None, - linear_interp=True, - storage_options={"anon": True}, - first_dimension=None, - group=None, - **kwargs, -): - # Load the data file in using NEXRADLevel2File Class - nfile = NEXRADLevel2File( - prepare_for_read(filename, storage_options=storage_options) - ) - - # Access the scan info and load in the available moments - scan_info = nfile.scan_info(scans) - available_moments = {m for scan in scan_info for m in scan["moments"]} - first_gate, gate_spacing, last_gate = _find_range_params(scan_info) - - # Interpolate to 360 degrees where neccessary - interpolate = _find_scans_to_interp(scan_info, first_gate, gate_spacing) - - # Deal with time - time_start, _time = nfile.get_times(scans) - time = get_time_attrs(time_start) - time["data"] = _time - - # range - _range = get_range_attrs() - first_gate, gate_spacing, last_gate = _find_range_params(scan_info) - _range["data"] = np.arange(first_gate, last_gate, gate_spacing, "float32") - _range["meters_to_center_of_first_gate"] = float(first_gate) - _range["meters_between_gates"] = float(gate_spacing) - - # metadata - metadata = { - "Conventions": "CF/Radial instrument_parameters", - "version": "1.3", - "title": "", - "institution": "", - "references": "", - "source": "", - "history": "", - "comment": "", - "intrument_name": "", - "original_container": "NEXRAD Level II", - } - - vcp_pattern = nfile.get_vcp_pattern() - if vcp_pattern is not None: - metadata["vcp_pattern"] = vcp_pattern - if "icao" in nfile.volume_header.keys(): - metadata["instrument_name"] = nfile.volume_header["icao"].decode() - - # scan_type - - # latitude, longitude, altitude - latitude = get_latitude_attrs() - longitude = get_longitude_attrs() - altitude = get_altitude_attrs() - - if nfile._msg_type == "1" and station is not None: - lat, lon, alt = get_nexrad_location(station) - elif ( - "icao" in nfile.volume_header.keys() - and nfile.volume_header["icao"].decode()[0] == "T" - ): - lat, lon, alt = get_nexrad_location(nfile.volume_header["icao"].decode()) - else: - lat, lon, alt = nfile.location() - latitude["data"] = np.array(lat, dtype="float64") - longitude["data"] = np.array(lon, dtype="float64") - altitude["data"] = np.array(alt, dtype="float64") - - # Sweep information - sweep_number = { - "units": "count", - "standard_name": "sweep_number", - "long_name": "Sweep number", - } - - sweep_mode = { - "units": "unitless", - "standard_name": "sweep_mode", - "long_name": "Sweep mode", - "comment": 'Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi"', - } - sweep_start_ray_index = { - "long_name": "Index of first ray in sweep, 0-based", - "units": "count", - } - sweep_end_ray_index = { - "long_name": "Index of last ray in sweep, 0-based", - "units": "count", - } - - if scans is None: - nsweeps = int(nfile.nscans) - else: - nsweeps = len(scans) - sweep_number["data"] = np.arange(nsweeps, dtype="int32") - sweep_mode["data"] = np.array(nsweeps * ["azimuth_surveillance"], dtype="S") - - rays_per_scan = [s["nrays"] for s in scan_info] - sweep_end_ray_index["data"] = np.cumsum(rays_per_scan, dtype="int32") - 1 - - rays_per_scan.insert(0, 0) - sweep_start_ray_index["data"] = np.cumsum(rays_per_scan[:-1], dtype="int32") - - # azimuth, elevation, fixed_angle - azimuth = get_azimuth_attrs() - elevation = get_elevation_attrs() - fixed_angle = { - "long_name": "Target angle for sweep", - "units": "degrees", - "standard_name": "target_fixed_angle", - } - - azimuth["data"] = nfile.get_azimuth_angles(scans) - elevation["data"] = nfile.get_elevation_angles(scans).astype("float32") - fixed_agl = [] - for i in nfile.get_target_angles(scans): - if i > 180: - i = i - 360.0 - warnings.warn( - "Fixed_angle(s) greater than 180 degrees present." - " Assuming angle to be negative so subtrating 360", - UserWarning, - ) - else: - i = i - fixed_agl.append(i) - fixed_angles = np.array(fixed_agl, dtype="float32") - fixed_angle["data"] = fixed_angles - - # fields - max_ngates = len(_range["data"]) - available_moments = {m for scan in scan_info for m in scan["moments"]} - interpolate = _find_scans_to_interp( - scan_info, - first_gate, - gate_spacing, - ) - - fields = {} - for moment in available_moments: - dic = get_moment_attrs(nexrad_mapping[moment]) - dic["_FillValue"] = -9999 - if delay_field_loading and moment not in interpolate: - dic = LazyLoadDict(dic) - data_call = _NEXRADLevel2StagedField(nfile, moment, max_ngates, scans) - dic.set_lazy("data", data_call) - else: - mdata = nfile.get_data(moment, max_ngates, scans=scans) - if moment in interpolate: - interp_scans = interpolate[moment] - warnings.warn( - "Gate spacing is not constant, interpolating data in " - + f"scans {interp_scans} for moment {moment}.", - UserWarning, - ) - for scan in interp_scans: - idx = scan_info[scan]["moments"].index(moment) - moment_ngates = scan_info[scan]["ngates"][idx] - start = sweep_start_ray_index["data"][scan] - end = sweep_end_ray_index["data"][scan] - if interpolate["multiplier"] == "4": - multiplier = "4" - else: - multiplier = "2" - _interpolate_scan( - mdata, start, end, moment_ngates, multiplier, linear_interp - ) - dic["data"] = mdata - fields[nexrad_mapping[moment]] = dic - instrument_parameters = {} - - return create_dataset_from_fields( - time, - _range, - fields, - metadata, - latitude, - longitude, - altitude, - sweep_number, - sweep_mode, - fixed_angle, - sweep_start_ray_index, - sweep_end_ray_index, - azimuth, - elevation, - instrument_parameters, - first_dimension=first_dimension, - group=group, - ) - - -def create_dataset_from_fields( - time, - _range, - fields, - metadata, - latitude, - longitude, - altitude, - sweep_number, - sweep_mode, - fixed_angle, - sweep_start_ray_index, - sweep_end_ray_index, - azimuth, - elevation, - instrument_parameters, - first_dimension=None, - group=None, -): - """ - Creates an xarray dataset given a selection of fields defined as a collection of dictionaries - """ - if first_dimension is None: - first_dimension = "azimuth" - - if group is not None: - encoding = {"group": f"sweep_{group}"} - else: - encoding = {} - group = 0 - - group_slice = slice( - sweep_start_ray_index["data"][group], sweep_end_ray_index["data"][group] - ) - - # Track the number of sweeps - nsweeps = len(sweep_number["data"]) - - # Handle the coordinates - azimuth["data"] = indexing.LazilyOuterIndexedArray(azimuth["data"][group_slice]) - elevation["data"] = indexing.LazilyOuterIndexedArray(elevation["data"][group_slice]) - time["data"] = indexing.LazilyOuterIndexedArray(time["data"][group_slice]) - sweep_mode["data"] = sweep_mode["data"][group] - sweep_number["data"] = sweep_number["data"][group] - fixed_angle["data"] = fixed_angle["data"][group] - - coords = { - "azimuth": dictionary_to_xarray_variable(azimuth, encoding, first_dimension), - "elevation": dictionary_to_xarray_variable( - elevation, encoding, first_dimension - ), - "time": dictionary_to_xarray_variable(time, encoding, first_dimension), - "range": dictionary_to_xarray_variable(_range, encoding, "range"), - "longitude": dictionary_to_xarray_variable(longitude), - "latitude": dictionary_to_xarray_variable(latitude), - "altitude": dictionary_to_xarray_variable(altitude), - "sweep_mode": dictionary_to_xarray_variable(sweep_mode), - "sweep_number": dictionary_to_xarray_variable(sweep_number), - "sweep_fixed_angle": dictionary_to_xarray_variable(fixed_angle), - } - data_vars = { - field: dictionary_to_xarray_variable( - data, dims=(first_dimension, "range"), subset=group_slice - ) - for (field, data) in fields.items() - } - - ds = xr.Dataset(data_vars=data_vars, coords=coords, attrs=metadata) - - ds.attrs["nsweeps"] = nsweeps - - return ds - - -def dictionary_to_xarray_variable(variable_dict, encoding=None, dims=None, subset=None): - attrs = variable_dict.copy() - data = variable_dict.copy() - attrs.pop("data") - if dims is None: - dims = () - if subset is not None: - data["data"] = data["data"][subset] - return Variable(dims, data["data"], encoding, attrs) - - -class NexradBackendEntrypoint(BackendEntrypoint): - """ - Xarray BackendEntrypoint for NEXRAD Level2 Data - """ - - description = "Open NEXRAD Level2 files in Xarray" - url = "tbd" - - def open_dataset( - self, - filename_or_obj, - *, - mask_and_scale=True, - decode_times=True, - concat_characters=True, - decode_coords=True, - drop_variables=None, - use_cftime=None, - decode_timedelta=None, - group=None, - first_dim="auto", - reindex_angle=False, - fix_second_angle=False, - site_coords=True, - optional=True, - ): - ds = format_nexrad_data(filename_or_obj, group=group) - - # reassign azimuth/elevation/time coordinates - ds = ds.assign_coords({"azimuth": ds.azimuth}) - ds = ds.assign_coords({"elevation": ds.elevation}) - ds = ds.assign_coords({"time": ds.time}) - - ds.encoding["engine"] = "NexradLevel2" - - # handle duplicates and reindex - if decode_coords and reindex_angle is not False: - ds = ds.pipe(util.remove_duplicate_rays) - ds = ds.pipe(util.reindex_angle, **reindex_angle) - ds = ds.pipe(util.ipol_time, **reindex_angle) - - dim0 = "elevation" if ds.sweep_mode.load() == "rhi" else "azimuth" - - # todo: could be optimized - if first_dim == "time": - ds = ds.swap_dims({dim0: "time"}) - ds = ds.sortby("time") - else: - ds = ds.sortby(dim0) - - return ds diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index e44bc14..5f32fc8 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -1,18 +1,47 @@ -""" -Functions for reading NEXRAD level 2 files. -""" - - -# This file has been adapted from is Py-ART, the Python ARM Radar Toolkit -# https://github.com/ARM-DOE/pyart - - import bz2 import struct import warnings from datetime import datetime, timedelta import numpy as np +import xarray as xr +from datatree import DataTree +from xarray.backends.common import BackendEntrypoint +from xarray.core import indexing +from xarray.core.variable import Variable + +from xradar import util +from xradar.io.backends.common import ( + LazyLoadDict, + _assign_root, + _attach_sweep_groups, + prepare_for_read, +) +from xradar.io.backends.nexrad_common import get_nexrad_location +from xradar.io.backends.nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) +from xradar.model import ( + get_altitude_attrs, + get_azimuth_attrs, + get_elevation_attrs, + get_latitude_attrs, + get_longitude_attrs, + get_moment_attrs, + get_range_attrs, + get_time_attrs, +) + +nexrad_mapping = { + "REF": "DBZH", + "VEL": "VRADH", + "SW": "WRADH", + "ZDR": "ZDR", + "PHI": "PHIDP", + "RHO": "RHOHV", + "CFP": "CCORH", +} class _NEXRADLevel2StagedField: @@ -958,3 +987,472 @@ def _unpack_structure(string, structure): ("nyquist_vel", SINT2), ("spare", "2s"), ) + + +def _find_range_params(scan_info): + """Return range parameters, first_gate, gate_spacing, last_gate.""" + min_first_gate = 999999 + min_gate_spacing = 999999 + max_last_gate = 0 + for scan_params in scan_info: + ngates = scan_params["ngates"][0] + for i, moment in enumerate(scan_params["moments"]): + first_gate = scan_params["first_gate"][i] + gate_spacing = scan_params["gate_spacing"][i] + last_gate = first_gate + gate_spacing * (ngates - 0.5) + + min_first_gate = min(min_first_gate, first_gate) + min_gate_spacing = min(min_gate_spacing, gate_spacing) + max_last_gate = max(max_last_gate, last_gate) + return min_first_gate, min_gate_spacing, max_last_gate + + +def _find_scans_to_interp(scan_info, first_gate, gate_spacing): + """Return a dict indicating what moments/scans need interpolation.""" + moments = {m for scan in scan_info for m in scan["moments"]} + interpolate = {moment: [] for moment in moments} + for scan_num, scan in enumerate(scan_info): + for moment in moments: + if moment not in scan["moments"]: + continue + index = scan["moments"].index(moment) + first = scan["first_gate"][index] + spacing = scan["gate_spacing"][index] + if first != first_gate or spacing != gate_spacing: + interpolate[moment].append(scan_num) + # for proper interpolation the gate spacing of the scan to be + # interpolated should be 1/4th the spacing of the radar + if spacing == gate_spacing * 4: + interpolate["multiplier"] = "4" + elif spacing == gate_spacing * 2: + interpolate["multiplier"] = "2" + else: + raise ValueError("Gate spacing is neither 1/4 or 1/2") + # assert first_gate + 1.5 * gate_spacing == first + # remove moments with no scans needing interpolation + interpolate = {k: v for k, v in interpolate.items() if len(v) != 0} + return interpolate + + +def _interpolate_scan(mdata, start, end, moment_ngates, multiplier, linear_interp=True): + """Interpolate a single NEXRAD moment scan from 1000 m to 250 m.""" + fill_value = -9999 + data = mdata.filled(fill_value) + scratch_ray = np.empty((data.shape[1],), dtype=data.dtype) + if multiplier == "4": + _fast_interpolate_scan_4( + data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp + ) + else: + _fast_interpolate_scan_2( + data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp + ) + mdata[:] = np.ma.array(data, mask=(data == fill_value)) + + +def format_nexrad_level2_data( + filename, + delay_field_loading=False, + station=None, + scans=None, + linear_interp=True, + storage_options={"anon": True}, + first_dimension=None, + group=None, + **kwargs, +): + # Load the data file in using NEXRADLevel2File Class + nfile = NEXRADLevel2File( + prepare_for_read(filename, storage_options=storage_options) + ) + + # Access the scan info and load in the available moments + scan_info = nfile.scan_info(scans) + available_moments = {m for scan in scan_info for m in scan["moments"]} + first_gate, gate_spacing, last_gate = _find_range_params(scan_info) + + # Interpolate to 360 degrees where neccessary + interpolate = _find_scans_to_interp(scan_info, first_gate, gate_spacing) + + # Deal with time + time_start, _time = nfile.get_times(scans) + time = get_time_attrs(time_start) + time["data"] = _time + + # range + _range = get_range_attrs() + first_gate, gate_spacing, last_gate = _find_range_params(scan_info) + _range["data"] = np.arange(first_gate, last_gate, gate_spacing, "float32") + _range["meters_to_center_of_first_gate"] = float(first_gate) + _range["meters_between_gates"] = float(gate_spacing) + + # metadata + metadata = { + "Conventions": "CF/Radial instrument_parameters", + "version": "1.3", + "title": "", + "institution": "", + "references": "", + "source": "", + "history": "", + "comment": "", + "intrument_name": "", + "original_container": "NEXRAD Level II", + } + + vcp_pattern = nfile.get_vcp_pattern() + if vcp_pattern is not None: + metadata["vcp_pattern"] = vcp_pattern + if "icao" in nfile.volume_header.keys(): + metadata["instrument_name"] = nfile.volume_header["icao"].decode() + + # scan_type + + # latitude, longitude, altitude + latitude = get_latitude_attrs() + longitude = get_longitude_attrs() + altitude = get_altitude_attrs() + + if nfile._msg_type == "1" and station is not None: + lat, lon, alt = get_nexrad_location(station) + elif ( + "icao" in nfile.volume_header.keys() + and nfile.volume_header["icao"].decode()[0] == "T" + ): + lat, lon, alt = get_nexrad_location(nfile.volume_header["icao"].decode()) + else: + lat, lon, alt = nfile.location() + latitude["data"] = np.array(lat, dtype="float64") + longitude["data"] = np.array(lon, dtype="float64") + altitude["data"] = np.array(alt, dtype="float64") + + # Sweep information + sweep_number = { + "units": "count", + "standard_name": "sweep_number", + "long_name": "Sweep number", + } + + sweep_mode = { + "units": "unitless", + "standard_name": "sweep_mode", + "long_name": "Sweep mode", + "comment": 'Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi"', + } + sweep_start_ray_index = { + "long_name": "Index of first ray in sweep, 0-based", + "units": "count", + } + sweep_end_ray_index = { + "long_name": "Index of last ray in sweep, 0-based", + "units": "count", + } + + if scans is None: + nsweeps = int(nfile.nscans) + else: + nsweeps = len(scans) + sweep_number["data"] = np.arange(nsweeps, dtype="int32") + sweep_mode["data"] = np.array(nsweeps * ["azimuth_surveillance"], dtype="S") + + rays_per_scan = [s["nrays"] for s in scan_info] + sweep_end_ray_index["data"] = np.cumsum(rays_per_scan, dtype="int32") - 1 + + rays_per_scan.insert(0, 0) + sweep_start_ray_index["data"] = np.cumsum(rays_per_scan[:-1], dtype="int32") + + # azimuth, elevation, fixed_angle + azimuth = get_azimuth_attrs() + elevation = get_elevation_attrs() + fixed_angle = { + "long_name": "Target angle for sweep", + "units": "degrees", + "standard_name": "target_fixed_angle", + } + + azimuth["data"] = nfile.get_azimuth_angles(scans) + elevation["data"] = nfile.get_elevation_angles(scans).astype("float32") + fixed_agl = [] + for i in nfile.get_target_angles(scans): + if i > 180: + i = i - 360.0 + warnings.warn( + "Fixed_angle(s) greater than 180 degrees present." + " Assuming angle to be negative so subtrating 360", + UserWarning, + ) + else: + i = i + fixed_agl.append(i) + fixed_angles = np.array(fixed_agl, dtype="float32") + fixed_angle["data"] = fixed_angles + + # fields + max_ngates = len(_range["data"]) + available_moments = {m for scan in scan_info for m in scan["moments"]} + interpolate = _find_scans_to_interp( + scan_info, + first_gate, + gate_spacing, + ) + + fields = {} + for moment in available_moments: + dic = get_moment_attrs(nexrad_mapping[moment]) + dic["_FillValue"] = -9999 + if delay_field_loading and moment not in interpolate: + dic = LazyLoadDict(dic) + data_call = _NEXRADLevel2StagedField(nfile, moment, max_ngates, scans) + dic.set_lazy("data", data_call) + else: + mdata = nfile.get_data(moment, max_ngates, scans=scans) + if moment in interpolate: + interp_scans = interpolate[moment] + warnings.warn( + "Gate spacing is not constant, interpolating data in " + + f"scans {interp_scans} for moment {moment}.", + UserWarning, + ) + for scan in interp_scans: + idx = scan_info[scan]["moments"].index(moment) + moment_ngates = scan_info[scan]["ngates"][idx] + start = sweep_start_ray_index["data"][scan] + end = sweep_end_ray_index["data"][scan] + if interpolate["multiplier"] == "4": + multiplier = "4" + else: + multiplier = "2" + _interpolate_scan( + mdata, start, end, moment_ngates, multiplier, linear_interp + ) + dic["data"] = mdata + fields[nexrad_mapping[moment]] = dic + instrument_parameters = {} + + return create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters, + first_dimension=first_dimension, + group=group, + ) + + +def create_dataset_from_fields( + time, + _range, + fields, + metadata, + latitude, + longitude, + altitude, + sweep_number, + sweep_mode, + fixed_angle, + sweep_start_ray_index, + sweep_end_ray_index, + azimuth, + elevation, + instrument_parameters={}, + first_dimension=None, + group=None, +): + """ + Creates an xarray dataset given a selection of fields defined as a collection of dictionaries + """ + if first_dimension is None: + first_dimension = "azimuth" + + if group is not None: + encoding = {"group": group} + group = int(group[6:]) + else: + encoding = {} + group = 0 + + group_slice = slice( + sweep_start_ray_index["data"][group], sweep_end_ray_index["data"][group] + ) + + # Track the number of sweeps + nsweeps = len(sweep_number["data"]) + + # Handle the coordinates + azimuth["data"] = indexing.LazilyOuterIndexedArray(azimuth["data"][group_slice]) + elevation["data"] = indexing.LazilyOuterIndexedArray(elevation["data"][group_slice]) + time["data"] = indexing.LazilyOuterIndexedArray(time["data"][group_slice]) + sweep_mode["data"] = sweep_mode["data"][group] + sweep_number["data"] = sweep_number["data"][group] + fixed_angle["data"] = fixed_angle["data"][group] + + coords = { + "azimuth": dictionary_to_xarray_variable(azimuth, encoding, first_dimension), + "elevation": dictionary_to_xarray_variable( + elevation, encoding, first_dimension + ), + "time": dictionary_to_xarray_variable(time, encoding, first_dimension), + "range": dictionary_to_xarray_variable(_range, encoding, "range"), + "longitude": dictionary_to_xarray_variable(longitude), + "latitude": dictionary_to_xarray_variable(latitude), + "altitude": dictionary_to_xarray_variable(altitude), + "sweep_mode": dictionary_to_xarray_variable(sweep_mode), + "sweep_number": dictionary_to_xarray_variable(sweep_number), + "sweep_fixed_angle": dictionary_to_xarray_variable(fixed_angle), + } + data_vars = { + field: dictionary_to_xarray_variable( + data, dims=(first_dimension, "range"), subset=group_slice + ) + for (field, data) in fields.items() + } + + ds = xr.Dataset(data_vars=data_vars, coords=coords, attrs=metadata) + + ds.attrs["nsweeps"] = nsweeps + + return ds + + +def dictionary_to_xarray_variable(variable_dict, encoding=None, dims=None, subset=None): + attrs = variable_dict.copy() + data = variable_dict.copy() + attrs.pop("data") + if dims is None: + dims = () + if subset is not None: + data["data"] = data["data"][subset] + return Variable(dims, data["data"], encoding, attrs) + + +class NexradLevel2BackendEntrypoint(BackendEntrypoint): + """ + Xarray BackendEntrypoint for NEXRAD Level2 Data + """ + + description = "Open NEXRAD Level2 files in Xarray" + url = "tbd" + + def open_dataset( + self, + filename_or_obj, + *, + decode_coords=True, + mask_and_scale=True, + decode_times=True, + concat_characters=True, + drop_variables=None, + use_cftime=None, + decode_timedelta=None, + group=None, + first_dim="auto", + reindex_angle=False, + ): + ds = format_nexrad_level2_data(filename_or_obj, group=group) + + # reassign azimuth/elevation/time coordinates + ds = ds.assign_coords({"azimuth": ds.azimuth}) + ds = ds.assign_coords({"elevation": ds.elevation}) + ds = ds.assign_coords({"time": ds.time}) + + ds.encoding["engine"] = "NexradLevel2" + + # handle duplicates and reindex + if decode_coords and reindex_angle is not False: + ds = ds.pipe(util.remove_duplicate_rays) + ds = ds.pipe(util.reindex_angle, **reindex_angle) + ds = ds.pipe(util.ipol_time, **reindex_angle) + + dim0 = "elevation" if ds.sweep_mode.load() == "rhi" else "azimuth" + + # todo: could be optimized + if first_dim == "time": + ds = ds.swap_dims({dim0: "time"}) + ds = ds.sortby("time") + else: + ds = ds.sortby(dim0) + + return ds + + +def _get_nexradlevel2_number_of_sweeps(filename_or_obj, **kwargs): + number_of_sweeps = xr.open_dataset( + filename_or_obj, engine="nexradlevel2", **kwargs + ).nsweeps + return [f"sweep_{sweep_number}" for sweep_number in np.arange(0, number_of_sweeps)] + + +def open_nexradlevel2_datatree(filename_or_obj, **kwargs): + """Open NEXRAD Level2 dataset as :py:class:`datatree.DataTree`. + + Parameters + ---------- + filename_or_obj : str, Path, file-like or DataStore + Strings and Path objects are interpreted as a path to a local or remote + radar file + + Keyword Arguments + ----------------- + sweep : int, list of int, optional + Sweep number(s) to extract, default to first sweep. If None, all sweeps are + extracted into a list. + first_dim : str + Can be ``time`` or ``auto`` first dimension. If set to ``auto``, + first dimension will be either ``azimuth`` or ``elevation`` depending on + type of sweep. Defaults to ``auto``. + reindex_angle : bool or dict + Defaults to False, no reindexing. Given dict should contain the kwargs to + reindex_angle. Only invoked if `decode_coord=True`. + fix_second_angle : bool + If True, fixes erroneous second angle data. Defaults to ``False``. + site_coords : bool + Attach radar site-coordinates to Dataset, defaults to ``True``. + kwargs : dict + Additional kwargs are fed to :py:func:`xarray.open_dataset`. + + Returns + ------- + dtree: datatree.DataTree + DataTree + """ + # handle kwargs, extract first_dim + backend_kwargs = kwargs.pop("backend_kwargs", {}) + # first_dim = backend_kwargs.pop("first_dim", None) + sweep = kwargs.pop("sweep", None) + sweeps = [] + kwargs["backend_kwargs"] = backend_kwargs + + if isinstance(sweep, str): + sweeps = [sweep] + elif isinstance(sweep, int): + sweeps = [f"sweep_{sweep}"] + elif isinstance(sweep, list): + if isinstance(sweep[0], int): + sweeps = [f"sweep_{i}" for i in sweep] + else: + sweeps.extend(sweep) + else: + sweeps = _get_nexradlevel2_number_of_sweeps(filename_or_obj) + + ds = [ + xr.open_dataset(filename_or_obj, group=swp, engine="nexradlevel2", **kwargs) + for swp in sweeps + ] + + ds.insert(0, xr.Dataset()) + + # create datatree root node with required data + dtree = DataTree(data=_assign_root(ds), name="root") + # return datatree with attached sweep child nodes + return _attach_sweep_groups(dtree, ds[1:]) From 3f8f8cfdd74e506486663aef42c3f9dbe7f52972 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 14:54:45 -0600 Subject: [PATCH 06/48] ADD: Add testing suite --- tests/conftest.py | 5 +++++ tests/io/test_nexrad_level2.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) create mode 100644 tests/io/test_nexrad_level2.py diff --git a/tests/conftest.py b/tests/conftest.py index b1a0ca1..e17a954 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -58,3 +58,8 @@ def iris0_file(): @pytest.fixture(scope="session") def iris1_file(): return DATASETS.fetch("SUR210819000227.RAWKPJV") + + +@pytest.fixture(scope="session") +def nexradlevel2_file(): + return DATASETS.fetch("KATX20130717_195021_V06") diff --git a/tests/io/test_nexrad_level2.py b/tests/io/test_nexrad_level2.py new file mode 100644 index 0000000..2a476e5 --- /dev/null +++ b/tests/io/test_nexrad_level2.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python +# Copyright (c) 2024, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +"""Tests for `xradar.io.nexrad_archive` module.""" + +import xarray as xr + +from xradar.io.backends import open_nexradlevel2_datatree + + +def test_open_nexradlevel2_datatree(nexradlevel2_file): + dtree = open_nexradlevel2_datatree(nexradlevel2_file) + ds = dtree["sweep_0"] + assert ds.attrs["instrument_name"] == "KATX" + assert ds.attrs["nsweeps"] == 16 + assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters" + assert ds["DBZH"].shape == (719, 1832) + assert ds["DBZH"].dims == ("azimuth", "range") + assert int(ds.sweep_number.values) == 0 + + +def test_open_nexrad_level2_backend(nexradlevel2_file): + ds = xr.open_dataset(nexradlevel2_file, engine="nexradlevel2") + assert ds.attrs["instrument_name"] == "KATX" + assert ds.attrs["nsweeps"] == 16 + assert ds.attrs["Conventions"] == "CF/Radial instrument_parameters" + assert ds["DBZH"].shape == (719, 1832) + assert ds["DBZH"].dims == ("azimuth", "range") + assert int(ds.sweep_number.values) == 0 From e042f1fa6a4df62aa5dab192967cf6a3b6736124 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:13:39 -0600 Subject: [PATCH 07/48] FIX: Fix the manifest --- MANIFEST.in | 2 +- xradar/io/backends/__init__.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/MANIFEST.in b/MANIFEST.in index 109f8ea..b533f51 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -6,6 +6,6 @@ include README.md recursive-include tests * recursive-exclude * __pycache__ -recursive-exclude * *.py[co] +recursive-exclude * *.py[o] recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 7eded4a..8c1b2e1 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -26,6 +26,5 @@ from .odim import * # noqa from .rainbow import * # noqa from .nexrad_level2 import * # noqa -from .nexrad_interpolate import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] From 0c2d0dcfdbb00d533dd394c0f3b06325911a3d06 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:21:54 -0600 Subject: [PATCH 08/48] FIX: Fix local import nexrad_interpolate --- xradar/io/backends/nexrad_level2.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 5f32fc8..7a79d8d 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -18,10 +18,6 @@ prepare_for_read, ) from xradar.io.backends.nexrad_common import get_nexrad_location -from xradar.io.backends.nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) from xradar.model import ( get_altitude_attrs, get_azimuth_attrs, @@ -33,6 +29,11 @@ get_time_attrs, ) +from .nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) + nexrad_mapping = { "REF": "DBZH", "VEL": "VRADH", @@ -1366,7 +1367,7 @@ def open_dataset( ds = ds.assign_coords({"elevation": ds.elevation}) ds = ds.assign_coords({"time": ds.time}) - ds.encoding["engine"] = "NexradLevel2" + ds.encoding["engine"] = "nexradlevel2" # handle duplicates and reindex if decode_coords and reindex_angle is not False: From 4f1d06c013dae8a6b7a1f21ed324f53dfeb622fb Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:26:31 -0600 Subject: [PATCH 09/48] FIX: Fix import of interpolation --- xradar/io/backends/nexrad_level2.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 7a79d8d..98806f9 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -18,6 +18,10 @@ prepare_for_read, ) from xradar.io.backends.nexrad_common import get_nexrad_location +from xradar.io.backends.nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) from xradar.model import ( get_altitude_attrs, get_azimuth_attrs, @@ -29,11 +33,6 @@ get_time_attrs, ) -from .nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) - nexrad_mapping = { "REF": "DBZH", "VEL": "VRADH", From 1af4649077167fde13b4dde06ea3bba34b56c43e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:33:33 -0600 Subject: [PATCH 10/48] ADD: Add missing cython depend --- ci/notebooktests.yml | 1 + ci/unittests.yml | 1 + 2 files changed, 2 insertions(+) diff --git a/ci/notebooktests.yml b/ci/notebooktests.yml index a0587a0..e0628db 100644 --- a/ci/notebooktests.yml +++ b/ci/notebooktests.yml @@ -7,6 +7,7 @@ dependencies: - cmweather - codecov - coverage + - cython - dask - h5netcdf - h5py diff --git a/ci/unittests.yml b/ci/unittests.yml index 9291f8e..d345200 100644 --- a/ci/unittests.yml +++ b/ci/unittests.yml @@ -6,6 +6,7 @@ dependencies: - codecov - cmweather - coverage + - cython - dask - h5netcdf - h5py From 2e3713026eaac69194718626ef2e6d979e6e167f Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:46:59 -0600 Subject: [PATCH 11/48] ADD: Add updated manifest --- MANIFEST.in | 1 + 1 file changed, 1 insertion(+) diff --git a/MANIFEST.in b/MANIFEST.in index b533f51..15c11fa 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -5,6 +5,7 @@ include LICENSE include README.md recursive-include tests * +recursive-include xradar *.pyx recursive-exclude * __pycache__ recursive-exclude * *.py[o] From 7ae79bdd32c93791be16623fce0a280d0244e2ad Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:49:10 -0600 Subject: [PATCH 12/48] ADD: Ensure cython is packaged --- MANIFEST.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MANIFEST.in b/MANIFEST.in index 15c11fa..3423154 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -5,7 +5,7 @@ include LICENSE include README.md recursive-include tests * -recursive-include xradar *.pyx +recursive-include xradar *.c *.pyx *.pxd recursive-exclude * __pycache__ recursive-exclude * *.py[o] From ec2261c3aa2cdd5fc87d6e049657870d2a142d87 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:52:31 -0600 Subject: [PATCH 13/48] ADD: Add specific submodules --- xradar/io/backends/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 8c1b2e1..81b0278 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -25,6 +25,9 @@ from .iris import * # noqa from .odim import * # noqa from .rainbow import * # noqa -from .nexrad_level2 import * # noqa +from .nexrad_level2 import ( + NexradLevel2BackendEntrypoint, # noqa + open_nexradlevel2_datatree, # noqa +) __all__ = [s for s in dir() if not s.startswith("_")] From 568e65b1d39c82cbc501acea980034a7d9a89d97 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 15:58:21 -0600 Subject: [PATCH 14/48] force reinstall --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index daa2f2b..8e09852 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,7 +58,7 @@ jobs: conda - name: Install xradar run: | - python -m pip install . --no-deps + python -m pip install . --no-deps --force-reinstall - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" From 0a56163e7ce001f492381c4a15e7695f5df1e3d4 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 19 Jan 2024 16:02:59 -0600 Subject: [PATCH 15/48] make sure more files are included --- MANIFEST.in | 1 - 1 file changed, 1 deletion(-) diff --git a/MANIFEST.in b/MANIFEST.in index 3423154..89afb7b 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -7,6 +7,5 @@ include README.md recursive-include tests * recursive-include xradar *.c *.pyx *.pxd recursive-exclude * __pycache__ -recursive-exclude * *.py[o] recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif From a904c0312b32152e94db71b9aa04bca656c77a65 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 14:28:48 -0600 Subject: [PATCH 16/48] FIX: Fix build of cython extensions --- setup.py | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.py b/setup.py index db3bcee..5a48f96 100644 --- a/setup.py +++ b/setup.py @@ -19,7 +19,6 @@ # This is the function that is executed setup( - name="xradar", # Required # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From 14aa01cd2989e462c2c727ec93f71f45db0f6b7b Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 14:31:27 -0600 Subject: [PATCH 17/48] FIX: Fix lowercase letter --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 9ac1069..6deab73 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,7 @@ requires = [ "setuptools>=45", "wheel", "setuptools_scm[toml]>=7.0", - "cython", + "Cython", "numpy" ] build-backend = "setuptools.build_meta" From 888fb1b35d1c8764bfb4bb473db2f3e248eaaf6b Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 14:44:26 -0600 Subject: [PATCH 18/48] revert couple of settings --- pyproject.toml | 2 +- setup.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 6deab73..9ac1069 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,7 @@ requires = [ "setuptools>=45", "wheel", "setuptools_scm[toml]>=7.0", - "Cython", + "cython", "numpy" ] build-backend = "setuptools.build_meta" diff --git a/setup.py b/setup.py index 5a48f96..3c60099 100644 --- a/setup.py +++ b/setup.py @@ -19,6 +19,7 @@ # This is the function that is executed setup( + name="xradar", # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From 44647acfba5317fde44f894eba8dccdb8865072e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 14:47:25 -0600 Subject: [PATCH 19/48] fix installation line --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8e09852..fb44926 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,7 +58,7 @@ jobs: conda - name: Install xradar run: | - python -m pip install . --no-deps --force-reinstall + python -m pip install . --force-reinstall - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" From 2b8d9b988a6cf6668f37edc98584c212b2737976 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 14:49:05 -0600 Subject: [PATCH 20/48] ADD: Add proper import --- xradar/io/backends/nexrad_level2.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 98806f9..7a79d8d 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -18,10 +18,6 @@ prepare_for_read, ) from xradar.io.backends.nexrad_common import get_nexrad_location -from xradar.io.backends.nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) from xradar.model import ( get_altitude_attrs, get_azimuth_attrs, @@ -33,6 +29,11 @@ get_time_attrs, ) +from .nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) + nexrad_mapping = { "REF": "DBZH", "VEL": "VRADH", From 0550e2ca45d7d6233b1fbcda018c85bdaf972949 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:00:57 -0600 Subject: [PATCH 21/48] move interpolation to its own submodule --- xradar/interpolate/__init__.py | 11 +++++++++++ .../backends => interpolate}/nexrad_interpolate.pyx | 0 xradar/io/backends/nexrad_level2.py | 9 ++++----- 3 files changed, 15 insertions(+), 5 deletions(-) create mode 100644 xradar/interpolate/__init__.py rename xradar/{io/backends => interpolate}/nexrad_interpolate.pyx (100%) diff --git a/xradar/interpolate/__init__.py b/xradar/interpolate/__init__.py new file mode 100644 index 0000000..6e1a1c9 --- /dev/null +++ b/xradar/interpolate/__init__.py @@ -0,0 +1,11 @@ +#!/usr/bin/env python +# Copyright (c) 2024, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +""" +XRadar Interpolation +==================== + +""" + +__all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/nexrad_interpolate.pyx b/xradar/interpolate/nexrad_interpolate.pyx similarity index 100% rename from xradar/io/backends/nexrad_interpolate.pyx rename to xradar/interpolate/nexrad_interpolate.pyx diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index 7a79d8d..e7af09c 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -11,6 +11,10 @@ from xarray.core.variable import Variable from xradar import util +from xradar.interpolate.nexrad_interpolate import ( + _fast_interpolate_scan_2, + _fast_interpolate_scan_4, +) from xradar.io.backends.common import ( LazyLoadDict, _assign_root, @@ -29,11 +33,6 @@ get_time_attrs, ) -from .nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) - nexrad_mapping = { "REF": "DBZH", "VEL": "VRADH", From 78fada349b97a4715fd65ef03109ba7845f1acb1 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:03:40 -0600 Subject: [PATCH 22/48] ADD: Update setup --- setup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 3c60099..b0bf9fa 100644 --- a/setup.py +++ b/setup.py @@ -10,8 +10,8 @@ # Modules to be compiled and include_dirs when necessary extensions = [ Extension( - "xradar.io.backends.nexrad_interpolate", - sources=["xradar/io/backends/nexrad_interpolate.pyx"], + "xradar.interpolate.nexrad_interpolate", + sources=["xradar/interpolate/nexrad_interpolate.pyx"], include_dirs=[numpy.get_include()], ), ] From cb9435001913a4b9d21ba7ff02d5857f9f100b78 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:08:44 -0600 Subject: [PATCH 23/48] Move to hidden module --- setup.py | 4 +- xradar/interpolate/_nexrad_interpolate.pyx | 110 +++++++++++++++++++++ xradar/io/backends/nexrad_level2.py | 2 +- 3 files changed, 113 insertions(+), 3 deletions(-) create mode 100644 xradar/interpolate/_nexrad_interpolate.pyx diff --git a/setup.py b/setup.py index b0bf9fa..2de08f2 100644 --- a/setup.py +++ b/setup.py @@ -10,8 +10,8 @@ # Modules to be compiled and include_dirs when necessary extensions = [ Extension( - "xradar.interpolate.nexrad_interpolate", - sources=["xradar/interpolate/nexrad_interpolate.pyx"], + "xradar.interpolate._nexrad_interpolate", + sources=["xradar/interpolate/_nexrad_interpolate.pyx"], include_dirs=[numpy.get_include()], ), ] diff --git a/xradar/interpolate/_nexrad_interpolate.pyx b/xradar/interpolate/_nexrad_interpolate.pyx new file mode 100644 index 0000000..2930d29 --- /dev/null +++ b/xradar/interpolate/_nexrad_interpolate.pyx @@ -0,0 +1,110 @@ +""" +Interpolation of NEXRAD moments from 1000 meter to 250 meter gate spacing. + +""" + +def _fast_interpolate_scan_4( + float[:, :] data, float[:] scratch_ray, float fill_value, + int start, int end, int moment_ngates, int linear_interp): + """ Interpolate a single NEXRAD moment scan from 1000 m to 250 m. """ + # This interpolation scheme is only valid for NEXRAD data where a 4:1 + # (1000 m : 250 m) interpolation is needed. + # + # The scheme here performs a linear interpolation between pairs of gates + # in a ray when the both of the gates are not masked (below threshold). + # When one of the gates is masked the interpolation changes to a nearest + # neighbor interpolation. Nearest neighbor is also performed at the end + # points until the new range bin would be centered beyond half of the range + # spacing of the original range. + # + # Nearest neighbor interpolation is performed when linear_interp is False, + # this is equivalent to repeating each gate four times in each ray. + # + # No transformation of the raw data is performed prior to interpolation, so + # reflectivity will be interpolated in dB units, velocity in m/s, etc, + # this may not be the best method for interpolation. + # + # This method was adapted from Radx and Py-ART + cdef int ray_num, i, interp_ngates + cdef float gate_val, next_val, delta + + interp_ngates = 4 * moment_ngates # number of gates interpolated + + for ray_num in range(start, end+1): + + # repeat each gate value 4 times + for i in range(moment_ngates): + gate_val = data[ray_num, i] + scratch_ray[i*4 + 0] = gate_val + scratch_ray[i*4 + 1] = gate_val + scratch_ray[i*4 + 2] = gate_val + scratch_ray[i*4 + 3] = gate_val + + if linear_interp: + # linear interpolate + for i in range(2, interp_ngates - 4, 4): + gate_val = scratch_ray[i] + next_val = scratch_ray[i+4] + if gate_val == fill_value or next_val == fill_value: + continue + delta = (next_val - gate_val) / 4. + scratch_ray[i+0] = gate_val + delta * 0.5 + scratch_ray[i+1] = gate_val + delta * 1.5 + scratch_ray[i+2] = gate_val + delta * 2.5 + scratch_ray[i+3] = gate_val + delta * 3.5 + + for i in range(interp_ngates): + data[ray_num, i] = scratch_ray[i] + + +def _fast_interpolate_scan_2( + float[:, :] data, float[:] scratch_ray, float fill_value, + int start, int end, int moment_ngates, int linear_interp): + """ Interpolate a single NEXRAD moment scan from 300 m to 150 m. """ + # This interpolation scheme is only valid for NEXRAD TWDR data where a 2:1 + # (300 m : 150 m) interpolation is needed. + # + # The scheme here performs a linear interpolation between pairs of gates + # in a ray when the both of the gates are not masked (below threshold). + # When one of the gates is masked the interpolation changes to a nearest + # neighbor interpolation. Nearest neighbor is also performed at the end + # points until the new range bin would be centered beyond half of the range + # spacing of the original range. + # + # Nearest neighbor interpolation is performed when linear_interp is False, + # this is equivalent to repeating each gate four times in each ray. + # + # No transformation of the raw data is performed prior to interpolation, so + # reflectivity will be interpolated in dB units, velocity in m/s, etc, + # this may not be the best method for interpolation. + # + # This method was adapted from Radx + cdef int ray_num, i, interp_ngates + cdef float gate_val, next_val, delta + + interp_ngates = 2 * moment_ngates - 1 # number of gates interpolated + + for ray_num in range(start, end+1): + + # repeat each gate value 4 times + for i in range(moment_ngates): + gate_val = data[ray_num, i] + if i == moment_ngates - 1: + scratch_ray[i*2 + 0] = gate_val + else: + scratch_ray[i*2 + 0] = gate_val + scratch_ray[i*2 + 1] = gate_val + + if linear_interp: + # linear interpolate + for i in range(1, interp_ngates - 2, 2): + gate_val = scratch_ray[i] + next_val = scratch_ray[i+2] + if gate_val == fill_value or next_val == fill_value: + continue + delta = (next_val - gate_val) / 2. + scratch_ray[i+0] = gate_val + delta * 0.5 + scratch_ray[i+1] = gate_val + delta * 1.5 + + for i in range(interp_ngates): + data[ray_num, i] = scratch_ray[i] diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index e7af09c..eebd901 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -11,7 +11,7 @@ from xarray.core.variable import Variable from xradar import util -from xradar.interpolate.nexrad_interpolate import ( +from xradar.interpolate._nexrad_interpolate import ( _fast_interpolate_scan_2, _fast_interpolate_scan_4, ) From c3824d64b81838aa0b5bc44dc6a40a23f6780024 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:14:45 -0600 Subject: [PATCH 24/48] FIX: Fix all imports in backends --- xradar/io/backends/__init__.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 81b0278..18fbbad 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -29,5 +29,3 @@ NexradLevel2BackendEntrypoint, # noqa open_nexradlevel2_datatree, # noqa ) - -__all__ = [s for s in dir() if not s.startswith("_")] From 5e0da73a14c0118cbd866cece7977e2ac2efbb85 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:19:24 -0600 Subject: [PATCH 25/48] fix manifest --- MANIFEST.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MANIFEST.in b/MANIFEST.in index 89afb7b..109f8ea 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -5,7 +5,7 @@ include LICENSE include README.md recursive-include tests * -recursive-include xradar *.c *.pyx *.pxd recursive-exclude * __pycache__ +recursive-exclude * *.py[co] recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif From d16e2e3fab82cb9e1333fc6bcdc4b597956d5e4d Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:20:02 -0600 Subject: [PATCH 26/48] include all cython --- MANIFEST.in | 1 + 1 file changed, 1 insertion(+) diff --git a/MANIFEST.in b/MANIFEST.in index 109f8ea..2473f4a 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -8,4 +8,5 @@ recursive-include tests * recursive-exclude * __pycache__ recursive-exclude * *.py[co] +global-include *.pyx recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif From 3a0a76f5710e2e044b6019fe8e7627c62e61e483 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:20:48 -0600 Subject: [PATCH 27/48] Add other types of cython files --- MANIFEST.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MANIFEST.in b/MANIFEST.in index 2473f4a..030399f 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -8,5 +8,5 @@ recursive-include tests * recursive-exclude * __pycache__ recursive-exclude * *.py[co] -global-include *.pyx +global-include *.pyx *pxd recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif From 62db2598f9b2639b0a5d6d25d09427ec577e28ba Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:35:33 -0600 Subject: [PATCH 28/48] debug submodule --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index fb44926..f98c456 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -59,6 +59,7 @@ jobs: - name: Install xradar run: | python -m pip install . --force-reinstall + ls xradar/interpolate/ - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" From 2edb91cfa94bdfdb278bc6328e6df751658440d3 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:47:02 -0600 Subject: [PATCH 29/48] ADD: add this back in --- xradar/interpolate/__init__.pxd | 11 +++ xradar/interpolate/nexrad_interpolate.pyx | 110 ---------------------- 2 files changed, 11 insertions(+), 110 deletions(-) create mode 100644 xradar/interpolate/__init__.pxd delete mode 100644 xradar/interpolate/nexrad_interpolate.pyx diff --git a/xradar/interpolate/__init__.pxd b/xradar/interpolate/__init__.pxd new file mode 100644 index 0000000..7e598c8 --- /dev/null +++ b/xradar/interpolate/__init__.pxd @@ -0,0 +1,11 @@ +#!/usr/bin/env python +# Copyright (c) 2024, openradar developers. +# Distributed under the MIT License. See LICENSE for more info. + +""" +XRadar Interpolation +==================== + +""" + +from xradar.interpolate._nexrad_interpolate cimport * # noqa diff --git a/xradar/interpolate/nexrad_interpolate.pyx b/xradar/interpolate/nexrad_interpolate.pyx deleted file mode 100644 index 2930d29..0000000 --- a/xradar/interpolate/nexrad_interpolate.pyx +++ /dev/null @@ -1,110 +0,0 @@ -""" -Interpolation of NEXRAD moments from 1000 meter to 250 meter gate spacing. - -""" - -def _fast_interpolate_scan_4( - float[:, :] data, float[:] scratch_ray, float fill_value, - int start, int end, int moment_ngates, int linear_interp): - """ Interpolate a single NEXRAD moment scan from 1000 m to 250 m. """ - # This interpolation scheme is only valid for NEXRAD data where a 4:1 - # (1000 m : 250 m) interpolation is needed. - # - # The scheme here performs a linear interpolation between pairs of gates - # in a ray when the both of the gates are not masked (below threshold). - # When one of the gates is masked the interpolation changes to a nearest - # neighbor interpolation. Nearest neighbor is also performed at the end - # points until the new range bin would be centered beyond half of the range - # spacing of the original range. - # - # Nearest neighbor interpolation is performed when linear_interp is False, - # this is equivalent to repeating each gate four times in each ray. - # - # No transformation of the raw data is performed prior to interpolation, so - # reflectivity will be interpolated in dB units, velocity in m/s, etc, - # this may not be the best method for interpolation. - # - # This method was adapted from Radx and Py-ART - cdef int ray_num, i, interp_ngates - cdef float gate_val, next_val, delta - - interp_ngates = 4 * moment_ngates # number of gates interpolated - - for ray_num in range(start, end+1): - - # repeat each gate value 4 times - for i in range(moment_ngates): - gate_val = data[ray_num, i] - scratch_ray[i*4 + 0] = gate_val - scratch_ray[i*4 + 1] = gate_val - scratch_ray[i*4 + 2] = gate_val - scratch_ray[i*4 + 3] = gate_val - - if linear_interp: - # linear interpolate - for i in range(2, interp_ngates - 4, 4): - gate_val = scratch_ray[i] - next_val = scratch_ray[i+4] - if gate_val == fill_value or next_val == fill_value: - continue - delta = (next_val - gate_val) / 4. - scratch_ray[i+0] = gate_val + delta * 0.5 - scratch_ray[i+1] = gate_val + delta * 1.5 - scratch_ray[i+2] = gate_val + delta * 2.5 - scratch_ray[i+3] = gate_val + delta * 3.5 - - for i in range(interp_ngates): - data[ray_num, i] = scratch_ray[i] - - -def _fast_interpolate_scan_2( - float[:, :] data, float[:] scratch_ray, float fill_value, - int start, int end, int moment_ngates, int linear_interp): - """ Interpolate a single NEXRAD moment scan from 300 m to 150 m. """ - # This interpolation scheme is only valid for NEXRAD TWDR data where a 2:1 - # (300 m : 150 m) interpolation is needed. - # - # The scheme here performs a linear interpolation between pairs of gates - # in a ray when the both of the gates are not masked (below threshold). - # When one of the gates is masked the interpolation changes to a nearest - # neighbor interpolation. Nearest neighbor is also performed at the end - # points until the new range bin would be centered beyond half of the range - # spacing of the original range. - # - # Nearest neighbor interpolation is performed when linear_interp is False, - # this is equivalent to repeating each gate four times in each ray. - # - # No transformation of the raw data is performed prior to interpolation, so - # reflectivity will be interpolated in dB units, velocity in m/s, etc, - # this may not be the best method for interpolation. - # - # This method was adapted from Radx - cdef int ray_num, i, interp_ngates - cdef float gate_val, next_val, delta - - interp_ngates = 2 * moment_ngates - 1 # number of gates interpolated - - for ray_num in range(start, end+1): - - # repeat each gate value 4 times - for i in range(moment_ngates): - gate_val = data[ray_num, i] - if i == moment_ngates - 1: - scratch_ray[i*2 + 0] = gate_val - else: - scratch_ray[i*2 + 0] = gate_val - scratch_ray[i*2 + 1] = gate_val - - if linear_interp: - # linear interpolate - for i in range(1, interp_ngates - 2, 2): - gate_val = scratch_ray[i] - next_val = scratch_ray[i+2] - if gate_val == fill_value or next_val == fill_value: - continue - delta = (next_val - gate_val) / 2. - scratch_ray[i+0] = gate_val + delta * 0.5 - scratch_ray[i+1] = gate_val + delta * 1.5 - - for i in range(interp_ngates): - data[ray_num, i] = scratch_ray[i] From f0ec0e7237c96c7de6befd11f17d33217085c87d Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 15:58:57 -0600 Subject: [PATCH 30/48] Be explicit about submodules --- xradar/io/__init__.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/xradar/io/__init__.py b/xradar/io/__init__.py index 560175e..061cad1 100644 --- a/xradar/io/__init__.py +++ b/xradar/io/__init__.py @@ -13,7 +13,24 @@ .. automodule:: xradar.io.export """ -from .backends import * # noqa +from .backends import ( + CfRadial1BackendEntrypoint, # noqa + FurunoBackendEntrypoint, # noqa + GamicBackendEntrypoint, # noqa + IrisBackendEntrypoint, # noqa + NexradLevel2BackendEntrypoint, # noqa + OdimBackendEntrypoint, # noqa + RainbowBackendEntrypoint, # noqa + open_cfradial1_datatree, # noqa + open_furuno_datatree, # noqa + open_gamic_datatree, # noqa + open_iris_datatree, # noqa + open_nexradlevel2_datatree, # noqa + open_odim_datatree, # noqa + open_rainbow_datatree, # noqa +) + +# noqa from .export import * # noqa __all__ = [s for s in dir() if not s.startswith("_")] From f1af9cd4b02d7f9a3c736634d00a3296d5e97d0f Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:01:22 -0600 Subject: [PATCH 31/48] ADD: Add clear submodule --- xradar/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xradar/__init__.py b/xradar/__init__.py index f87acfb..30d7c38 100644 --- a/xradar/__init__.py +++ b/xradar/__init__.py @@ -20,6 +20,7 @@ # import subpackages from . import accessors # noqa from . import georeference # noqa +from . import interpolate # noqa from . import io # noqa from . import model # noqa from . import util # noqa From 6c7b824ce8605a7d1154d1517afd05bb0c952017 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:07:57 -0600 Subject: [PATCH 32/48] ADD: Add imports --- xradar/interpolate/__init__.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/xradar/interpolate/__init__.py b/xradar/interpolate/__init__.py index 6e1a1c9..a53824b 100644 --- a/xradar/interpolate/__init__.py +++ b/xradar/interpolate/__init__.py @@ -5,7 +5,11 @@ """ XRadar Interpolation ==================== - """ +from ._nexrad_interpolate import ( + _fast_interpolate_scan_2, # noqa + _fast_interpolate_scan_4, # noqa +) + __all__ = [s for s in dir() if not s.startswith("_")] From f162928051a1431846c3b0673e4e7e16b67c5656 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:17:15 -0600 Subject: [PATCH 33/48] remove name in setup --- setup.py | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.py b/setup.py index 2de08f2..63dedc1 100644 --- a/setup.py +++ b/setup.py @@ -19,7 +19,6 @@ # This is the function that is executed setup( - name="xradar", # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From 9592148496e77769f70bde7c5308a46ee90e6767 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:21:10 -0600 Subject: [PATCH 34/48] remove check on version --- .github/workflows/ci.yml | 3 --- 1 file changed, 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f98c456..cf60010 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -60,9 +60,6 @@ jobs: run: | python -m pip install . --force-reinstall ls xradar/interpolate/ - - name: Version Info - run: | - python -c "import xradar; print(xradar.version.version)" - name: Test with pytest run: | pytest -n auto --dist loadfile --verbose --durations=15 --cov-report xml:coverage_unit.xml --cov=xradar --pyargs tests From 045ba9962cd8e20a403d43ea99349a72128273a2 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:28:03 -0600 Subject: [PATCH 35/48] make sure submodule is blank --- xradar/interpolate/__init__.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/xradar/interpolate/__init__.py b/xradar/interpolate/__init__.py index a53824b..70e6149 100644 --- a/xradar/interpolate/__init__.py +++ b/xradar/interpolate/__init__.py @@ -6,10 +6,3 @@ XRadar Interpolation ==================== """ - -from ._nexrad_interpolate import ( - _fast_interpolate_scan_2, # noqa - _fast_interpolate_scan_4, # noqa -) - -__all__ = [s for s in dir() if not s.startswith("_")] From 2f265658354be95be9d1e5534ab07a85a13efe33 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:29:09 -0600 Subject: [PATCH 36/48] ADD: Add extra line --- setup.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 63dedc1..3b4eb5a 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,3 @@ -import numpy from Cython.Build import cythonize from Cython.Compiler import Options from setuptools import Extension, setup @@ -12,13 +11,13 @@ Extension( "xradar.interpolate._nexrad_interpolate", sources=["xradar/interpolate/_nexrad_interpolate.pyx"], - include_dirs=[numpy.get_include()], ), ] # This is the function that is executed setup( + install_requires=[""], # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From d7d65b206924417a6f80112c46c004898cded823 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 16:30:22 -0600 Subject: [PATCH 37/48] be more explicit --- setup.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 3b4eb5a..d983abe 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,13 @@ # This is the function that is executed setup( - install_requires=[""], + install_requires=[ + "setuptools>=45", + "wheel", + "setuptools_scm[toml]>=7.0", + "cython", + "numpy", + ], # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From 3b49bc3584676fd94c9fccbf2dad2e9c45f027f1 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 19:20:34 -0600 Subject: [PATCH 38/48] add setup.py run --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index cf60010..ebc17c5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -59,7 +59,7 @@ jobs: - name: Install xradar run: | python -m pip install . --force-reinstall - ls xradar/interpolate/ + python setup.py build_ext --inplace - name: Test with pytest run: | pytest -n auto --dist loadfile --verbose --durations=15 --cov-report xml:coverage_unit.xml --cov=xradar --pyargs tests From 9c11f0b56e5ff9617d2f4afeeb00ab7742af0f96 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 19:38:14 -0600 Subject: [PATCH 39/48] ADD: Add proper installation to other parts --- .github/workflows/ci.yml | 6 +++++- .github/workflows/upload_pypi.yml | 1 + setup.py | 7 ------- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ebc17c5..11c1048 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,8 +58,11 @@ jobs: conda - name: Install xradar run: | - python -m pip install . --force-reinstall + python -m pip install . --no-deps python setup.py build_ext --inplace + - name: Version Info + run: | + python -c "import xradar; print(xradar.version.version)" - name: Test with pytest run: | pytest -n auto --dist loadfile --verbose --durations=15 --cov-report xml:coverage_unit.xml --cov=xradar --pyargs tests @@ -102,6 +105,7 @@ jobs: - name: Install xradar run: | python -m pip install . --no-deps + python setup.py build_ext --inplace - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" diff --git a/.github/workflows/upload_pypi.yml b/.github/workflows/upload_pypi.yml index 0a14fad..ce0d5f6 100644 --- a/.github/workflows/upload_pypi.yml +++ b/.github/workflows/upload_pypi.yml @@ -30,4 +30,5 @@ jobs: TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }} run: | python -m build + python setup.py build_ext --inplace twine upload dist/* diff --git a/setup.py b/setup.py index d983abe..c2bce61 100644 --- a/setup.py +++ b/setup.py @@ -17,13 +17,6 @@ # This is the function that is executed setup( - install_requires=[ - "setuptools>=45", - "wheel", - "setuptools_scm[toml]>=7.0", - "cython", - "numpy", - ], # external to be compiled ext_modules=cythonize( extensions, compiler_directives={"language_level": "3", "cpow": True} From e1a59dfd1228e5f9b90e93cd0857267ef6a01164 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 20:37:26 -0600 Subject: [PATCH 40/48] ADD: add test for lazy dict --- tests/io/backends/test_common.py | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 tests/io/backends/test_common.py diff --git a/tests/io/backends/test_common.py b/tests/io/backends/test_common.py new file mode 100644 index 0000000..0a0afb6 --- /dev/null +++ b/tests/io/backends/test_common.py @@ -0,0 +1,9 @@ +from xradar.io.backends import common + + +def test_lazy_dict(): + d = common.LazyLoadDict({"key1": "value1", "key2": "value2"}) + assert d["key1"] == "value1" + lazy_func = lambda: 999 + d.set_lazy("lazykey1", lazy_func) + assert d["lazykey1"] == 999 From 98090adbe737d06cfddc34a4f09128c15ea8798c Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 25 Jan 2024 20:49:19 -0600 Subject: [PATCH 41/48] DEL: Remove unused sections of common --- xradar/io/backends/common.py | 42 ------------------------------------ 1 file changed, 42 deletions(-) diff --git a/xradar/io/backends/common.py b/xradar/io/backends/common.py index b338498..32a36ea 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -22,7 +22,6 @@ import fsspec import h5netcdf -import netCDF4 import numpy as np import xarray as xr from datatree import DataTree @@ -281,47 +280,6 @@ def prepare_for_read(filename, storage_options={"anon": True}): ).open() -def stringarray_to_chararray(arr, numchars=None): - """ - Convert an string array to a character array with one extra dimension. - - Parameters - ---------- - arr : array - Array with numpy dtype 'SN', where N is the number of characters - in the string. - - numchars : int - Number of characters used to represent the string. If numchar > N - the results will be padded on the right with blanks. The default, - None will use N. - - Returns - ------- - chararr : array - Array with dtype 'S1' and shape = arr.shape + (numchars, ). - - """ - carr = netCDF4.stringtochar(arr) - if numchars is None: - return carr - - arr_numchars = carr.shape[-1] - if numchars <= arr_numchars: - raise ValueError("numchars must be >= %i" % (arr_numchars)) - chararr = np.zeros(arr.shape + (numchars,), dtype="S1") - chararr[..., :arr_numchars] = carr[:] - return chararr - - -def _test_arguments(dic): - """Issue a warning if receive non-empty argument dict.""" - if dic: - import warnings - - warnings.warn("Unexpected arguments: %s" % dic.keys()) - - def make_time_unit_str(dtobj): """Return a time unit string from a datetime object.""" return "seconds since " + dtobj.strftime("%Y-%m-%dT%H:%M:%SZ") From 528a16d469c463b4de3f86906afc8e91682019e2 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 14:46:36 -0500 Subject: [PATCH 42/48] DEL: Remove the interpolation step --- ci/notebooktests.yml | 1 - ci/unittests.yml | 1 - xradar/interpolate/__init__.py | 8 -- xradar/interpolate/_nexrad_interpolate.pyx | 110 --------------------- xradar/io/backends/__init__.py | 2 + xradar/io/backends/nexrad_level2.py | 79 +-------------- 6 files changed, 3 insertions(+), 198 deletions(-) delete mode 100644 xradar/interpolate/__init__.py delete mode 100644 xradar/interpolate/_nexrad_interpolate.pyx diff --git a/ci/notebooktests.yml b/ci/notebooktests.yml index e0628db..a0587a0 100644 --- a/ci/notebooktests.yml +++ b/ci/notebooktests.yml @@ -7,7 +7,6 @@ dependencies: - cmweather - codecov - coverage - - cython - dask - h5netcdf - h5py diff --git a/ci/unittests.yml b/ci/unittests.yml index d345200..9291f8e 100644 --- a/ci/unittests.yml +++ b/ci/unittests.yml @@ -6,7 +6,6 @@ dependencies: - codecov - cmweather - coverage - - cython - dask - h5netcdf - h5py diff --git a/xradar/interpolate/__init__.py b/xradar/interpolate/__init__.py deleted file mode 100644 index 70e6149..0000000 --- a/xradar/interpolate/__init__.py +++ /dev/null @@ -1,8 +0,0 @@ -#!/usr/bin/env python -# Copyright (c) 2024, openradar developers. -# Distributed under the MIT License. See LICENSE for more info. - -""" -XRadar Interpolation -==================== -""" diff --git a/xradar/interpolate/_nexrad_interpolate.pyx b/xradar/interpolate/_nexrad_interpolate.pyx deleted file mode 100644 index 2930d29..0000000 --- a/xradar/interpolate/_nexrad_interpolate.pyx +++ /dev/null @@ -1,110 +0,0 @@ -""" -Interpolation of NEXRAD moments from 1000 meter to 250 meter gate spacing. - -""" - -def _fast_interpolate_scan_4( - float[:, :] data, float[:] scratch_ray, float fill_value, - int start, int end, int moment_ngates, int linear_interp): - """ Interpolate a single NEXRAD moment scan from 1000 m to 250 m. """ - # This interpolation scheme is only valid for NEXRAD data where a 4:1 - # (1000 m : 250 m) interpolation is needed. - # - # The scheme here performs a linear interpolation between pairs of gates - # in a ray when the both of the gates are not masked (below threshold). - # When one of the gates is masked the interpolation changes to a nearest - # neighbor interpolation. Nearest neighbor is also performed at the end - # points until the new range bin would be centered beyond half of the range - # spacing of the original range. - # - # Nearest neighbor interpolation is performed when linear_interp is False, - # this is equivalent to repeating each gate four times in each ray. - # - # No transformation of the raw data is performed prior to interpolation, so - # reflectivity will be interpolated in dB units, velocity in m/s, etc, - # this may not be the best method for interpolation. - # - # This method was adapted from Radx and Py-ART - cdef int ray_num, i, interp_ngates - cdef float gate_val, next_val, delta - - interp_ngates = 4 * moment_ngates # number of gates interpolated - - for ray_num in range(start, end+1): - - # repeat each gate value 4 times - for i in range(moment_ngates): - gate_val = data[ray_num, i] - scratch_ray[i*4 + 0] = gate_val - scratch_ray[i*4 + 1] = gate_val - scratch_ray[i*4 + 2] = gate_val - scratch_ray[i*4 + 3] = gate_val - - if linear_interp: - # linear interpolate - for i in range(2, interp_ngates - 4, 4): - gate_val = scratch_ray[i] - next_val = scratch_ray[i+4] - if gate_val == fill_value or next_val == fill_value: - continue - delta = (next_val - gate_val) / 4. - scratch_ray[i+0] = gate_val + delta * 0.5 - scratch_ray[i+1] = gate_val + delta * 1.5 - scratch_ray[i+2] = gate_val + delta * 2.5 - scratch_ray[i+3] = gate_val + delta * 3.5 - - for i in range(interp_ngates): - data[ray_num, i] = scratch_ray[i] - - -def _fast_interpolate_scan_2( - float[:, :] data, float[:] scratch_ray, float fill_value, - int start, int end, int moment_ngates, int linear_interp): - """ Interpolate a single NEXRAD moment scan from 300 m to 150 m. """ - # This interpolation scheme is only valid for NEXRAD TWDR data where a 2:1 - # (300 m : 150 m) interpolation is needed. - # - # The scheme here performs a linear interpolation between pairs of gates - # in a ray when the both of the gates are not masked (below threshold). - # When one of the gates is masked the interpolation changes to a nearest - # neighbor interpolation. Nearest neighbor is also performed at the end - # points until the new range bin would be centered beyond half of the range - # spacing of the original range. - # - # Nearest neighbor interpolation is performed when linear_interp is False, - # this is equivalent to repeating each gate four times in each ray. - # - # No transformation of the raw data is performed prior to interpolation, so - # reflectivity will be interpolated in dB units, velocity in m/s, etc, - # this may not be the best method for interpolation. - # - # This method was adapted from Radx - cdef int ray_num, i, interp_ngates - cdef float gate_val, next_val, delta - - interp_ngates = 2 * moment_ngates - 1 # number of gates interpolated - - for ray_num in range(start, end+1): - - # repeat each gate value 4 times - for i in range(moment_ngates): - gate_val = data[ray_num, i] - if i == moment_ngates - 1: - scratch_ray[i*2 + 0] = gate_val - else: - scratch_ray[i*2 + 0] = gate_val - scratch_ray[i*2 + 1] = gate_val - - if linear_interp: - # linear interpolate - for i in range(1, interp_ngates - 2, 2): - gate_val = scratch_ray[i] - next_val = scratch_ray[i+2] - if gate_val == fill_value or next_val == fill_value: - continue - delta = (next_val - gate_val) / 2. - scratch_ray[i+0] = gate_val + delta * 0.5 - scratch_ray[i+1] = gate_val + delta * 1.5 - - for i in range(interp_ngates): - data[ray_num, i] = scratch_ray[i] diff --git a/xradar/io/backends/__init__.py b/xradar/io/backends/__init__.py index 18fbbad..81b0278 100644 --- a/xradar/io/backends/__init__.py +++ b/xradar/io/backends/__init__.py @@ -29,3 +29,5 @@ NexradLevel2BackendEntrypoint, # noqa open_nexradlevel2_datatree, # noqa ) + +__all__ = [s for s in dir() if not s.startswith("_")] diff --git a/xradar/io/backends/nexrad_level2.py b/xradar/io/backends/nexrad_level2.py index eebd901..c4ae52b 100644 --- a/xradar/io/backends/nexrad_level2.py +++ b/xradar/io/backends/nexrad_level2.py @@ -11,10 +11,6 @@ from xarray.core.variable import Variable from xradar import util -from xradar.interpolate._nexrad_interpolate import ( - _fast_interpolate_scan_2, - _fast_interpolate_scan_4, -) from xradar.io.backends.common import ( LazyLoadDict, _assign_root, @@ -1007,55 +1003,11 @@ def _find_range_params(scan_info): return min_first_gate, min_gate_spacing, max_last_gate -def _find_scans_to_interp(scan_info, first_gate, gate_spacing): - """Return a dict indicating what moments/scans need interpolation.""" - moments = {m for scan in scan_info for m in scan["moments"]} - interpolate = {moment: [] for moment in moments} - for scan_num, scan in enumerate(scan_info): - for moment in moments: - if moment not in scan["moments"]: - continue - index = scan["moments"].index(moment) - first = scan["first_gate"][index] - spacing = scan["gate_spacing"][index] - if first != first_gate or spacing != gate_spacing: - interpolate[moment].append(scan_num) - # for proper interpolation the gate spacing of the scan to be - # interpolated should be 1/4th the spacing of the radar - if spacing == gate_spacing * 4: - interpolate["multiplier"] = "4" - elif spacing == gate_spacing * 2: - interpolate["multiplier"] = "2" - else: - raise ValueError("Gate spacing is neither 1/4 or 1/2") - # assert first_gate + 1.5 * gate_spacing == first - # remove moments with no scans needing interpolation - interpolate = {k: v for k, v in interpolate.items() if len(v) != 0} - return interpolate - - -def _interpolate_scan(mdata, start, end, moment_ngates, multiplier, linear_interp=True): - """Interpolate a single NEXRAD moment scan from 1000 m to 250 m.""" - fill_value = -9999 - data = mdata.filled(fill_value) - scratch_ray = np.empty((data.shape[1],), dtype=data.dtype) - if multiplier == "4": - _fast_interpolate_scan_4( - data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp - ) - else: - _fast_interpolate_scan_2( - data, scratch_ray, fill_value, start, end, moment_ngates, linear_interp - ) - mdata[:] = np.ma.array(data, mask=(data == fill_value)) - - def format_nexrad_level2_data( filename, delay_field_loading=False, station=None, scans=None, - linear_interp=True, storage_options={"anon": True}, first_dimension=None, group=None, @@ -1069,10 +1021,6 @@ def format_nexrad_level2_data( # Access the scan info and load in the available moments scan_info = nfile.scan_info(scans) available_moments = {m for scan in scan_info for m in scan["moments"]} - first_gate, gate_spacing, last_gate = _find_range_params(scan_info) - - # Interpolate to 360 degrees where neccessary - interpolate = _find_scans_to_interp(scan_info, first_gate, gate_spacing) # Deal with time time_start, _time = nfile.get_times(scans) @@ -1189,42 +1137,17 @@ def format_nexrad_level2_data( # fields max_ngates = len(_range["data"]) - available_moments = {m for scan in scan_info for m in scan["moments"]} - interpolate = _find_scans_to_interp( - scan_info, - first_gate, - gate_spacing, - ) fields = {} for moment in available_moments: dic = get_moment_attrs(nexrad_mapping[moment]) dic["_FillValue"] = -9999 - if delay_field_loading and moment not in interpolate: + if delay_field_loading: dic = LazyLoadDict(dic) data_call = _NEXRADLevel2StagedField(nfile, moment, max_ngates, scans) dic.set_lazy("data", data_call) else: mdata = nfile.get_data(moment, max_ngates, scans=scans) - if moment in interpolate: - interp_scans = interpolate[moment] - warnings.warn( - "Gate spacing is not constant, interpolating data in " - + f"scans {interp_scans} for moment {moment}.", - UserWarning, - ) - for scan in interp_scans: - idx = scan_info[scan]["moments"].index(moment) - moment_ngates = scan_info[scan]["ngates"][idx] - start = sweep_start_ray_index["data"][scan] - end = sweep_end_ray_index["data"][scan] - if interpolate["multiplier"] == "4": - multiplier = "4" - else: - multiplier = "2" - _interpolate_scan( - mdata, start, end, moment_ngates, multiplier, linear_interp - ) dic["data"] = mdata fields[nexrad_mapping[moment]] = dic instrument_parameters = {} From 1b0d97bbce32aab726b5f044d447bc2bc0f9ae0e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 14:49:50 -0500 Subject: [PATCH 43/48] FIX: Fix linting --- setup.py | 24 ------------------------ xradar/interpolate/__init__.pxd | 11 ----------- 2 files changed, 35 deletions(-) delete mode 100644 setup.py delete mode 100644 xradar/interpolate/__init__.pxd diff --git a/setup.py b/setup.py deleted file mode 100644 index c2bce61..0000000 --- a/setup.py +++ /dev/null @@ -1,24 +0,0 @@ -from Cython.Build import cythonize -from Cython.Compiler import Options -from setuptools import Extension, setup - -# These are optional -Options.docstrings = True -Options.annotate = False - -# Modules to be compiled and include_dirs when necessary -extensions = [ - Extension( - "xradar.interpolate._nexrad_interpolate", - sources=["xradar/interpolate/_nexrad_interpolate.pyx"], - ), -] - - -# This is the function that is executed -setup( - # external to be compiled - ext_modules=cythonize( - extensions, compiler_directives={"language_level": "3", "cpow": True} - ), -) diff --git a/xradar/interpolate/__init__.pxd b/xradar/interpolate/__init__.pxd deleted file mode 100644 index 7e598c8..0000000 --- a/xradar/interpolate/__init__.pxd +++ /dev/null @@ -1,11 +0,0 @@ -#!/usr/bin/env python -# Copyright (c) 2024, openradar developers. -# Distributed under the MIT License. See LICENSE for more info. - -""" -XRadar Interpolation -==================== - -""" - -from xradar.interpolate._nexrad_interpolate cimport * # noqa From bae3d194802101bfa93d77299450c064921d040e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 15:05:40 -0500 Subject: [PATCH 44/48] Only use ruff for linting --- .github/workflows/ci.yml | 6 +----- xradar/io/backends/nexrad_common.py | 1 + 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 11c1048..8fb95f7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,9 +22,6 @@ jobs: run: | python -m pip install --upgrade pip pip install black black[jupyter] ruff - - name: Black style check - run: | - black --check . - name: Lint with ruff run: | ruff . @@ -104,8 +101,7 @@ jobs: conda - name: Install xradar run: | - python -m pip install . --no-deps - python setup.py build_ext --inplace + python setup.py build_ext --inplace --no-deps - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" diff --git a/xradar/io/backends/nexrad_common.py b/xradar/io/backends/nexrad_common.py index e34c277..13f3c6a 100644 --- a/xradar/io/backends/nexrad_common.py +++ b/xradar/io/backends/nexrad_common.py @@ -2,6 +2,7 @@ Data and functions common to all types of NEXRAD files. """ + # The functions in this module are intended to be used in other # nexrad related modules. The functions are not and should not be # exported into the pyart.io namespace. From aee82ed8a8b5bc26435e261c4965989b2fc10d75 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 15:07:58 -0500 Subject: [PATCH 45/48] DOC: Add addition to history file --- docs/history.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/history.md b/docs/history.md index ad84b90..b85dc5e 100644 --- a/docs/history.md +++ b/docs/history.md @@ -1,5 +1,8 @@ # History +## Development Version +* ENH: Add a reader for nexrad level2 files ({pull}`147`) by [@mgrover1](https://github.com/mgrover1) + ## 0.4.2 (2023-11-02) * FIX: Fix handling of sweep_mode attributes ({pull}`143`) by [@mgrover1](https://github.com/mgrover1) From fc03120f51281f2549fddc51e014f676786fde62 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 15:12:43 -0500 Subject: [PATCH 46/48] ADD: add original ci back in --- .github/workflows/ci.yml | 3 +-- pyproject.toml | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8fb95f7..246d77c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -56,7 +56,6 @@ jobs: - name: Install xradar run: | python -m pip install . --no-deps - python setup.py build_ext --inplace - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" @@ -101,7 +100,7 @@ jobs: conda - name: Install xradar run: | - python setup.py build_ext --inplace --no-deps + python -m pip install . --no-deps - name: Version Info run: | python -c "import xradar; print(xradar.version.version)" diff --git a/pyproject.toml b/pyproject.toml index 9ac1069..85665f7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,6 @@ requires = [ "setuptools>=45", "wheel", "setuptools_scm[toml]>=7.0", - "cython", "numpy" ] build-backend = "setuptools.build_meta" From c33d52153dd5507c7fc56aee27a761678370ed69 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 30 Jan 2024 17:29:46 -0500 Subject: [PATCH 47/48] DEL: Delete extra init --- xradar/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/xradar/__init__.py b/xradar/__init__.py index 30d7c38..f87acfb 100644 --- a/xradar/__init__.py +++ b/xradar/__init__.py @@ -20,7 +20,6 @@ # import subpackages from . import accessors # noqa from . import georeference # noqa -from . import interpolate # noqa from . import io # noqa from . import model # noqa from . import util # noqa From 5c7bd7e61ce6ac4bc9bb61251429850275469f6e Mon Sep 17 00:00:00 2001 From: Max Grover Date: Wed, 31 Jan 2024 07:36:27 -0500 Subject: [PATCH 48/48] Update xradar/io/backends/common.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Kai Mühlbauer --- xradar/io/backends/common.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xradar/io/backends/common.py b/xradar/io/backends/common.py index 32a36ea..f094b11 100644 --- a/xradar/io/backends/common.py +++ b/xradar/io/backends/common.py @@ -302,7 +302,7 @@ class LazyLoadDict(MutableMapping): The comparison methods, __cmp__, __ge__, __gt__, __le__, __lt__, __ne__, nor the view methods, viewitems, viewkeys, viewvalues, are implemented. - Neither is the the fromkeys method. + Neither is the fromkeys method. This is from Py-ART.