Source code for pyart.io.rsl

"""
Python wrapper around the RSL library.

"""

import datetime
import warnings

import numpy as np

from ..config import FileMetadata, get_fillvalue

try:
    from . import _rsl_interface

    _RSL_AVAILABLE = True
except ImportError:
    _RSL_AVAILABLE = False
from ..core.radar import Radar
from ..exceptions import MissingOptionalDependency
from ..lazydict import LazyLoadDict
from .common import make_time_unit_str


[docs]def read_rsl( filename, field_names=None, additional_metadata=None, file_field_names=False, exclude_fields=None, delay_field_loading=False, include_fields=None, radar_format=None, callid=None, skip_range_check=False, ): """ Read a file supported by RSL. Parameters ---------- filename : str or RSL_radar Name of file whose format is supported by RSL. field_names : dict, optional Dictionary mapping RSL data type names to radar field names. If a data type found in the file does not appear in this dictionary or has a value of None it will not be placed in the radar.fields dictionary. A value of None, the default, will use the mapping defined in the Py-ART configuration file. additional_metadata : dict of dicts, optional Dictionary of dictionaries to retrieve metadata from during this read. This metadata is not used during any successive file reads unless explicitly included. A value of None, the default, will not introduct any addition metadata and the file specific or default metadata as specified by the Py-ART configuration file will be used. file_field_names : bool, optional True to use the RSL data type names for the field names. If this case the field_names parameter is ignored. The field dictionary will likely only have a 'data' key, unless the fields are defined in `additional_metadata`. exclude_fields : list or None, optional List of fields to exclude from the radar object. This is applied after the `file_field_names` and `field_names` parameters. Set to None to include all fields specified by include_fields. include_fields : list or None, optional List of fields to include from the radar object. This is applied after the `file_field_names` and `field_names` parameters. Set to None to include all fields not specified by exclude_fields. delay_field_loading : bool True to delay loading of field data from the file until the 'data' key in a particular field dictionary is accessed. In this case the field attribute of the returned Radar object will contain LazyLoadDict objects not dict objects. radar_format : str or None Format of the radar file. Must be 'wsr88d' or None. callid : str or None Four letter NEXRAD radar Call ID, only used when radar_format is 'wsr88d'. skip_range_check : bool, optional True to skip check for uniform range bin location, the reported range locations will only be verified true for the first ray. False will perform the check and raise a IOError when the locations of the gates change between rays. Returns ------- radar : Radar Radar object. """ # check that TRMM RSL is available if not _RSL_AVAILABLE: raise MissingOptionalDependency( "Py-ART must be build with support for TRMM RSL to use the " "read_rsl function." ) # create metadata retrieval object filemetadata = FileMetadata( "rsl", field_names, additional_metadata, file_field_names, exclude_fields, include_fields, ) # read the file, determine common parameters fillvalue = get_fillvalue() rslfile = _rsl_interface.RslFile(filename.encode("ascii"), radar_format, callid) available_vols = rslfile.available_moments() first_volume = rslfile.get_volume(available_vols[0]) first_sweep = first_volume.get_sweep(0) first_ray = first_sweep.get_ray(0) nsweeps = first_volume.nsweeps # scan_type, sweep_mode, fixed_angle fixed_angle = filemetadata("fixed_angle") fdata = first_volume.get_sweep_fix_angles() fdata[fdata < 0.0] = 360.0 + fdata[fdata < 0.0] # all angles positive fixed_angle["data"] = fdata sweep_mode = filemetadata("sweep_mode") if rslfile.scan_mode == 0: # 0 = PPI, 1 = RHI scan_type = "ppi" sweep_mode["data"] = np.array(nsweeps * ["azimuth_surveillance"], dtype="S") else: scan_type = "rhi" sweep_mode["data"] = np.array(nsweeps * ["rhi"], dtype="S") # time time = filemetadata("time") t_start = datetime.datetime( rslfile.year, rslfile.month, rslfile.day, rslfile.hour, rslfile.minute, int(rslfile.sec), int(1e6 * (rslfile.sec % 1)), ) datetimes = [] for i in range(nsweeps): sweep = first_volume.get_sweep(i) for j in range(sweep.nrays): datetimes.append(sweep.get_ray(j).get_datetime()) t_delta = [t - t_start for t in datetimes] sec_since_start = [ t.seconds + t.days * 3600 * 24 + t.microseconds / 1.0e6 for t in t_delta ] time["data"] = np.array(sec_since_start, dtype=np.float64) time["units"] = make_time_unit_str(t_start) # range if (not skip_range_check) and (not first_volume.is_range_bins_uniform()): message = ( "Range bin locations change between rays. File cannot be read " "with with correct range locations for all rays. " "To read in data reporting the ranges from the first ray set the " "**skip_range_check** parameter to True." ) raise OSError(message) _range = filemetadata("range") gate0 = first_ray.range_bin1 gate_size = first_ray.gate_size ngates = first_ray.nbins _range["data"] = gate0 + gate_size * np.arange(ngates, dtype="float32") _range["meters_to_center_of_first_gate"] = _range["data"][0] _range["meters_between_gates"] = np.array(gate_size, dtype="float32") # fields # transfer only those which are available fields = {} for volume_num in available_vols: # if a volume number is not recognized it is skipped. if volume_num not in VOLUMENUM2RSLNAME: warnings.warn("Unknown Volume Number %d" % volume_num) continue rsl_field_name = VOLUMENUM2RSLNAME[volume_num] field_name = filemetadata.get_field_name(rsl_field_name) if field_name is None: continue # create the field dictionary field_dic = filemetadata(field_name) field_dic["_FillValue"] = fillvalue data_extractor = _RslVolumeDataExtractor(rslfile, volume_num, fillvalue) if delay_field_loading: field_dic = LazyLoadDict(field_dic) field_dic.set_lazy("data", data_extractor) else: field_dic["data"] = data_extractor() fields[field_name] = field_dic # metadata metadata = filemetadata("metadata") metadata["original_container"] = "rsl" rsl_dict = rslfile.get_radar_header() need_from_rsl_header = { "name": "instrument_name", "project": "project", "state": "state", "country": "country", } # rsl_name : radar_metadata_name for rsl_key, metadata_key in need_from_rsl_header.items(): metadata[metadata_key] = rsl_dict[rsl_key] # latitude latitude = filemetadata("latitude") lat = _dms_to_d((rsl_dict["latd"], rsl_dict["latm"], rsl_dict["lats"])) latitude["data"] = np.array([lat], dtype="float64") # longitude longitude = filemetadata("longitude") lon = _dms_to_d((rsl_dict["lond"], rsl_dict["lonm"], rsl_dict["lons"])) longitude["data"] = np.array([lon], dtype="float64") # altitude altitude = filemetadata("altitude") altitude["data"] = np.array([rsl_dict["height"]], dtype="float64") # sweep_number, sweep_mode, sweep_start_ray_index, sweep_end_ray_index sweep_number = filemetadata("sweep_number") sweep_start_ray_index = filemetadata("sweep_start_ray_index") sweep_end_ray_index = filemetadata("sweep_end_ray_index") sweep_number["data"] = np.arange(nsweeps, dtype="int32") ray_count = first_volume.get_nray_array() # array of rays in each sweep ssri = np.cumsum(np.append([0], ray_count[:-1])).astype("int32") sweep_start_ray_index["data"] = ssri sweep_end_ray_index["data"] = np.cumsum(ray_count).astype("int32") - 1 # azimuth, elevation azimuth = filemetadata("azimuth") elevation = filemetadata("elevation") _azimuth, _elevation = first_volume.get_azimuth_and_elev_array() azimuth["data"] = _azimuth elevation["data"] = _elevation # instrument_parameters prt = filemetadata("prt") prt_mode = filemetadata("prt_mode") nyquist_velocity = filemetadata("nyquist_velocity") unambiguous_range = filemetadata("unambiguous_range") beam_width_h = filemetadata("radar_beam_width_h") beam_width_v = filemetadata("radar_beam_width_v") pm_data, nv_data, pr_data, ur_data = first_volume.get_instr_params() prt["data"] = pr_data prt_mode["data"] = pm_data.astype("S") nyquist_velocity["data"] = nv_data unambiguous_range["data"] = ur_data beam_width_h["data"] = np.array([first_sweep.horz_half_bw * 2.0], dtype="float32") beam_width_v["data"] = np.array([first_sweep.vert_half_bw * 2.0], dtype="float32") instrument_parameters = { "unambiguous_range": unambiguous_range, "prt_mode": prt_mode, "prt": prt, "nyquist_velocity": nyquist_velocity, "radar_beam_width_h": beam_width_h, "radar_beam_width_v": beam_width_v, } return Radar( time, _range, fields, metadata, scan_type, latitude, longitude, altitude, sweep_number, sweep_mode, fixed_angle, sweep_start_ray_index, sweep_end_ray_index, azimuth, elevation, instrument_parameters=instrument_parameters, )
def _dms_to_d(dms): """Degrees, minutes, seconds to degrees""" return dms[0] + (dms[1] + dms[2] / 60.0) / 60.0 class _RslVolumeDataExtractor: """ Class facilitating on demand extraction of data from a RSL file. Parameters ---------- rslfile : RslFile Open RslFile object to extract data from. volume_num : int Volume number of data to be extracted. fillvalue : int Value used to fill masked values in the returned array. """ def __init__(self, rslfile, volume_num, fillvalue): """initialize the object.""" self.rslfile = rslfile self.volume_num = volume_num self.fillvalue = fillvalue def __call__(self): """Return an array containing data from the referenced volume.""" # extract the field and mask data = self.rslfile.get_volume_array(self.volume_num) data[np.where(np.isnan(data))] = self.fillvalue data[np.where(data == 131072)] = self.fillvalue return np.ma.masked_equal(data, self.fillvalue) VOLUMENUM2RSLNAME = { 0: "DZ", 1: "VR", 2: "SW", 3: "CZ", 4: "ZT", 5: "DR", 6: "LR", 7: "ZD", 8: "DM", 9: "RH", 10: "PH", 11: "XZ", 12: "CD", 13: "MZ", 14: "MD", 15: "ZE", 16: "VE", 17: "KD", 18: "TI", 19: "DX", 20: "CH", 21: "AH", 22: "CV", 23: "AV", 24: "SQ", 25: "VS", 26: "VL", 27: "VG", 28: "VT", 29: "NP", 30: "HC", 31: "VC", 32: "V2", 33: "S2", 34: "V3", 35: "S3", 36: "CR", 37: "CC", 38: "PR", 39: "SD", 40: "ZZ", 41: "RD", 42: "ET", 43: "EZ", } RSLNAME2VOLUMENUM = {v: k for k, v in VOLUMENUM2RSLNAME.items()}