diff --git a/pyart/util/__init__.py b/pyart/util/__init__.py index db1ea66831..a52ffa3928 100644 --- a/pyart/util/__init__.py +++ b/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('_')] diff --git a/pyart/util/columnsect.py b/pyart/util/columnsect.py new file mode 100755 index 0000000000..4662f7a4e5 --- /dev/null +++ b/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") diff --git a/pyart/util/tests/test_columnsect.py b/pyart/util/tests/test_columnsect.py new file mode 100644 index 0000000000..ba80f2e1ec --- /dev/null +++ b/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)