In this notebook I show how to find a segment of an orbit using start and stop times, and make some plots that are rough and ready -- not as pretty as the ones in orbit_tools, but good for quick looks at the data
It uses Cloudsat and MODIS subset data for Day 2006303 (Oct. 30, 2006), which can be downloaded from
http://clouds.eos.ubc.ca/~phil/Downloads/a301
and
In [0]:
from __future__ import division
from __future__ import print_function
import site
site.addsitedir('../lib/cloudsat')
site.addsitedir('../lib')
import glob
import h5py
import numpy as np
import datetime
import dateutil.tz as tz
import importlib
#
from matplotlib import pyplot as plt
import orbit_tool as OT
import cloudsat_tool as CT
importlib.reload(OT)
importlib.reload(CT)
#help(plt)
In [0]:
%matplotlib inline
First open the ECMWF reanalysis file that contains pressure, temperature etc. along the swath from the European Center forecast model
In [0]:
hdf5_EC=glob.glob('../data/20063032*ECMWF-AUX*.h5')[0]
print("ECMWF-AUX file found\n{}".format(hdf5_EC))
It's a good idea to check to see what the missing values are, and whether the data is integer or float (in this case they are 4 byte (32 bit) floats)
In [0]:
with h5py.File(hdf5_EC) as obj:
P=obj['ECMWF-AUX/Data Fields/Pressure']
missing_val=P.attrs['missing']
print("missing value: ",missing_val)
print("pressure dtype: ",P.value.dtype)
P=P.value
print("dimensions of P are: {}".format(P.shape))
Check to see what fraction of the pressure field is missing
In [0]:
hit=(P == missing_val)
print('found {} missing values out of {} total'.format(np.sum(hit),P.size))
#
# mask out the missing values
#
P=np.ma.masked_where(P == missing_val, P)
use the "~" symbol (logical not, like Matlab) to get the good values that are note masked and histogram them. Note that indexing by the mask flattens the array to 1 dimension
In [0]:
tryit=P[~P.mask]
print(P.shape,tryit.shape)
out=plt.hist(P[~P.mask])
In [0]:
P, SLP, T, T2m, SKT, q, O3 = CT.read_ecmwf(hdf5_EC, maskid=2)
To compare do the same thing for the return values of the read_ecmwf function. The pressure is converted from Pa to hPa
In [0]:
out=plt.hist(P[~P.mask])
In [0]:
lat, lon, date_day, prof_sec, dem_elev = CT.get_geo(hdf5_EC, monotonic_id=1)
Note in the next cell that I resize the figure
In [0]:
proj=OT.draw_orbit_pha(hdf5_EC, lonlim=None, latlim=None, time_step=625, label_step=2)
proj.ax.figure.set_size_inches([10,10])
From the orbit plot we can see that we are near Vancouver at about 22:24 UCT on October 30, 2006. Create start and stop datetimes between 22:21 and 22:24 UCT and use those to segment the rest of the data
In [0]:
start=datetime.datetime(2006,10,30,22,21,0,tzinfo=tz.tzutc())
stop=datetime.datetime(2006,10,30,22,24,0,tzinfo=tz.tzutc())
hit=np.logical_and(date_day > start,date_day < stop)
In [0]:
out=plt.plot(lon[hit],lat[hit])
title=plt.gca().set_title('quick plot to check to see lat/lon are correct')
In [0]:
hdf5_radar=glob.glob('../data/20063032*CS_2B-GEOPROF_GRANULE*.h5')[0]
print("geoprof radar file found\n{}".format(hdf5_radar))
Now plot the radar reflectivity using the pcolor (pseudocolor) routine
In [0]:
height, reflect = CT.read_radar(hdf5_radar, maskid=2)
color_map=plt.cm.jet
height.shape
lon[hit].shape
print(reflect.shape,height[hit,:].shape)
pcolor requires a set of heights and a set of lons that are the same shape as the reflectivity: 1125 values by 125 vertical levels. Use np.meshgrid for this. Note that the heights don't change along the swath, so just take the values for the first radar pulse (row 0)
In [0]:
height_seg=height[hit,:]
lon_seg=lon[hit]
lat_seg=lat[hit]
time_seg=date_day[hit]
height_mesh, lon_seg_mesh=np.meshgrid(height_seg[0,:], lon_seg)
fig=plt.figure(1,figsize=[10,4])
ax=fig.add_subplot(111)
image=ax.pcolor(lon_seg_mesh,height_mesh,reflect[hit,:],cmap=color_map)
ax.set_ylim([0,10])
cbar=plt.colorbar(image)
cbar.set_label('reflectivity (dbZ)')
ax.set_ylabel('height (km)')
out=ax.set_xlabel('longitude (deg W)')
In [0]:
hdf5_rain=glob.glob('../data/20063032*2C-RAIN-PROFILE*.h5')[0]
rain, precli, precice, clw = CT.read_rain(hdf5_rain, maskid=2)
In [0]:
rain_seg=rain[hit]
fig=plt.figure(1,figsize=[10,4])
ax=fig.add_subplot(111)
from matplotlib.dates import DateFormatter
ax.xaxis.set_major_formatter(DateFormatter('%H:%M'))
ax.plot(time_seg,rain_seg)
ax.set_ylabel('rain rate (mm/hour)')
Finally, plot temperature soundings at the beginning and end of the segment
In [0]:
fig=plt.figure(1,figsize=[8,6])
ax=fig.add_subplot(111)
the_temps=T[hit,:]
z=height[0,:]
start=date_day[hit][0].strftime("%H:%M UCT")
stop=date_day[hit][-1].strftime("%H:%M UCT")
ax.plot(the_temps[0,:],z,'b',label=start)
ax.plot(the_temps[-1,:],z,'r',label=stop)
ax.legend(loc='best')
ax.set_ylabel('height (km)')
ax.set_ylim([0,25])
ax.set_title('ECMWF start/stop soundings')
out=ax.set_xlabel('ECMWF temperature (deg C)')
I converted the MODIS level1b and geom file for the granule starting at 22:20 UCT using subset.output_h5
In [18]: from subset import output_h5
In [19]: level1b=glob.glob('MYD021KM*2006*h5')[0]
In [20]: geom=glob.glob('MYD03*2006*h5')[0]
In [21]: output='A2006303_subset.h5'
In [22]: output_h5(level1b,geom,output)
You can download this file: http://clouds.eos.ubc.ca/~phil/Downloads/a301/A2006303_subset.h5
I'm assuming it is in the folder ../data
In [0]:
subset=glob.glob('../data/A2006303_subset.h5')[0]
In [0]:
from reproject import reproj_L1B
from planck import planckInvert
In [0]:
with h5py.File(subset) as h5_file:
chan1=h5_file['channel1'][...]
chan31=h5_file['channel31'][...]
lats=h5_file['lattitude'][...]
lons=h5_file['longitude'][...]
chan31=chan31*1.e6 #planckInvert needs W/m^2/m/sr, Modis is W/m^2/micron/sr
In [0]:
Tb=planckInvert(11.05*1.e-6,chan31) #needs wavelength in meters
res=0.05
xlim=[np.min(lons), np.max(lons)]
ylim=[np.min(lats), np.max(lats)]
Tb_grid, longitude, latitude, bin_count = reproj_L1B(Tb, lons, lats, xlim, ylim, res)
In [0]:
out=plt.hist(Tb_grid[~np.isnan(Tb_grid)])
In [0]:
from mpl_toolkits.basemap import Basemap
lcc_values=dict(resolution='l',projection='lcc',
lat_1=30,lat_2=50,lat_0=45,lon_0=-135,
llcrnrlon=-150,llcrnrlat=30,
urcrnrlon=-120,urcrnrlat=50)
proj=Basemap(**lcc_values)
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
ax=fig.add_subplot(111)
## define parallels and meridians to draw.
parallels=np.arange(-90, 90, 5)
meridians=np.arange(0, 360, 5)
proj.drawparallels(parallels, labels=[1, 0, 0, 0],\
fontsize=10, latmax=90)
proj.drawmeridians(meridians, labels=[0, 0, 0, 1],\
fontsize=10, latmax=90)
# draw coast & fill continents
#map.fillcontinents(color=[0.25, 0.25, 0.25], lake_color=None) # coral
out=proj.drawcoastlines(linewidth=1.5, linestyle='solid', color='k')
x, y=proj(longitude, latitude)
CS=proj.pcolor(x, y, Tb_grid, cmap=plt.cm.hot)
print(Tb_grid.shape)
#
# now add the radar ground track to the image
#
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('Channel 31 brightness temperature (K)')
radarx,radary=proj(lon_seg,lat_seg)
out=proj.plot(radarx,radary,'b-',lw=10)
title=ax.set_title('MODIS Brightness Temperature (K) for Day 2006303 at 22:20 UCT')
In [0]: