azim_wind_latlon.py

# Function to convert wind from (time,lat,lon) to the tangential
#   wind as f(time,z,x,y) based on a TC track.
# 
# INPUTS:
#       track_file - file path and name for netCDF file containing track as clon,clat
#           (generated by the program run_tracking.py)
#       u - zonal wind as (time,x,y) (note, the order of x,y is arbitrary)
#       v - meridional wind " "
#       lon - array of longitude points as (x,y) [deg]
#       lat - " " latitude points
#       t0, t1 - bounding time indices of var, assuming shorter in time than clon/clat
# 
# RETURNS:
#       returns a time,x,y series masked array of identical shape to var.
# 
# James Ruppert  
# jruppert@ou.edu  
# September 2022

# from netCDF4 import Dataset
# import numpy as np

from netCDF4 import Dataset
import numpy as np
import sys

def azim_wind_latlon(track_file, u, v, lon, lat, t0, t1):

    # Settings
    # rmax = 3 # degrees

    # Input dimensions
    nt,nz,nx1,nx2 = u.shape

    # Function to account for crossing of the Intl Date Line
    def dateline_lon_shift(lon_in, reverse):
        if reverse == 0:
            lon_offset = np.zeros(lon_in.shape)
            lon_offset[np.where(lon_in < 0)] += 360
        else:
            lon_offset = np.zeros(lon_in.shape)
            lon_offset[np.where(lon_in > 180)] -= 360
        # return lon_in + lon_offset
        return lon_offset
    
    # Read TC track
    ncfile = Dataset(track_file)
    clon = ncfile.variables['clon'][:] # deg
    clat = ncfile.variables['clat'][:] # deg
    ncfile.close()

    # Check for crossing Date Line
    if (lon.min() < 0) and (lon.max() > 0):
        lon_offset = dateline_lon_shift(lon, reverse=0)
        clon_offset = dateline_lon_shift(clon, reverse=0)
    else:
        lon_offset = 0
        clon_offset = 0

    # Center grid on track
    lon4d = np.repeat((lon+lon_offset)[np.newaxis,np.newaxis,:,:], nt, axis=0)
    lat4d = np.repeat(lat[np.newaxis,np.newaxis,:,:], nt, axis=0)
    lon4d -= (clon+clon_offset)[t0:t1,np.newaxis,np.newaxis,np.newaxis]
    lat4d -= clat[t0:t1,np.newaxis,np.newaxis,np.newaxis]

    # Get radius and azimuth
    radius4d = np.sqrt( lon4d**2 + lat4d**2 )
    azim = np.arctan(lat4d/lon4d)
    azim[(lon4d < 0 )] += np.pi

    # Subtract storm motion
    deg2x = 6371e3*np.pi/180 # m / deg
    # Assume hourly time steps
    dt = 3600. # s/time step
    mpsec = deg2x/dt # m / deg / (s / time step)
    u_track = np.gradient(clon+clon_offset) * np.cos(clat*np.pi/180) * mpsec # deg / time step --> m / s
    v_track = np.gradient(clat) * mpsec # deg / time step --> m / s
    # print('U = ',u_track)
    # print('V = ',v_track)

    # Get tangential wind
    vtx = (u - u_track[t0:t1,np.newaxis,np.newaxis,np.newaxis]) * np.sin(azim)
    vty = (v - v_track[t0:t1,np.newaxis,np.newaxis,np.newaxis]) * np.cos(azim)
    # vtx = u * np.sin(azim)
    # vty = v * np.cos(azim)
    v_tan = vty - vtx

    return v_tan