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.
In [0]:
from __future__ import print_function
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/test_output.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]:
with h5py.File
chan1=h5_file['channel1'][...]
chan31=h5_file['channel31'][...]
lats=h5_file['lattitude'][...]
lons=h5_file['longitude'][...]
turn the repoj_L1B notebook function into a module:
In [0]:
from reproject import reproj_L1B
In [0]:
xlim=[np.min(lons), np.max(lons)]
ylim=[np.min(lats), np.max(lats)]
C31_grid, longitude, latitude, bin_count = reproj_L1B(chan31,lons, lats, xlim, ylim, 0.1)
C1_grid, longitude, latitude, bin_count = reproj_L1B(chan1,lons, lats, xlim, ylim, 0.1)
In [0]:
C31_grid=np.ma.masked_where(np.isnan(C31_grid), C31_grid)
bin_count=np.ma.masked_where(np.isnan(bin_count), bin_count)
longitude=np.ma.masked_where(np.isnan(longitude), longitude)
latitude=np.ma.masked_where(np.isnan(latitude), latitude)
longitude.shape
**`numpy.histogram2d`** is the main function we use here to create a 2-D histogram, it partitions two 1-D arrays into two 1-D binned arrays, and returns 2-D counts in a matrix formed by sqauares consisting of the intersection of the two sets of bins, as well as the 2-D bin edges.
The I/O format of numpy.histogram2d is a little obscure -- the arguments look like:
H, y_edges, x_edges = np.histogram2d(y, x, bins=(y_bins, x_bins))
X, Y = np.meshgrid(x_edges[:-1], y_edges[:-1]) # '-1' because number_bins=number_data-1
numpy.histogram2d is different from
**`numpy.digitize`** used for regridding. numpy.digitize does not returns counts in each bin and we have to do this in a for loop (as we did in our function reproj_L1B).
There is a counterpart of numpy.histogram2d for 1-D histogram named
**`numpy.histogram`**. One can also make histograms through **`pyplot.hist`** and **`pyplot.hist2d`**.
In [0]:
# create bins for channel-31
C31_bins = 100
C31_lim=[np.nanmin(C31_grid), np.nanmax(C31_grid)]
C31_bins=np.linspace(C31_lim[0], C31_lim[1], C31_bins, dtype=np.float)
# and channel-1
C1_bins = 150
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)
In [0]:
fig=plt.figure(figsize=(10.5, 9.5))
ax=plt.gca()
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
image=ax.pcolormesh(longitude, latitude, C31_grid)
In [0]:
#got to 0.05 degrees for higher resolution
res=0.05
xlim=[np.min(lons), np.max(lons)]
ylim=[np.min(lats), np.max(lats)]
C31_grid, longitude, latitude, bin_count = reproj_L1B(chan31, lons, lats, xlim, ylim, res)
C1_grid, longitude, latitude, bin_count = reproj_L1B(chan1, lons, lats, xlim, ylim, res)
As a first try: bin the longwave radiances into 100 bins and the shortwave reflectances into 150 bins
In [0]:
# create bins for channel-31
C31_bins = 100
C31_lim=[np.nanmin(C31_grid), np.nanmax(C31_grid)]
C31_bins=np.linspace(C31_lim[0], C31_lim[1], C31_bins, dtype=np.float)
# and channel-1
C1_bins = 150
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=C31_grid.flat[:]; y_bins=C31_bins # x: C31
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 radiance ($W/m^2/\mu m/sr$)', 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
#
axMain.plot(x_edges, x_edges*4.5+7.1, \
color='k', linestyle='--', linewidth=5)
axMain.text(0.4, 6.25, 'Cloud', fontsize=16, fontweight='bold', \
ha='center', va='center', color='k')
axMain.text(0.125, 8.0, 'Ocean', 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]: