For Wednesday's assignment you wrote a script subset.py that defined the output_h5 function and wrote test_output.h5 Now read channel1 and channel 31 back in and use them to calculate the cloud fraction. (Recall that MODIS channels 1 and 31 are centered at about 0.64 microns and 11 microns, respectively.) This notebook assumes that you have converted your MYD021KM and MYD03 files to test_output. To reproduce my results go to http://clouds.eos.ubc.ca/~phil/Downloads/a301 and get the file A2005188.h5 and put it in a folder called data, at the same level a lib and notebooks
In [0]:
from __future__ import print_function
from __future__ import division
import os,site
import glob
import h5py
import numpy as np
from matplotlib import pyplot as plt
#
# add the lib folder to the path assuming it is on the same
# level as the notebooks folder
#
libdir=os.path.abspath('../lib')
site.addsitedir(libdir)
from modismeta_h5 import parseMeta
from h5dump import dumph5
In [0]:
%matplotlib inline
the glob function finds a file using a wildcard to save typing (google: python glob wildcard)
In [0]:
h5_filename=glob.glob('../data/A2005188.h5')[0]
print("found subsetted file {}".format(h5_filename))
In [0]:
h5_file=h5py.File(h5_filename)
Read back the arrays I wrote with output_h5
In [0]:
chan1=h5_file['channel1'][...]
chan31=h5_file['channel31'][...]
lats=h5_file['lattitude'][...]
lons=h5_file['longitude'][...]
In [0]:
c=2.99792458e+08 #m/s -- speed of light in vacumn
h=6.62606876e-34 #J s -- Planck's constant
kb=1.3806503e-23 # J/K -- Boltzman's constant
c1=2.*h*c**2.
c2=h*c/kb
def planckInvert(wavel,Llambda):
"""input wavelength in microns and Llambda in W/m^2/micron/sr, output
output brightness temperature in K (note that we've remove the factor
of pi because we are working with radiances, not fluxes)
"""
Llambda=Llambda*1.e6 #convert to W/m^2/m/sr
wavel=wavel*1.e-6 #convert wavelength to m
Tbright=c2/(wavel*np.log(c1/(wavel**5.*Llambda) + 1.))
return Tbright
In [0]:
Tb=planckInvert(11.05,chan31)
print(Tb)
turn the repoj_L1B notebook function into a module:
In [0]:
from reproject import reproj_L1B
reprojecting onto a regular lat/lon grid takes a long time. To speed things up, we save the gridded variables in an npz file if write_grid=True. If we don't need to change the grid, we can set write_grid=False and just read the gridded files in.
In [0]:
write_grid=True
if os.path.exists('keepgrid.npz'):
print('found keepgrid.npz, reusing gridded fields')
write_grid=False
if write_grid:
print('writing new keepgrid.npz file')
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)
C1_grid, longitude, latitude, bin_count = reproj_L1B(chan1, lons, lats, xlim, ylim, res)
#
# save these arrays in a dictionary and write out to keepgrid.npz
#
keep_grids=dict(Tb_grid=Tb_grid,C1_grid=C1_grid,longitude=longitude,latitude=latitude,bin_count=bin_count)
np.savez('keepgrid.npz',**keep_grids)
else:
grid_dict=np.load('keepgrid.npz')
latitude=grid_dict['latitude'];longitude=grid_dict['longitude'];bin_count=grid_dict['bin_count'];
Tb_grid=grid_dict['Tb_grid'];C1_grid=grid_dict['C1_grid']
As a first try: bin the longwave radiances into 100 bins and the shortwave reflectances into 150 bins
In [0]:
# create bins for brightness temp
Tb_bins = 50
Tb_lim=[np.nanmin(Tb_grid), np.nanmax(Tb_grid)]
Tb_bins=np.linspace(Tb_lim[0], Tb_lim[1], Tb_bins, dtype=np.float)
# and channel-1 reflectance
C1_bins = 50
C1_lim=[np.nanmin(C1_grid), np.nanmax(C1_grid)]
C1_bins=np.linspace(C1_lim[0], C1_lim[1], C1_bins, dtype=np.float)
Here, we define channel-1 data on x-axis and call np.histogram2d to get the bin_count value x_edges and y_edges.
Use some print statments to find the dimensions of H, X, and Y. In particular, make sure you understand
how np.meshgrid is turning two 1-dimensional arrays into two 2-dimensional arrays.
Try something like:
np.meshgrid([1,2,3,4],[-1,-2,-3,-4,-5,-6]) to see what it does
In [0]:
y=Tb_grid.flat[:]; y_bins=Tb_bins # x: Tb
x=C1_grid.flat[:]; x_bins=C1_bins # y: C1
H, y_edges, x_edges = np.histogram2d(y, x, bins=(y_bins, x_bins))
X, Y = np.meshgrid(x_edges[:-1], y_edges[:-1])
Then we make 2-D histogram to see the difference between clouds and ocean, the core idea is:
# 2-D histogram
ax.contourf(X, Y, H/np.max(H)) # use percentage, because H depends on the resolution 'res' we used before.
# try to distinguish clouds from ocean through linear function
# x is channel-1
We will display the 2-D histogram as the central plot, with 1-d histograms positioned to the left and below. We want an approximately 9-color colormap showing percentages in increments of 10% between our min and max
In [0]:
# make_axes_locatable ---> for axis control
from mpl_toolkits.axes_grid1 import make_axes_locatable
# set axis locations for three different plots
# a main contour plot and two histograms below and on
# the left hand side
#
left=0.1; width = 0.8; bottom=0.1; height = 0.65
gap=0.02; hist_len=0.2; cbar_len=0.12
# location of the three axes -- located by trial and error
#
#
rect_main = [left+hist_len+gap, bottom, width, height]
rect_histx = [left+hist_len+gap, left+height+gap, width-cbar_len, hist_len]
rect_histy = [left, bottom, hist_len, height]
# clev
#clevs=range(40, 281, 40)
#
# set up a color map for the histogram percentages, again
# trial and erro shows that about 10 bins between 3 and 30%
# looks attractive
#
# also assign white to percentages below 3%
#
#
clevs=np.arange(3, 31, 3)
CMap=plt.cm.PuBu
CMap.set_over(CMap(np.arange(256))[-1, 0:3])
CMap.set_under('w')
#
xlim_bin=[np.min(X), np.max(X)]
ylim_bin=[np.min(Y), np.max(Y)]
#
# ========== figure ========== #
fig=plt.figure(figsize=(9, 9))
# ========== Central histogram here ========== #
# axis
axMain=plt.axes(rect_main)
axMain.yaxis.tick_right()
axMain.yaxis.set_label_position('right')
axMain.set_xlim(xlim_bin)
axMain.set_ylim(ylim_bin)
axMain.set_xlabel('Channel-1 reflectivity', fontsize=12)
axMain.set_ylabel('Channel-31 brightness temperature (K)', fontsize=12)
axMain.set_title('2-D Histogram -- MODIS 0.6 and 11 microns', fontsize=16, fontweight='bold', x=1.15, y=1.15)
#
# divider is a container used locate the colorbar to the right of the 2-d histogram
# it knows where axMain is located on the page
#
divider=make_axes_locatable(axMain)
# grid and frame
plt.grid() # grid on
#
# make the axes thicker than default
#
[i.set_linewidth(2) for i in axMain.spines.itervalues()] # a bold frame
#
# plot the main histogram with the counts converted to percent
#
CS=axMain.contourf(X, Y, H/np.max(H)*100, clevs, cmap=CMap, extend='both') # 2-D histgram
CAx=divider.append_axes('right', size='5%', pad=0.75)
CBar=plt.colorbar(CS, cax=CAx)
CBar.set_label('Percentage ( % )', fontsize=10)
CBar.ax.tick_params(axis='y', length=22.5)
# draw line showing the boundary between cloudy and clear pixels
#
# slope and intercept chosen by eye
# the line Tbright=reflectivity*20./0.2 + 275. seems to separate the
# main warm cluster from the colder clouds
#
axMain.plot(x_edges, x_edges*20./0.2 + 275., \
color='k', linestyle='--', linewidth=5)
axMain.text(0.2, 295., 'Ocean', fontsize=16, fontweight='bold', \
ha='center', va='center', color='k')
axMain.text(0.12, 260, 'Cloud', fontsize=16, fontweight='bold', \
ha='center', va='center', color='k')
#
# ========== Hist-x on bottom axis========== #
# place a 1-d histogram below the 2-d histogram to show
# channel 1 reflectivities
#
axHistx=plt.axes(rect_histx)
axHistx.hist(x, bins=x_bins, color=[0.3, 0.6, 0.8])
axHistx.set_xlim(xlim_bin)
axHistx.axes.get_xaxis().set_visible(False)
# scientific notation for x, y-axis
plt.ticklabel_format(style='sci', axis='both', scilimits=(0,0))
[i.set_linewidth(2) for i in axHistx.spines.itervalues()]
#
# ========== Hist-y ========== #
#
# place a 1-d histogram to the left of the 2-d histogram to show channel 31
#
axHisty = plt.axes(rect_histy)
axHisty.hist(y, bins=y_bins, color=[0.3, 0.6, 0.8], orientation='horizontal')
axHisty.set_ylim(ylim_bin)
axHisty.invert_xaxis()
axHisty.axes.get_yaxis().set_visible(False)
plt.ticklabel_format(style='sci', axis='both', scilimits=(0,0))
out=[i.set_linewidth(2) for i in axHisty.spines.itervalues()]
# savefig
#plt.savefig('01_MODIS_L1B_histgram.png', dpi=450, facecolor='w', edgecolor='w',
# orientation='portrait', papertype='letter', format='png',
# transparent=True, bbox_inches='tight', pad_inches=0,
# frameon=None)
In [0]:
clear_cloudy=np.empty_like(Tb_grid)
clear_cloudy[...]=np.nan
for i in range(Tb_grid.shape[0]):
for j in range(Tb_grid.shape[1]):
#dashed line: Tbright=reflectivity*20./0.2 + 275.
if(Tb_grid[i, j] > C1_grid[i, j]*20/0.2 + 275.):
#print(i,j,Tb_grid[i,j])
clear_cloudy[i, j]=Tb_grid[i, j]
open_water=clear_cloudy[~np.isnan(clear_cloudy)]
open_fraction=open_water.size/clear_cloudy.size
print("Fraction of scene that is open water is: {:5.3f}".format(open_fraction))
out=plt.hist(clear_cloudy.flat[~np.isnan(clear_cloudy.flat)])
In [0]:
from mpl_toolkits.basemap import Basemap
lcc_values=dict(resolution='l',projection='lcc',
lat_1=20,lat_2=40,lat_0=30,lon_0=135,
llcrnrlon=120,llcrnrlat=20,
urcrnrlon=150,urcrnrlat=42)
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)
#
# the cloudy pixels have values of np.nan.
# find these and set them to 200 K for
# visualization
#
hit=np.isnan(clear_cloudy)
clear_cloudy[hit]=200.
CS=proj.pcolor(x, y, clear_cloudy, cmap=plt.cm.hot)
# colorbar
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('Sea surface temperature (K)', fontsize=10)
CBar.ax.tick_params(axis='y', length=0)
In [0]: