getting the cloud fraction from channels 1 and 31 -- brightness temperature

This is a modification of cloud_fraction.ipynb to use brightness temperature instead of Channel 31 radiance

For Wednesday's assignment you wrote a script 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 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
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]:
print("found subsetted file {}".format(h5_filename))

In [0]:

Read back the arrays I wrote with output_h5

In [0]:

find the brightness temperature for the chan31 radiances at 11.05 microns

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

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]:

turn the repoj_L1B notebook function into a module:

In [0]:
from reproject import reproj_L1B

Two dimensional histograms using histogram2d

writing the gridded files out to be reused in subsequent runs (memoization)

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]:
if os.path.exists('keepgrid.npz'):
    print('found keepgrid.npz, reusing gridded fields')
if write_grid:
    print('writing new keepgrid.npz file')
    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

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.set_over(CMap(np.arange(256))[-1, 0:3])
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.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
# 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)'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.hist(x, bins=x_bins, color=[0.3, 0.6, 0.8])
# 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')
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)

Now find all the pixels that lie above our "decision boundary" and count them to get open water fraction of scene

In [0]:
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.):
           clear_cloudy[i, j]=Tb_grid[i, j]
print("Fraction of scene that is open water is: {:5.3f}".format(open_fraction))

Finally, make a map of the open water, setting the cloud temp to 200 K for a "don't care" value

In [0]:
from mpl_toolkits.basemap import Basemap
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
## 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
CS=proj.pcolor(x, y, clear_cloudy,
# colorbar
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('Sea surface temperature (K)', fontsize=10)'y', length=0)

In [0]: