This notebook plots vertical cross sections through a cyclone of $Q_R$ the longwave and shortwave heating rate in K/hour. It uses the level 2B product FLXHR, which is described at the cloudsat website as follows:
"This algorithm derives estimates of broadband fluxes and heating rates consistent with liquid and ice water content estimates from the CloudSat Profiling Radar (CPR). For each radar profile, a broadband radiative transfer model is used to calculate upwelling and downwelling longwave and shortwave fluxes at each CPR range gate from the surface to the lower stratosphere. Profiles of cloud ice and liquid water content and cloud particle effective radii are defined based on the CloudSat 2B-LWC and 2B-IWC products while precipitation properties are defined using the CloudSat 2C-PRECIP-COLUMN dataset. Ancillary atmospheric state variables are interpolated from ECMWF analyses and surface albedos are assigned based on seasonally-varying maps of surface reflectance properties in combination with daily snow and sea ice cover maps from passive microwave instruments. Equivalent clear sky radiative flux profiles are generated by removing all clouds and repeating the calculations. Corresponding profiles of atmospheric heating are inferred from the vertical derivative of these fluxes."
Format is described on the Cloudsat web site
Units: K/hr
Variable name: Qr
Shape: Qr[2, 37082, 125] where Qr[0,37082,125] is the shortwave heating rate and Qr[1,37082,125] is the longwave heating rate. The other two dimensions are the same as the radar reflectivity: there are 37082 radar measurements in an orbit, binned into 125 vertical height bins
In [1]:
import h5py
import numpy as np
import datetime as dt
from datetime import timezone as tz
import matplotlib
from matplotlib import pyplot as plt
import pyproj
from numpy import ma
from a301utils.a301_readfile import download
from a301lib.cloudsat import get_geo
from IPython.display import Image, display
from datetime import datetime,timezone
flx_file='2008082060027_10105_CS_2B-FLXHR_GRANULE_P2_R04_E02.npz'
download(flx_file)
with np.load(flx_file) as npz:
lons=npz['lons']
lats=npz['lats']
height=npz['height']
shortwave_hr=npz['shortwave_hr']
longwave_hr=npz['longwave_hr']
date_times=npz['date_times']
prof_times=npz['prof_times']
date_times=[datetime.fromtimestamp(item,timezone.utc) for item in date_times]
date_times=np.array(date_times)
meters2km=1.e3
In [2]:
long_wave_hr=np.ma.masked_invalid(longwave_hr)
short_wave_hr=np.ma.masked_invalid(shortwave_hr)
In [3]:
first_time=date_times[0]
print('orbit start: {}'.format(first_time))
start_hour=6
start_minute=45
storm_start=starttime=dt.datetime(first_time.year,first_time.month,first_time.day,
start_hour,start_minute,0,tzinfo=tz.utc)
#
# get 3 minutes of data from the storm_start
#
storm_stop=storm_start + dt.timedelta(minutes=3)
print('storm start: {}'.format(storm_start))
time_hit = np.logical_and(date_times > storm_start,date_times < storm_stop)
storm_lats = lats[time_hit]
storm_lons=lons[time_hit]
storm_prof_times=prof_times[time_hit]
storm_sw_hr=shortwave_hr[time_hit,:]
storm_lw_hr=longwave_hr[time_hit,:]
storm_height=height[time_hit,:]
storm_date_times=date_times[time_hit]
len(date_times)
Out[3]:
In [4]:
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
i.e. assume that height[0,:] = height[1,:] = ...
in reality, the bin heights are depend on the details of the radar returns, so we would need to historgram the heights into a uniform set of bins -- ignore that for this qualitative picture
In [5]:
%matplotlib inline
plt.close('all')
from matplotlib import cm
from matplotlib.colors import Normalize
vmin=-30
vmax=30
the_norm=Normalize(vmin=vmin,vmax=vmax,clip=False)
cmap_ref=cm.RdBu_r
cmap_ref.set_over('pink')
cmap_ref.set_under('k')
cmap_ref.set_bad('0.75') #75% grey
#
# Q-1: What is the difference between the distance,height,field,ax
# and cmap,norm arguments to this function? Why do I structure
# the function signature this way?
#
def plot_field(distance,height,field,ax,cmap=None,norm=None):
if cmap is None:
cmap=cm.jet
col=ax.pcolormesh(distance,height,field,cmap=cmap,
norm=the_norm)
cax=fig.colorbar(col,extend='both',ax=ax,pad= 0.01)
return ax,cax
fig,(ax1,ax2)=plt.subplots(2,1,figsize=(20,10))
cloud_height_km=height[0,:]/meters2km
ax1,cax1=plot_field(distance,cloud_height_km,storm_sw_hr.T,ax1,cmap=cmap_ref,
norm=the_norm)
ax2,cax2=plot_field(distance,cloud_height_km,storm_lw_hr.T,ax2,cmap=cmap_ref,
norm=the_norm)
for colorbar in [cax1,cax2]:
text=colorbar.set_label('heating rate (K/hr)',rotation=-90,verticalalignment='bottom')
for ax in [ax1,ax2]:
ax.set(ylim=[0,17],xlim=(0,1200))
ax.set_xlabel('distance (km)',fontsize=15)
ax.set_ylabel('height (km)',fontsize=15)
text=fig.suptitle('storm radiative heating rates: shortwave (top), longwave (bottom)',size=25)
fig.savefig('heating_rates.png',dpi=100)
In [ ]: