Calculating rain rate from radar reflectivity of Hurricane Julia

Skills - python, hdf files, rasterio, ECMWF, weather forecasting, geoproj, numpy, pyproj, cloudsat

Introduction

This research project was carried out to see radar reflectivity values in storms, and calculate the rain rate from the radar pulses. The raw data files are in the hdf format and obtained from Cloudsat’s Data Center.

Related projects: Hurricane Karl
This work finds the the point at which the ECMWF tempearture= 0 deg C for each radar pulse and overlays that on the reflectivity plot to check to see whether the bright band occurs at the freezing level. It will also plot the rain rate for the storm track along with the corresponding radar reflectivity to study the eye of the storm

Storm Details

Name: Hurricane Julia
Year: 2010
Region: Atlantic Basin (did not make landfall)
Wiki Reference and Storm Details

The full notebook with this work is at the github project for this analysis.



1. After reading the height and reflectivity raw data from the hdf files, I plot the rain rates:

rain_rate = HDFvd_read(r_file,'rain_rate',vgroup='Data Fields')
invalid = (rain_rate == -9999.)
rain_rate[invalid] = np.nan #creates runtime warning due to nan value
hit = rain_rate < 0.
rain_rate[hit] = np.abs(rain_rate[hit])
plt.plot(rain_rate)

2. Creating a masked array of the reflectivity so that pcolormesh will plot it

hit=(radar_reflectivity == radar_missing)
radar_reflectivity=radar_reflectivity.astype(np.float)
radar_reflectivity[hit]=np.nan
zvals = radar_reflectivity/radar_scale
zvals=ma.masked_invalid(zvals)

3. Find the 3 minutes of the pass that includes the storm, then convert time to distance to get the greatcircle distance between shots

great_circle=pyproj.Geod(ellps='WGS84')
distance=[0]
start=(storm_lons[0],storm_lats[0])
for index in np.arange(1,len(storm_lons)):
    azi12,azi21,step= great_circle.inv(storm_lons[index-1],storm_lats[index-1],
                                       storm_lons[index],storm_lats[index])
    distance.append(distance[index-1] + step)
distance=np.array(distance)/meters2km

4. Make the plot for rain rate and reflectivity


%matplotlib inline

from matplotlib import cm
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import gridspec

def plot_field2(distance,height,field,fig,cmap=None,norm=None):
    """
    draw a 2 panel plot with different panel sizes.  Put the radar reflectivity in the top panel 
    with a colorbar along the bottom, and pass the second axis back to be filled in later.  Uses 
    the sharex keyword to give both plots the same x axis (distance) and the gridspec class to 
    layout the grid
    """
    
    gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1]) 
    ax1 = fig.add_subplot(gs[0])
    ax2 = fig.add_subplot(gs[1],sharex=ax1)
    if cmap is None:
        cmap=cm.inferno
    col=ax1.pcolormesh(distance,height,field,cmap=cmap, norm=the_norm)
    divider = make_axes_locatable(ax1)
    cax = divider.append_axes("bottom", size="5%", pad=0.55)
    ax1.figure.colorbar(col,extend='both',cax=cax,orientation='horizontal')
    return ax1, ax2

vmin=-30
vmax=20
the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
cmap_ref=cm.plasma
cmap_ref.set_over('w')
cmap_ref.set_under('b',alpha=0.2)
cmap_ref.set_bad('0.75') #75% grey

cloud_height_km=radar_height[0,:]/meters2km
fig = plt.figure(figsize=(15, 8)) 
ax1, ax2 = plot_field2(distance,cloud_height_km,storm_zvals.T,fig,cmap=cmap_ref,norm=the_norm)
ax1.set(ylim=[0,17],xlim=(0,1200))
ax1.set(xlabel='distance (km)',ylabel='height (km)',title='equivalent radar reflectivity in dbZe')

ax2.plot(distance,rain_rate)
ax2.set(xlabel='distance (km)',ylabel='rain rate (mm/hour)')

5. Get ECMWF temps

z_file= list(a301.data_dir.glob('*ECMWF-AUX_GRANULE*hdf'))[2]
ec_height=HDFvd_read(z_file,'EC_height')
ec_temps, temps_attributes = HDFsd_read(z_file,'Temperature')
ec_missing = temps_attributes['missing']

6. Subsetting the ECMWF data

bad_temps = (ec_temps == ec_missing)
ec_temps[bad_temps]=np.nan
ec_temps=np.ma.masked_invalid(ec_temps)
ec_temps = ec_temps - 273.15
ec_temps=ec_temps[time_hit,:]

7. ECMWF temperatures for the segment

def plot_field(distance,height,field,ax,cmap=None,norm=None):
    """
    given an axis, draw a cloudsat cross section
    """
    if cmap is None:
        cmap=cm.inferno
    col=ax.pcolormesh(distance,height,field,cmap=cmap,
                  norm=the_norm)
    ax.figure.colorbar(col,extend='both',ax=ax)
    return ax

fig, ax =plt.subplots(1,1,figsize=(15,4))
vmin=-30
vmax=30
the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
cmap_ec= cm.bwr
cmap_ec.set_over('w')
cmap_ec.set_under('b',alpha=0.2)
cmap_ec.set_bad('0.75') #75% grey

ec_height_km=ec_height/meters2km
ax=plot_field(distance,ec_height_km,ec_temps.T,ax,cmap=cmap_ec,
                  norm=the_norm)
ax.set(ylim=[0,10],xlim=(0,1200))
ax.set(xlabel='distance (km)',ylabel='height (km)',title='ECMWF temps in deg C')

8. Read in heating rate


hr_file= list(a301.data_dir.glob('*FLXHR*hdf'))[2]
lats,lons,date_times,prof_times,dem_elevation=get_geo(hr_file)
lats=lats.squeeze()
lons=lons.squeeze()
qr, qr_attrs = HDFsd_read(hr_file,'QR')
qr_height, height_attrs = HDFsd_read(hr_file,'Height')
factor = HDFvd_read(hr_file,'QR.factor',vgroup='Swath Attributes')[0][0]
missing = HDFvd_read(hr_file,'QR.missing',vgroup='Swath Attributes')[0][0]
units = HDFvd_read(hr_file,'QR.units',vgroup='Swath Attributes')[0][0]
#set_trace()
hit = (qr == missing)
qr = qr.astype(np.float64)/factor
qr[hit]=np.nan


storm_qr=qr[:,time_hit,:]
storm_height=qr_height[time_hit,:]

9. Split out the long and shortwave heating rates

shortwave_qr=storm_qr[0,:,:]
longwave_qr=storm_qr[1,:,:]

10. Make a plots for shortwave and longwave heating rates which look like so:




More about Cloudsat

A detailed overview of the Cloudsat mission is given in this Bulletin of the AMS article article on cloudsat.

From the Wikipedia Cloudsat page :

Cloudsat is a NASA Earth observation satellite, which was launched on a Delta II rocket on April 28, 2006. It uses radar to measure the altitude and properties of clouds, adding to information on the relationship between clouds and climate in order to help resolve questions about global warming.

From the Cloudsat Project Page at the Dept. of Atmospheric Science, Colorado State U:

CloudSat flies in formation in the A Train, with several other satellites: Aqua, Aura, CALIPSO and the French PARASOL. CloudSat was selected as a NASA Earth System Science Pathfinder satellite mission in 1999 to provide observations necessary to advance our understanding of cloud abundance, distribution, structure, and radiative properties. Since 2006, CloudSat has flown the first satellite-based millimeter-wavelength cloud radar—a radar that is more than 1000 times more sensitive than existing weather radars.

Unlike ground-based weather radars that use centimeter wavelengths to detect raindrop-sized particles, CloudSat’s radar allows us to detect the much smaller particles of liquid water and ice that constitute the large cloud masses that make our weather.

CloudSat was co-manifested with the CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations) satellite aboard a Delta II rocket for its launch on 28 April 2006. The satellites fly in a nearly circular orbit with an equatorial altitude of approximately 705 km. The orbit is sun-synchronous, maintaining a roughly fixed angle between the orbital plane and the mean solar meridian. CloudSat maintains a close formation with Aqua and a particularly close formation with CALIPSO, providing near-simultaneous and collocated observations with the instruments on these two platforms.