Merge pull request #38 from wx4stg/lmainterceptrhi

deeplycloudy · web-flow · commit e4fcc0f1950f · 2024-03-25T10:25:11.000-05:00
Add LMA Intercept RHI functionality originally by @cssjessica; ported by @wx4stg from TRACER-PAWS-NEXRAD-LMA
diff --git a/pyxlma/lmalib/lma_intercept_rhi.py b/pyxlma/lmalib/lma_intercept_rhi.py
@@ -0,0 +1,247 @@
+from pyxlma.coords import RadarCoordinateSystem, TangentPlaneCartesianSystem, GeographicSystem
+import numpy as np  
+
+
+def rcs_to_tps(radar_latitude, radar_longitude, radar_altitude, radar_azimuth):
+    """
+    Find the unit vector coordinates (east, north, up) of the plane of a radar RHI scan.
+    
+    Creates a azimuth, elevation, range and tangent plane cartesian system at the radar's latitude and longitude,
+    and converts the RHI azimuth direction to the tangent plane coordinate system.
+
+    Parameters
+    ----------
+    radar_latitude : `float`
+        Latitude of the radar in degrees.
+    radar_longitude : `float`
+        Longitude of the radar in degrees.
+    radar_altitude : `float`
+        Altitude of the radar in meters.
+    radar_azimuth : `float`
+        Azimuth of the RHI scan in degrees.
+
+    Returns
+    ----------
+    X : `numpy.ndarray`
+        A 1x2 array representing the start and end points eastward component of the RHI scan.
+    Y : `numpy.ndarray`
+        A 1x2 array representing the start and end points northward component of the RHI scan.
+    Z : `numpy.ndarray`
+        A 1x2 array representing the start and end points upward component of the RHI scan.
+    """
+
+    # Coordinates Systems
+    rcs = RadarCoordinateSystem(radar_latitude, radar_longitude, radar_altitude)
+    tps = TangentPlaneCartesianSystem(radar_latitude, radar_longitude, radar_altitude)
+
+    # - Elevations, azimuth, range
+    r = np.array([0, 1])
+
+    els = np.array([0])
+    els = np.tensordot(els, np.ones_like(r), axes=0)
+
+    azi = np.array([radar_azimuth])
+    az = np.tensordot(azi, np.ones_like(r), axes=0)
+
+    a, b, c = rcs.toECEF(r,az,els)
+    abc = np.vstack((a,b,c))
+    # ECEF to TPS
+    n = tps.toLocal(abc)
+    X = n[0,:]
+    Y = n[1,:]
+    Z = n[2,:]
+
+    X = np.reshape(X,  (1, 2))
+    Y = np.reshape(Y,  (1, 2))
+    Z = np.reshape(Z,  (1, 2))
+    
+    return X, Y, Z
+
+
+def geo_to_tps(lma_file, tps_latitude, tps_longitude, tps_altitude):
+    """
+    Convert the latitude, longitude, and altitude of LMA VHF sources to x/y/z distances (in meters) from an arbitrary latitude, longitude, and altitude.
+
+    Creates a tangent plane cartesian system at the latitude, longitude, and altitude provided, and converts the LMA VHF sources to the tangent plane coordinate system.
+
+    Parameters
+    ----------
+    lma_file : `xarray.Dataset`
+        A pyxlma dataset containing latitude, longitude, and altitude of LMA VHF sources.
+    tps_latitude : `float`
+        Latitude of the tangent plane in degrees.
+    tps_longitude : `float`
+        Longitude of the tangent plane in degrees.
+    tps_altitude : `float`
+        Altitude of the tangent plane in meters.
+    
+    Returns
+    ----------
+    Xlma : `numpy.ndarray`
+        A 1xN array representing the eastward distance (in meters) of the tangent plane center to the LMA VHF sources.
+    Ylma : `numpy.ndarray`
+        A 1xN array representing the northward distance (in meters) of the tangent plane center to the LMA VHF sources.
+    Zlma : `numpy.ndarray`
+        A 1xN array representing the upward distance (in meters) of the tangent plane center to the LMA VHF sources.
+    """
+    # GeographicSystem GEO - Lat, lon, alt
+    geo = GeographicSystem()
+    # Tangent Plane Cartesian System TPS - 
+    tps = TangentPlaneCartesianSystem(tps_latitude, tps_longitude, tps_altitude)
+
+    # GEO to TPS
+
+    d, e, h = geo.toECEF(lma_file.event_longitude.data, lma_file.event_latitude.data, lma_file.event_altitude.data)
+    
+
+    deh = np.vstack((d,e,h))
+    m = tps.toLocal(deh)
+    
+    Xlma = m[0]
+    Ylma = m[1]
+    Zlma = m[2]
+    
+    return Xlma,Ylma,Zlma
+
+
+def ortho_proj_lma(lma_file, radar_latitude, radar_longitude, radar_altitude, radar_azimuth):
+    """
+    Convert the latitude, longitude, and altitude of LMA VHF sources to distance along, distance from, and height above the ground of a radar RHI scan.
+    
+    Creates tangent plane at the radar's location and converts the LMA VHF sources to the tangent plane coordinate system, then rotates the tangent plane in the direction of the scan azimuth.
+
+    
+    Parameters
+    ----------
+    lma_file : `xarray.Dataset`
+        A pyxlma dataset containing latitude, longitude, and altitude of N number of LMA VHF sources.
+    radar_latitude : `float`
+        Latitude of the radar in degrees.
+    radar_longitude : `float`
+        Longitude of the radar in degrees.
+    radar_altitude : `float`
+        Altitude of the radar in meters.
+    radar_azimuth : `float`
+        Azimuth of the RHI scan in degrees.
+
+    Returns
+    ----------
+    lma_file_loc : `numpy.ndarray`
+        A Nx3 array representing the distance along, distance from, and height above the ground (in m) of the LMA VHF sources.
+
+    """
+    Xlma,Ylma,Zlma = geo_to_tps(lma_file, radar_latitude, radar_longitude, radar_altitude)
+    X, Y, Z = rcs_to_tps(radar_latitude, radar_longitude, radar_altitude, radar_azimuth)
+
+    lon_ini1 = X[0,0]
+    lat_ini1 = Y[0,0]
+    lon_fin1 = X[0,-1]
+    lat_fin1 = Y[0,-1]
+
+
+    dlon1 = lon_fin1 - lon_ini1 # dx
+    dlat1 = lat_fin1 - lat_ini1 # dy
+    ds1 = np.array((dlon1,dlat1))
+    norm_ds1 = np.linalg.norm(ds1)
+    cross_ds1 = np.tensordot(ds1, ds1, (0,0))
+
+    # LMA
+    lma_file_n = np.column_stack((Xlma, Ylma))
+
+    lma_file_loc_par = np.zeros(shape=(len(lma_file_n), 2))
+    lma_file_loc_perp = np.zeros(shape=(len(lma_file_n), 2))
+    lma_file_loc  = np.zeros(shape=(len(lma_file_n), 3))
+
+    #
+    # ##################################
+    #
+    #   (Xlma[i],Ylma[i]).ds1   .ds1
+    #   ----------------------
+    #        ds1 . ds1
+    #
+    # ##################################
+    #
+    lma_file_loc_tensor_x = np.tensordot(ds1,lma_file_n,(0,1))
+    lma_file_loc_par = np.tensordot((lma_file_loc_tensor_x / cross_ds1 ),ds1,0)
+
+    #
+    # #######################################################################
+    #
+    #     (Xlma[i],Ylma[i])     _     (Xlma[i],Ylma[i]).ds1   .ds1
+    #                                ----------------------
+    #                                       ds1 . ds1
+    #
+    ##########################################################################
+    #
+    lma_file_loc_perp = lma_file_n - lma_file_loc_par
+
+    #
+    lma_file_loc[:,0] = np.sqrt(lma_file_loc_par[:,0]**2 + lma_file_loc_par[:,1]**2)
+    if radar_azimuth <= 90 or radar_azimuth >= 270:
+        lma_file_loc[:,0][lma_file_loc_par[:,1] < 0] = -lma_file_loc[:,0][lma_file_loc_par[:,1] < 0]
+    elif radar_azimuth > 180 and radar_azimuth < 270:
+        lma_file_loc[:,0][lma_file_loc_par[:,0] > 0] = -lma_file_loc[:,0][lma_file_loc_par[:,0] > 0]
+    else:
+        lma_file_loc[:,0][lma_file_loc_par[:,0] < 0] = -lma_file_loc[:,0][lma_file_loc_par[:,0] < 0]
+    lma_file_loc[:,1] = np.sqrt(lma_file_loc_perp[:,0]**2 + lma_file_loc_perp[:,1]**2)
+    lma_file_loc[:,2] = Zlma
+
+    return lma_file_loc
+
+
+def find_points_near_rhi(lma_file, radar_latitude, radar_longitude, radar_altitude, radar_azimuth, radar_scan_time, distance_threshold=1000, time_threshold=30):
+    """
+    Find the LMA VHF sources near a radar RHI scan.
+    
+    Creates tangent plane at the radar's location and converts the LMA VHF sources to the tangent plane coordinate system, then rotates the tangent plane in the direction of the scan azimuth.
+    Filters RHI scan points based on distance and time thresholds.
+
+    
+    Parameters
+    ----------
+    lma_file : `xarray.Dataset`
+        A pyxlma dataset containing latitude, longitude, and altitude, and event_id of N number of LMA VHF sources.
+    radar_latitude : `float`
+        Latitude of the radar in degrees.
+    radar_longitude : `float`
+        Longitude of the radar in degrees.
+    radar_altitude : `float`
+        Altitude of the radar in meters.
+    radar_azimuth : `float`
+        Azimuth of the RHI scan in degrees.
+    radar_scan_time : `datetime.datetime` or `numpy.datetime64` or `pandas.Timestamp`
+        Time of the RHI scan.
+    distance_threshold : `float`
+        Maximum distance from the radar to the LMA VHF sources in meters. Default is 1000.
+    time_threshold : `float`
+        Number of seconds before and after the RHI scan time to include LMA VHF sources. Default is 30. (total length: 1 minute)
+
+
+    Returns
+    ----------
+    lma_range : `numpy.ndarray`
+        A 1D array representing the distance along the tangent plane in the direction of the RHI scan.
+    lma_dist : `numpy.ndarray`
+        A 1D array representing the distance from the radar RHI scan plane to each filtered LMA point.
+    lma_alt : `numpy.ndarray`
+        A 1D array representing the height above the tangent plane centered at radar level of each filtered LMA point.
+    point_mask : `numpy.ndarray`
+        A 1D array of booleans representing the VHF points that were included in the return.
+    """
+
+    radar_azimuth = radar_azimuth % 360
+
+    radar_scan_time = np.array([radar_scan_time]).astype('datetime64[s]')
+
+    projected_lma = ortho_proj_lma(lma_file, radar_latitude, radar_longitude, radar_altitude, radar_azimuth)
+    lma_range = projected_lma[:,0]
+    lma_dist = projected_lma[:,1]
+    lma_alt = projected_lma[:,2]
+
+    lma_times = lma_file.event_time.data.astype('datetime64[s]')
+    point_mask = (lma_dist < distance_threshold) & (np.abs(lma_times - radar_scan_time).astype(float) < time_threshold)
+    lma_range = lma_range[point_mask]
+    lma_dist = lma_dist[point_mask]
+    lma_alt = lma_alt[point_mask]
+
+    return lma_range, lma_dist, lma_alt, point_mask
diff --git a/tests/test_intercept_rhi.py b/tests/test_intercept_rhi.py
@@ -0,0 +1,69 @@
+from pyxlma.lmalib import lma_intercept_rhi
+from pyxlma.lmalib.io import read as lma_read
+from datetime import timedelta
+from os import listdir
+import xarray as xr
+import numpy as np
+
+
+def test_intercept_rhi():
+    """Test functionality of pyxlma.plot.lma_intercept_rhi"""
+    files_to_read = listdir('examples/data/')
+    files_to_read = ['examples/data/'+file for file in files_to_read]
+    dataset, start_date = lma_read.dataset(files_to_read)
+    lma_radar_range, lma_plane_distance, lma_arl, event_mask = lma_intercept_rhi.find_points_near_rhi(dataset, 33.56, -101.81, 1600.0, 225, start_date+timedelta(seconds=30), distance_threshold=500, time_threshold=15)
+    lma_indices = np.flatnonzero(event_mask)
+    true_radar_range = np.array([-24314.95714263, -24255.42134556, -23947.00309296, -24201.07999466,
+       -24216.22494555, -24209.0779143 , -24252.54770026, -24196.90100367,
+       -24261.77984816, -24002.07311165, -23977.02349918, -23922.33835871,
+       -24061.79138268, -24063.22463418, -24095.4865483 , -23629.07408446,
+       -23711.71702149, -23815.42007446, -23801.22424582, -23811.93829348,
+       -23880.1053328 , -23787.41575292, -23714.05514112, -23753.40801102,
+       -23697.50367439, -23808.79529204, -23702.79936852, -23733.26491339,
+       -23654.19223296, -23639.08173578, -23684.40998865, -22516.2312817 ,
+       -24956.31925763, -25087.31684706, -25278.91196761, -25284.01320286,
+       -25287.53041298, -25404.67473694, -25398.64205462, -25328.81651605,
+       -25326.06791898, -25236.99678994, -25662.94030104, -25703.02511737,
+       -25733.48865903, -25548.97709685, -25840.7709326 , -26074.465209  ,
+       -26135.94920909, -25718.73533927, -26070.13443815, -25742.36675196,
+       -26088.72774852, -26150.28911721, -22628.68886763])
+    true_plane_distance = np.array([476.82409578, 471.16428085, 324.56896919, 464.99685801,
+       443.51557595, 420.92058199, 399.02641209, 300.20108514,
+       330.2928808 , 284.48486544, 277.49727128, 228.24504655,
+       101.86785646,  46.09758876, 159.49307647, 422.34413364,
+       462.82523707, 479.04669236, 450.06522663, 495.67782222,
+       458.62470004, 397.71230587, 322.04159886, 358.54206568,
+       263.69453322, 215.27970084, 234.5031954 , 274.60201124,
+       184.85941046, 153.78011609, 203.25043091, 141.91307406,
+       291.47165215, 274.21197041, 378.27769326, 296.90514221,
+       280.02141178, 354.41770628, 327.55116845, 385.75338378,
+       177.16633758,  35.34312855, 287.23074474, 220.35513643,
+       134.09140466,  13.21239026,  82.51552308,   6.20537471,
+        46.62220121, 275.7219309 ,  67.74536475, 254.6569667 ,
+       175.35763511, 190.40871913, 100.09991281])
+    true_radar_arl = np.array([2427.10820635, 2622.80703538, 2642.20087632, 2740.2745904 ,
+       2930.63971071, 2752.31680521, 2927.29418501, 2859.50970951,
+       3625.92759368, 2942.6473481 , 2931.8120999 , 2970.00909331,
+       3111.63730722, 3062.02183736, 2859.75475096, 1394.13388557,
+       3537.06953818, 3000.85787932, 2971.55268989, 3552.81355203,
+       2836.86606787, 3226.92886031, 2460.49052586, 3200.61752493,
+       2790.06678895, 2551.85152503, 2572.98642808, 3500.51837957,
+       3063.97138001, 3031.9081397 , 3487.58191853, 6265.1033319 ,
+       2473.16902922, 2203.98364055, 2856.44738566, 2345.22649489,
+       2515.22474203, 2635.02689176, 2923.88449116, 2959.58956895,
+       2263.4423628 , 2731.91960772, 2772.38586661, 2878.88741945,
+       2746.43494287, 2584.35605103, 2341.75645168, 2613.90859694,
+       2919.30959546, 2614.69497102, 2473.25507194, 2782.82223582,
+       2461.79799166, 2377.74468257, 7332.69246886])
+    true_lma_ids = np.array([12731, 12733, 12734, 12736, 12741, 12743, 12745, 12752, 12756, 
+                             12768, 12769, 12775, 12787, 12788, 12797, 12991, 13002, 13005,
+                             13006, 13007, 13008, 13009, 13010, 13011, 13012, 13013, 13015,
+                             13016, 13023, 13029, 13030, 13201, 13388, 13390, 13395, 13396,
+                             13397, 13399, 13402, 13405, 13408, 13409, 13411, 13412, 13415,
+                             13419, 13428, 13437, 13439, 13442, 13443, 13444, 13448, 13451,
+                             13556], dtype='uint64')
+    assert np.allclose(lma_radar_range, true_radar_range, atol=1e-3)
+    assert np.allclose(lma_plane_distance, true_plane_distance, atol=1e-3)
+    assert np.allclose(lma_arl, true_radar_arl, atol=1e-3)
+    assert np.all(lma_indices == true_lma_ids)
+
