ADD: Utility to subset a radar volume for the column above a location (…

…#1113)

* ADD: Utility to subset a radar volume for the column above a location
given the lat/lon of the location. Return is a dictionary containing all radar object fields for the column above the location, as well as, x,y,z gates for the column. Tested with ppi (nexrad) and rhi (csapr) scans

* MOD: first corrections to initial PR for columnsect.py
Still working on transfering output to an xarry object

* ADD: Unit Tests for sphere_distance, for_azimuth, and get_field_location functions within
columnsect.py utility

* MOD: Added the requested changes from reviewers. Main modifications include changing the output of
get_field_location form a dictionary to xarray Datatset.

pyart/util/__init__.py

@@ -19,5 +19,7 @@
 from .simulated_vel import simulated_vel_from_profile
 from .sigmath import texture_along_ray, rolling_window
 from .sigmath import texture, angular_texture_2d
+from .columnsect import get_field_location, get_column_rays
+from .columnsect import for_azimuth, sphere_distance
 
 __all__ = [s for s in dir() if not s.startswith('_')]

pyart/util/columnsect.py

@@ -0,0 +1,322 @@
+"""
+Function for extracting the radar column above a target
+given position in latitude, longitude
+
+"""
+
+import numpy as np
+import xarray as xr
+
+from ..core.transforms import antenna_vectors_to_cartesian
+
+
+def sphere_distance(radar_latitude, target_latitude, radar_longitude,
+                    target_longitude):
+    """
+    Calculated of the great circle distance between radar and target
+
+    Assumptions
+    -----------
+    Radius of the Earth = 6371 km / 6371000 meters
+    Distance is calculated for a smooth sphere
+    Radar and Target are at the same altitude (need to check)
+
+    Parameter
+    ---------
+    radar_latitude : float, [degrees]
+        latitude of the radar in degrees
+    target_latitude : float, [degrees]
+        latitude of the target in degrees
+    radar_longitude : float, [degrees]
+        longitude of the radar in degrees
+    target_longitude : float, [degrees]
+        longitude of the target in degress
+
+    Returns
+    -------
+    distance : float, [meters]
+        Great-Circle Distance between radar and target in meters
+    """
+    # check if latitude, longitudes are valid
+    check_latitude(radar_latitude)
+    check_latitude(target_latitude)
+    check_longitude(radar_longitude)
+    check_longitude(target_longitude)
+
+    # convert latitudeitude/longitudegitudes to radians
+    radar_latitude = radar_latitude * (np.pi/180.)
+    target_latitude = target_latitude * (np.pi/180.)
+    radar_longitude = radar_longitude * (np.pi/180.)
+    target_longitude = target_longitude * (np.pi/180.)
+
+    # difference in latitudeitude
+    d_latitude = (target_latitude - radar_latitude)
+    # difference in longitudegitude
+    d_longitude = (target_longitude - radar_longitude)
+
+    # Haversine formula
+    numerator = ((np.sin(d_latitude/2.0)**2.0)
+                 + np.cos(radar_latitude) * np.cos(target_latitude)
+                 * (np.sin(d_longitude/2.0)**2.0))
+    distance = 2 * 6371000 * np.arcsin(np.sqrt(numerator))
+
+    # return the output
+    return distance
+
+
+def for_azimuth(radar_latitude, target_latitude, radar_longitude,
+                target_longitude):
+    """
+    Calculation of inital bearing alongitudeg a great-circle arc
+    Known as Forward Azimuth Angle.
+
+    Assumptions
+    -----------
+    Radius of the Earth = 6371 km / 6371000 meters
+    Distance is calculatitudeed for a smooth sphere
+    Radar and Target are at the same altitude (need to check)
+
+    Parameters
+    ----------
+    radar_latitude : float, [degrees]
+        latitude of the radar in degrees
+    target_latitude : float, [degrees]
+        latitude of the target in degrees
+    radar_longitude : float, [degrees]
+        longitude of the radar in degrees
+    target_longitude : float, [degrees]
+        longitude of the target in degress
+
+    Returns
+    -------
+    azimuth : float, [degrees]
+        azimuth angle from the radar where
+        target is located within the scan.
+        output is in degrees.
+    """
+    # check the input latitude/longitude
+    check_latitude(radar_latitude)
+    check_latitude(target_latitude)
+    check_longitude(radar_longitude)
+    check_longitude(target_longitude)
+
+    # convert latitudeitude/longitudegitudes to radians
+    radar_latitude = radar_latitude * (np.pi/180.)
+    target_latitude = target_latitude * (np.pi/180.)
+    radar_longitude = radar_longitude * (np.pi/180.)
+    target_longitude = target_longitude * (np.pi/180.)
+
+    # Differnce in longitudegitudes
+    d_longitude = target_longitude - radar_longitude
+
+    # Determine x,y coordinates for arc tangent function
+    corr_y = np.sin(d_longitude) * np.cos(target_latitude)
+    corr_x = ((np.cos(radar_latitude) * np.sin(target_latitude))
+              - (np.sin(radar_latitude) * (np.cos(target_latitude)
+                                           * np.cos(d_longitude))))
+
+    # Determine forward azimuth angle
+    azimuth = np.arctan2(corr_y, corr_x) * (180./np.pi)
+
+    # Return the output as a function of 0-360 degrees
+    if azimuth < 0:
+        azimuth += 360.
+
+    return azimuth
+
+
+def get_field_location(radar, latitude, longitude):
+
+    """
+    Given the location (in latitude, longitude) of a target, extract the
+    radar column above that point for further analysis.
+
+    Parameters
+    ----------
+    radar : pyart.core.Radar Object
+        Py-ART Radar Object from which distance to the target, along
+        with gates above the target, will be calculated.
+    latitude : float, [degrees]
+        Latitude, in degrees North, of the target.
+    longitude : float, [degrees]
+        Longitude, in degrees East, of the target.
+
+    Function Calls
+    --------------
+    sphere_distance
+    for_azimuth
+    get_column_rays
+
+    Returns
+    -------
+    column : xarray DataSet
+        Xarary Dataset containing the radar column above the target for
+        the various fields within the radar object.
+    """
+    # Make sure latitude, longitudes are valid
+    check_latitude(latitude)
+    check_longitude(longitude)
+
+    # initiate a dictionary to hold the gates above the radar.
+    zgate = []
+
+    # initiate a diciontary to hold the moment data.
+    moment = {key: [] for key in radar.fields.keys()}
+
+    # call the sphere_distance function
+    dis = sphere_distance(radar.latitude['data'][0],
+                          latitude,
+                          radar.longitude['data'][0],
+                          longitude)
+
+    # call the for_azimuth function
+    azim = for_azimuth(radar.latitude['data'][0],
+                       latitude,
+                       radar.longitude['data'][0],
+                       longitude)
+
+    # call the get_column_ray function
+    ray = get_column_rays(radar, azim)
+
+    # Determine the gates for the rays
+    (rhi_x,
+     rhi_y,
+     rhi_z) = antenna_vectors_to_cartesian(radar.range['data'],
+                                           radar.azimuth['data'][ray],
+                                           radar.elevation['data'][ray],
+                                           edges=True)
+    # Calculate distance from the x,y coordinates to target
+    rhidis = np.sqrt((rhi_x**2) + (rhi_y**2)) * np.sign(rhi_z)
+    for i in range(len(ray)):
+        tar_gate = np.argmin(abs(rhidis[i, 1:] - (dis)))
+        for key in moment:
+            if radar.fields[key]['data'][ray[i], tar_gate] is np.ma.masked:
+                moment[key].append(np.nan)
+            else:
+                moment[key].append(radar.fields[key]
+                                   ['data'][ray[i], tar_gate])
+        zgate.append(rhi_z[i, tar_gate])
+
+    # Create a blank list to hold the xarray DataArrays
+    ds_container = []
+    da_meta = ['units', 'standard_name', 'long_name', 'valid_max',
+               'valid_min']
+    # Convert the moment dictionary to xarray DataArray.
+    # Apply radar object meta data to DataArray attribute
+    for key in moment:
+        da = xr.DataArray(moment[key], coords=[zgate],
+                          name=key, dims=['height'])
+        for tag in da_meta:
+            if tag in radar.fields[key]:
+                da.attrs[tag] = radar.fields[key][tag]
+        # Append to ds container
+        ds_container.append(da.to_dataset(name=key))
+
+    # Create a xarray DataSet from the DataArrays
+    column = xr.merge(ds_container)
+    if len(radar.time['units'].split(' ')) > 3:
+        column.attrs['date'] = radar.time['units'].split(' ')[2]
+        column.attrs['time'] = radar.time['units'].split(' ')[3]
+    else:
+        column.attrs['date'] = radar.time['units'].split(' ')[2].split('T')[0]
+        column.attrs['time'] = radar.time['units'].split(' ')[2].split('T')[-1]
+    column.attrs['distance'] = (str(np.around(dis / 1000., 3)) + ' km')
+    column.attrs['azimuth'] = (str(np.around(azim, 3)) + ' degrees')
+    column.attrs['latitude'] = (str(latitude) + ' degrees')
+    column.attrs['longitude'] = (str(longitude) + ' degrees')
+
+    return column
+
+
+def get_column_rays(radar, azimuth):
+
+    """
+    Given the location (in latitude,longitude) of a target, return the rays
+    that correspond to radar column above the target.
+
+    Parameters
+    ----------
+    radar    : Radar Object
+               Py-ART Radar Object from which distance to the target, along
+               with gates above the target, will be calculated.
+    azimuth  : float,int
+               forward azimuth angle from radar to target in degrees.
+
+    Returns
+    -------
+    nrays : List
+            radar ray indices that correspond to the column above a
+            target location
+    """
+    if isinstance(azimuth, int) or isinstance(azimuth, float) is True:
+        if (azimuth <= 0) or (azimuth >= 360):
+            raise ValueError("azimuth not valid (not between 0-360 degrees)")
+    else:
+        raise TypeError("radar longitudegitude type not valid."
+                        + " Please convert input to be an int or float")
+    # define a list to hold the valid rays
+    rays = []
+    # check to see which radar scan
+    if radar.scan_type == 'rhi':
+        for i in range(radar.sweep_number['data'].shape[0]):
+            nstart = radar.sweep_start_ray_index['data'][i]
+            nstop = radar.sweep_end_ray_index['data'][i]
+            counter = 0
+            for j in range(nstart, nstop):
+                if (abs(radar.azimuth['data'][nstart + counter] -
+                        azimuth) < 1):
+                    rays.append(nstart+counter)
+                counter += 1
+    else:
+        # taken from pyart.graph.RadarDisplay.get_azimuth_rhi_data_x_y_z
+        for sweep in radar.iter_slice():
+            sweep_azi = radar.azimuth['data'][sweep]
+            nray = np.argmin(np.abs(sweep_azi - azimuth))
+            rays.append(nray + sweep.start)
+    # make sure rays were found
+    if len(rays) == 0:
+        raise ValueError("No rays were found between azimuth and target")
+
+    return rays
+
+
+def check_latitude(latitude):
+
+    """
+    Function to check if input latitude is valid for type and value
+
+    Parameters
+    ----------
+    latitude : int, float
+        Latitude of a location that should be between 90S and 90N
+
+    """
+    if (isinstance(latitude, int) or isinstance(latitude, float) or
+            isinstance(latitude, np.floating)) is True:
+        if (latitude <= -90) or (latitude >= 90):
+            raise ValueError("Latitude not between -90 and 90 degrees, need to"
+                             + "convert to values between -90 and 90")
+    else:
+        raise TypeError("Latitude type not valid, need to convert input to be"
+                        + " an int or float")
+
+
+def check_longitude(longitude):
+
+    """
+    Function to check if input latitude is valid for type and value
+
+    Parameters
+    ---------
+    longitude : int, float
+        Longitude of a location taht should be between 180W and 180E
+    """
+    if (isinstance(longitude, int) or isinstance(longitude, float) or
+            isinstance(longitude, np.floating)) is True:
+        if (longitude <= -180) or (longitude >= 180):
+            raise ValueError("Longitude not valid between -180 and 180"
+                             + " degrees, need to convert to values between"
+                             + "  -180 and 180")
+    else:
+        raise TypeError("Longitude type not valid, need to convert input to"
+                        + " be an int or float")

pyart/util/tests/test_columnsect.py

@@ -0,0 +1,40 @@
+""" Unit Tests for Py-ART's util/columnsect.py module. """
+
+import numpy as np
+import pyart
+
+# read in example file
+radar = pyart.io.read_nexrad_archive(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE)
+
+
+def test_get_azimuth():
+    # test to make sure azimuth is correct to Everett, WA
+    azimuth = pyart.util.columnsect.for_azimuth(radar.latitude['data'][0],
+                                                47.97,
+                                                radar.longitude['data'][0],
+                                                -122.20)
+    test_azi = abs(np.around(azimuth, 2) - 138.57)
+    assert (test_azi < 0.001)
+
+
+def test_sphere_distance():
+    # test to make sure sphere distance is correct to Everett, WA
+    distance = pyart.util.columnsect.sphere_distance(radar.latitude['data'][0],
+                                                     47.97,
+                                                     (radar.
+                                                      longitude['data'][0]),
+                                                     -122.20)
+    test_dis = abs(np.around(distance / 1000., 2) - 33.27)
+    assert (test_dis < 0.001)
+
+
+def test_get_field_location():
+    # test to make sure column above location is pulled correctly
+    column = pyart.util.columnsect.get_field_location(radar, 47.97, -122.20)
+    # check to make sure z-gate is pulled correctly.
+    test_height = abs(column.height[0] - 367)
+    assert (test_height < 0.001)
+
+    # check to make sure reflectivity value is minimum
+    test_z = abs(column.reflectivity[0] + 32)
+    assert (test_z < 0.001)