Function for storm-relative velocity

[EWolffWX]

Hi all,

I've been working to analyze some mobile radar data obtained during an event with fast storm motions and was looking for a way to subtract this out to better identify areas of rotation. From a discussion on Twitter, Dr. Tanamachi pointed my to [https://github.com/NCAR/lrose-core/blob/master/codebase/libs/Solo/src/Solo/RemoveSurface.cc](this function) within the python interface for Solo. Is there any way to integrate this directly within pyart?

Thanks in advance!!

[mgrover1]

@EWolffWX any updates here? I know you mentioned you were working on a prototype!

[EWolffWX]

Hi Max, yes, it's coming along finally! I have a few bugs to work out, but hopefully I can pass something along in the next week or two.

[wx4stg]

I'll second the desire for a function like this and throw in a "maybe I'll do this someday if I have time, etc.", unless it's already/still being worked on...

[EWolffWX]

Glad there's additional desire for this! Also just noticing my message that said "a week or two" almost a year ago (whoops, time really flew by). I actually do have a working script for this, I just need to take the time to format everything properly. If folks want to check for any extra bugs, feel free! I'll attach the draft here.

import numpy as np
import pyart 

def storm_relative_velocity(
        radar, direction, speed, field):
    """
    This function calculates storm-relative Doppler velocities.

    Created by: Edward C Wolff, 2024-06-17
    Modified by: Leanne Blind, 2024-09-12, added having a single sweep in a radar file that does not start at 0

    Parameters
    ----------
    radar: Radar
        Radar object used.
    direction: float or string
        Direction of the storm motion vector (where north equals 0 degrees). 
        Accepts a float or a string with the abbreviation of a cardinal or 
        ordinal/intercardinal direction (for example: N, SE, etc.).
    speed: string
        Speed of the storm motion vector.
        Units should be identical to those in the provided radar 
        object.
    field: string
        Velocity field to use for storm-relative calculation. A value of None 
        will use the default field name as defined in the Py-ART configuration
        file.

    Returns
    -------
    sr_data : dict
        Field dictionary containing storm-relative Doppler velocities in the
        same units as original velocities and the specified storm speed.
        Array is stored under the 'data' key.

    """
    # Parse the field parameter
    if field is None:
        field = get_field_name("velocity")
    
    # Obtain velocity data and copy the array
    sr_data = radar.fields[field]['data'].copy()

    # Specify cardinal directions that can be interpreted
    direction_dict = {
    "N": 0,
    "NE": 45,
    "E": 90,
    "SE": 135,
    "S": 180,
    "SW": 225,
    "W": 270,
    "NW": 315,
    }

    # Set the direction of the storm motion vector
    if isinstance(direction, int) or isinstance(direction, float):
        alpha = direction
    elif isinstance(direction, str):
        if direction in direction_dict.keys():
            alpha = direction_dict[direction]
        else:
            raise ValueError('Direction string must be cardinal/ordinal direction')
    else:
        raise ValueError('Direction must be an integer, float, or string')

    # Calculates the storm relative velocities 
    if len(radar.sweep_number['data']) == 1: # If the radar file only contains one sweep and not multiple (e.g. some research radars)
        sweep = 0
        start, end = radar.get_start_end(sweep)
        angle_array = radar.get_azimuth(sweep=sweep)
        ray_array = np.arange(start, end, 1)
        for count, ray in enumerate(ray_array):
            correction = speed*np.cos(np.deg2rad(alpha-angle_array[count]))
            sr_data[ray] = radar.fields[field]['data'][ray]-correction  
    else: # If the radar file contains several sweeps, one volume scan (e.g. NEXRAD)
        for sweep in radar.sweep_number['data']:
            start, end = radar.get_start_end(sweep)
            angle_array = radar.get_azimuth(sweep=sweep)
            ray_array = np.arange(start, end+1, 1)
            for count, ray in enumerate(ray_array):
                correction = speed*np.cos(np.deg2rad(alpha-angle_array[count]))
                sr_data[ray] = radar.fields[field]['data'][ray]-correction
                
    return sr_data

[zssherman]

@EWolffWX Thanks for sharing the code! Would you be interested in submitting a pull request into Py-ART with it? The documentation looks great!

[wx4stg]

haha, that's all good, I have soooo many "I'll finish that next week" projects, a few of which are over 5 years old now.

with that in mind, "I'll take a look at this next week!" 😄

also, hi Leanne! @leanneblind

[wx4stg]

I'd suggest maybe allowing u/v in addition to speed/dir as that seems to fit better within the "metpy ecosystem", but I'm not sure how best to implement that.

[EWolffWX]

@zssherman I'll definitely submit it! I'm not too well-versed in contributing, so I wasn't sure if there were any formatting things I had forgotten to account for

[mgrover1]

@zssherman I'll definitely submit it! I'm not too well-versed in contributing, so I wasn't sure if there were any formatting things I had forgotten to account for

We can help with that part as well :) go ahead and submit a pull request, and we can walk you through the process!

[scollis]

100% new contributors are more than welcome

[EWolffWX]

@wx4stg That's a good idea. I can work on adding that functionality too. Should be a simple line or two that converts to the speed and direction style the rest of the script uses... (hopefully)

[EWolffWX]

@mgrover1 and @scollis Great, thank you all! I'll get on that this afternoon