Plotting SAR Fixes
The radius of maximum wind (RMW) is a key strutural parameter of tropical cyclones required to understand the evolution of winds, rainfall, and storm surge. Observations of the RMW are sparse over the ocean, however synthetic aperature radar (SAR) observations have fine enough resolution to observe the RMW in tropical cyclones. The SAR is a c-band radar on polar orbitting satellites, so the RMW is not continously sampled. When an overpass does occur, a fix of the tropical cyclone is added to the f-deck for use by forecasters to estimate the wind structure in real-time. This code below will read in the fix files and save them to numpy arrays.
Code
import copy
import numpy as np
import datetime as dt
import subprocess
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
def get_fdeck_grep(OBSTYPE,STORMPATH,FIXTYPE=''):
'''
Function to get a dataset of f-decks as a ND-numpy array
The formatting of the f-decks can be found here:
https://www.nrlmry.navy.mil/atcf_web/docs/database/new/newfdeck.txt
Input:
OBSTYPE : The observation type to get using grep format for multiple types.
i.e. '"AIRC\|RDRD"' for aircraft and radar or 'ASCT' for just ASCAT
STORMPATH : Path(s) and filename(s) for all files to grep using wildcard for multiple storms
i.e. '../dev/fdeck/fal*2019.dat''
FIXTYPE : Type of fixes to include, such as 'C' for only center fixes. Flagged fixes will not be
output
Output:
O_DATE : Datetime of observation as a string in %Y%m%D%H%M
O_STID : Stormid in standard ATCF format
O_FXFORMAT : Fix format
O_FIXTTYPE : Fix type such as AIRC for aircraft,
O_CENTERI : What was sampled, 'C' for center fix, 'R' for wind radii fix
O_FLAGIND : Flagged indicator for suspect data, same as CENTERI but with 'F' for all
O_FLAT : Latitude
O_FLON : Longitude
O_HEIGHT : Height of the observation in meters
O_POSITCONF : Position confidence, 1=good, 3=poor
O_VMAX : Maximum wind speed in kts
O_VMAXCONF : Confidence of VMAX, 1=good, 3=poor
O_MSLP : Central Pressure in mb
O_MSLPCONF : Confidence in central pressure, 1=good, 3=poor
O_MSLPDERIV : Pressure derivation such as 'DVRK' or 'MEAS' for measured
O_WRADII : Wind radii ID such as 34, 50, or 64
O_WRNE : Radius in NE quadrant in n mi
O_WRSE : Radius in SE quadrant in n mi
O_WRSW : Radius in SW quadrant in n mi
O_WRNW : Radius in NW quadrant in n mi
O_WRNEMODEL : Radii mode, E=Edge of Pass, C=cut off by land, B=both
O_WRSEMODEL : Radii mode, E=Edge of Pass, C=cut off by land, B=both
O_WRSWMODEL : Radii mode, E=Edge of Pass, C=cut off by land, B=both
O_WRNWMODEL : Radii mode, E=Edge of Pass, C=cut off by land, B=both
O_WRCONF : Wind radii confidence, 1=good, 3=poor
O_RMW : Radius of maximum winds in n mi
O_RADEYE : Diameter of eye in n mi
A_RMW : Radius of maximum winds in n mi derived from Aircraft
Example:
DATE, STID, FXFORMAT, FIXTTYPE, CENTERI, FLAGIND, FLAT, FLON, HEIGHT, POSITCONF, VMAX, VMAXCONF, MSLP, \
MSLPCONF, MSLPDERIV, WRADII, WRNE, WRSE, WRSW, WRNW, WRNEMODEL,WRSEMODEL, WRSWMODEL, WRNWMODEL, WRCONF, RMW, \
RADEYE, ARMW = get_fdeck_grep('"ASCT\|AMSU\|RDRD"' ,'../dev/fdeck/fal092021.dat',FIXTYPE='C')
'''
# First read in the b-deck file
output = subprocess.check_output("grep -h "+OBSTYPE+" "+STORMPATH, shell=True).decode().split('\n')
# Since not all lines have the same number of variables loop through and only keep the first nmax
new_data = []
nmax = 56 # max number of collumns to read in f-deck
nout = 28 # number of output variables
for i in range(0,len(output)):
data = output[i].split(',') # Split the lines by commas
if (len(data)>2): #make sure we only use lines with nmax values
# Since older files don't have RMW in the operational best-tracks we need to add data to
# to make sure the file sizes are correct
if len(data)<nmax:
nzeros = nmax-len(data) #number of values to pad the list with
data = data+[' ' for n in range(0,nzeros)] # Pad the data to get to nmax
new_data.append(data[:nmax]) #only take the 20 values
new_data = np.array(new_data).astype(str) #convert back to numpy array
# If no data is found then return missing values
if len(new_data)<1:
return [-999]*nout
# Get the number of lines
nrows = len(new_data[:,2])
#Using the unique number of rows, initialize the arrays we will use
O_DATE = np.empty((nrows),dtype='|U12')
O_STID = np.empty((nrows),dtype='|U8')
O_FXFORMAT = np.zeros((nrows))
O_FIXTTYPE = np.empty((nrows),dtype='|U4')
O_CENTERI = np.empty((nrows),dtype='|U10')
O_FLAGIND = np.empty((nrows),dtype='|U1')
O_VMAX = np.zeros((nrows))
O_VMAXCONF = np.zeros((nrows))
O_MSLP = np.zeros((nrows))
O_MSLPCONF = np.zeros((nrows))
O_FLON = np.zeros((nrows))
O_FLAT = np.zeros((nrows))
O_HEIGHT = np.zeros((nrows))
O_POSITCONF = np.zeros((nrows))
O_MSLPDERIV = np.empty((nrows),dtype='|U4')
O_WRCODE = np.zeros((nrows))
O_WRADII = np.zeros((nrows))
O_WRCONF = np.zeros((nrows))
O_WRNE = np.zeros((nrows))
O_WRSE = np.zeros((nrows))
O_WRSW = np.zeros((nrows))
O_WRNW = np.zeros((nrows))
O_WRNEMODEL = np.empty((nrows),dtype='|U1')
O_WRSEMODEL = np.empty((nrows),dtype='|U1')
O_WRSWMODEL = np.empty((nrows),dtype='|U1')
O_WRNWMODEL = np.empty((nrows),dtype='|U1')
O_RADEYE = np.zeros((nrows))
O_RMW = np.zeros((nrows))
A_RMW = np.zeros((nrows))
# Now lets loop through and save the data to the specified variables
i_dat = 0
for i_dat in range(0,nrows): # Loop through the unique datetimes
line = new_data[i_dat].astype(str) #Get the data for the specific date time
O_DATE[i_dat] = str(int(line[2]))
# Since each date could have multiple rows for each wind radii, we need to loop through them
O_STID[i_dat] = str(line[0]).strip()+str(line[1]).strip() + str(int(line[2])).strip()[:4]
O_FXFORMAT[i_dat] = int(line[3])
O_FIXTTYPE[i_dat] = line[4].strip()
O_CENTERI[i_dat] = line[5].strip()
O_FLAGIND[i_dat] = line[6].strip()
# Convert latitudes to floats
if line[7][-1:]=='N':
O_FLAT[i_dat] = float(line[7][:-1])/100.
elif line[7][-1:]=='S':
O_FLAT[i_dat] = -float(line[7][:-1])/100.
# Convert longitudes to floats
if line[8][-1:]=='E':
O_FLON[i_dat]= float(line[8][:-1])/100.
elif line[8][-1:]=='W':
O_FLON[i_dat]= -float(line[8][:-1])/100.
O_HEIGHT[i_dat] = int(line[9]) if (line[9].strip() or None) else np.nan
O_POSITCONF[i_dat] = int(line[10]) if (line[10].strip() or None) else np.nan
O_VMAX[i_dat] = int(line[11]) if (line[11].strip() or None) else np.nan
O_VMAXCONF[i_dat] = int(line[12]) if (line[12].strip() or None) else np.nan
O_MSLP[i_dat] = int(line[13]) if (line[13].strip() or None) else np.nan
O_MSLPCONF[i_dat] = int(line[14]) if (line[14].strip() or None) else np.nan
O_MSLPDERIV[i_dat] = line[15].strip()
# Check if there are wind radii information
if (line[17].strip()=='NEQ'):
O_WRADII[i_dat] = int(line[16])
if line[18]=='*****':
O_WRNE[i_dat] = np.nan
else:
O_WRNE[i_dat] = int(line[18])
if line[19]=='*****':
O_WRSE[i_dat] = np.nan
else:
O_WRSE[i_dat] = int(line[19])
if line[20]=='*****':
O_WRSW[i_dat] = np.nan
else:
O_WRSW[i_dat] = int(line[20])
if line[21]=='*****':
O_WRNW[i_dat] = np.nan
else:
O_WRNW[i_dat] = int(line[21])
O_WRNEMODEL[i_dat] = line[22].strip()
O_WRSEMODEL[i_dat] = line[23].strip()
O_WRSWMODEL[i_dat] = line[24].strip()
O_WRNWMODEL[i_dat] = line[25].strip()
O_WRCONF[i_dat] = int(line[26]) if (line[26].strip() or None) else np.nan
else:
# If there are no wind radii information then set to missing values
O_WRADII[i_dat] = np.nan
O_WRNE[i_dat] = np.nan
O_WRSE[i_dat] = np.nan
O_WRSW[i_dat] = np.nan
O_WRNW[i_dat] = np.nan
O_WRNEMODEL[i_dat] = ' '
O_WRSEMODEL[i_dat] = ' '
O_WRSWMODEL[i_dat] = ' '
O_WRNWMODEL[i_dat] = ' '
O_WRCONF[i_dat] = np.nan
O_RMW[i_dat] = int(line[27]) if (line[27].strip() or None) else np.nan
O_RADEYE[i_dat] = int(line[28]) if (line[28].strip() or None) else np.nan
# Only for aircraft
if "AIRC" in line[4].strip():
A_RMW[i_dat] = int(line[37]) if (line[37].strip() or None) else np.nan
else:
A_RMW[i_dat] = np.nan
# Check if we need to apply a fix constraint
# If a contraint is applied then we also make sure that that constraint is not flagged as suspect
if FIXTYPE!='':
cond_fx = [i for (i,x,y) in zip(range(0,len(O_CENTERI)),O_CENTERI,O_FLAGIND) if (((FIXTYPE or 'F') not in y)) ]
else:
cond_fx = [i for i in range(0,len(O_CENTERI))]
return O_DATE[cond_fx], O_STID[cond_fx], O_FXFORMAT[cond_fx], O_FIXTTYPE[cond_fx], O_CENTERI[cond_fx], O_FLAGIND[cond_fx], O_FLAT[cond_fx], O_FLON[cond_fx], O_HEIGHT[cond_fx], O_POSITCONF[cond_fx], \
O_VMAX[cond_fx], O_VMAXCONF[cond_fx], O_MSLP[cond_fx], O_MSLPCONF[cond_fx], O_MSLPDERIV[cond_fx], O_WRADII[cond_fx], O_WRNE[cond_fx], O_WRSE[cond_fx], O_WRSW[cond_fx], O_WRNW[cond_fx], O_WRNEMODEL[cond_fx],\
O_WRSEMODEL[cond_fx], O_WRSWMODEL[cond_fx], O_WRNWMODEL[cond_fx], O_WRCONF[cond_fx], O_RMW[cond_fx], O_RADEYE[cond_fx],A_RMW[cond_fx]
if __name__=='__main__':
PATH_TO_FDECK_FILES = 'FIX/f*.dat'
# Read in the data from the function specifyin the path and the observation type, here we want SAR
DATE, STID, FXFORMAT, FIXTTYPE, CENTERI, FLAGIND, FLAT, FLON, HEIGHT, POSITCONF, VMAX, VMAXCONF, MSLP, \
MSLPCONF, MSLPDERIV, WRADII, WRNE, WRSE, WRSW, WRNW, WRNEMODEL,WRSEMODEL, WRSWMODEL, WRNWMODEL, WRCONF, \
RMW, RADEYE, ARMW = get_fdeck_grep('"SAR"' ,PATH_TO_FDECK_FILES)
#--------------------------------------------------------------
# Plot RMW v. Latitude
fig = plt.figure(figsize=(10,6))
ax = plt.subplot(111)
cmap1 = copy.copy(plt.cm.rainbow)
XVAR = RMW[WRADII==34] # Don't include duplicates from R50 and R64
YVAR = FLAT[WRADII==34] # Don't include duplicates from R50 and R64
im = ax.hexbin(XVAR,YVAR,gridsize=20,alpha=.9,linewidths=0,vmin=1,cmap=cmap1)
im.cmap.set_under('w')
ax.scatter(XVAR,YVAR,c='k',s=10)
# Limit the boundaries
ax.set_xlim(0,225)
ax.set_ylim(-55,55)
# Add labels
ax.set_title('SAR Observations (2016-2023)',fontsize=18)
ax.set_xlabel('RMW (n mi)', fontsize=14)
ax.set_ylabel('Latitude (degrees)', fontsize=14)
# Increase size of xticks
ax.tick_params(axis='y',labelsize=12)
ax.tick_params(axis='x',labelsize=12)
# Remove the top and right lines
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
cb = plt.colorbar(im,pad=.02,shrink=.7)
plt.savefig('SAR_RMWvLat_Allbasins.png',bbox_inches='tight')
plt.show()
The code above reads in the fix dataset and can be used to plot SAR data. To illustrate, we can plot the distribution of RMW vs. latitude.
A figure showing latitude vs. RMW for all SAR fixes.
Enjoy Reading This Article?
Here are some more articles you might like to read next: