• The objectives of this notebook is to provide examples about working with AQUA MODIS L1B files, including:
    • Convert HDF4 format data to HDF5 format;
    • Read *.HDF5 data through "h5py" module;
    • Correct the raw data by using scale_factor and offset
    • Reproject sinusoidal projection grid to a regular grid (cylindrical equidistant projection);
    • Identify cloud and ocean based on 2-D histograms;
    • Calculate brightness temperature on the top of the clouds and save the result as *.mat file;
    • Visualizing cloud data through mpl_toolkits.basemap.
  • The files required for this notebook including:
  • MYD021KM* are MODIS AQUA L1B files version 6, observed on Oct 30 2006. Which shows a extra-tropical cyclone approaching the Vancouver Island.
  • MYD03* is the geolocation file indecates the MODIS AQUA L1B geolocation field
  • *.hdf is the HDF4 file which can be downloaded from LAADS Web.
  • *.h5 is the HDF5 file which can be converted from HDF4 files.

In [0]:
__author__ = 'ATSC-301 UBC'


Convert HDF4 format file to HDF5 format

The hdf4 module in Python, pyhdf is not included in anaconda channel for both Windows and Mac. So here we plan to convert HDF4 into HDF5 or netCDF format and solve the problem.

But one can also try to install this module manually, the source code and documents can be get from here. Linux user can try to install python module mutirri , which includes pyhdf in it.

Convert HDF4 to HDF5 on Windows

Windows user can convert from HDF4 to HDF5 using h4toh5 which provided by hdf group.

One need to do the following steps:

  • Download h4h5tools-2.2.2-win32.zip or h4h5tools-2.2.2-win64.zip on local machine.
  • Uncompress the *.zip file, run installation file *.exe
  • Open DOS command line window
  • cd to bin folder (e.g. cd F:\HDF_Group\H4TOH5\2.2.2\bin)
  • use h4toh5 tool in the following ways:
    h4toh5convert input.hdf output.h5
          h4toh5convert -na input.hdf output.h5
          h4toh5convert -nc4 input.hdf output.h5`
  • You can add h4toh5 on your system's Environmental Variable PATH, or set PATH=F:\HDF_Group\H4TOH5\2.2.2\bin and then you can use h4toh5convert everywhere.
  • Here is a screen shot of h4toh5convert in cmd:

Convert HDF4 to HDF5 on Linux

h4toh5 has a Linux version, so we can use it by repeating the steps above and edit PATH in .bashrc file.

Import modules & tools

Before you run the cell below, you need to install h5py through anaconda:

conda install h5py

Here we use h5py to read HDF5 file. The h5py package is a open-source Pythonic interface to the HDF5 binary data format, one can find something useful via h5py's documentation and their website. Also, one can find more information about HDF5's file format from here.

In [0]:
import glob
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
from __future__ import print_function
% matplotlib inline

We import h5py to read HDF5 files:

In [0]:
import h5py

scipy.io for saving data in *.mat format

In [0]:
import scipy.io

For the map view of data, we need mpl_toolkits.basemap

In [0]:
from mpl_toolkits.basemap import Basemap

Read MODIS L1B data through h5py

Searching all *.h5 files in _data/MODIS_L1B/ directory, Make sure that you have the required HDF5 files exist in it.

In [0]:
print("MODIS L1B file found {}".format(hdf5_L1B))
print("MODIS Geolocation file found {}".format(hdf5_Geo))

Get HDF5 file object

In [0]:
hdf5_obj=h5py.File(hdf5_L1B[0], 'r')
geo_obj=h5py.File(hdf5_Geo[0], 'r')

h5py.File('filename.h5', 'r') returns the object of HDF5 file.

Here we read MODIS L1B channe-31 and channel-1 data from MYD021KM.A*.h5 file and read geolocation fields from a corresponding MYD02QKM.A*.h5 file.

The structure of HDF file

Get attribute and sub-attribute information

Something abou using h5py to discover the structure of HDF files:

  • Apply a_obj.keys() to get all the attributes in a file object.
  • HDF file usually has more than one layer attributes.
  • Call a_obj['an_attribute'].keys() to see the sub-attribute.
  • Use build-in function dir(a_obj['an_attribute']) to see if the option keys still available.
  • If you cannot find keys by dir(a_obj['an_attribute']['more_attr']), you reach the bottom layer of the file, and usually there are data in it.

Here we attempt to read MODIS L1B channe-31 data as an example. The data 'EV_1KM_Emissive' is in:

In [0]:
print('Attributes in {}'.format(hdf5_L1B))
print('hdf5_attr=hdf5_obj.keys() \n\n{}'.format(hdf5_attr))

print("\n\n\tSub-attributes in 'MODIS_SWATH_Type_L1B'")
print("\tsub_attr=hdf5_obj['MODIS_SWATH_Type_L1B'].keys() \n\n\t{}".format(sub_attr))

subsub_attr=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields'].keys()
print("\n\n\t\tSubsub-attributes in 'Data Fields'")
print("\t\tsubsub_attr=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields'].keys() \n\n\t\t{}".format(subsub_attr))

Some external tools

If you have tools like HDF Explorer or HDF Viewer. You may have a more direct view of HDF's file structure.

Read MODIS L1B channel-31 and channel-1 data

Read raw data

Based on chapter 3.B, we read channe-31 and channel-1 data.

  • For channel-31, it is in 'EV_1KM_Emissive' (from channel-20 to 36, 16 different channels)
      'MODIS_SWATH_Type_L1B' '/' 'Data Fields' '/' 'EV_1KM_Emissive'
  • For channel-1, it is in 'EV_250_Aggr1km_RefSB' (channel-1 and 2, 2 different channels)
      'MODIS_SWATH_Type_L1B' '/' 'Data Fields' '/' 'EV_250_Aggr1km_RefSB'

In [0]:
# Channel-31
L1B_emiss=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'][:];
print("Size of 'EV_1KM_Emissive':\n===========================\n{}".format(L1B_emiss.shape))
# Channel-1
L1B_ref=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'][:];
print("\nSize of 'EV_500_Aggr1km_RefSB':\n================================\n{}".format(L1B_ref.shape))

Here the file has a size of Channels * Longitude * Latitude.

Select the channel

The channel information of 'EV_1KM_Emissive' and 'EV_500_Aggr1km_RefSB' can be found in:

    'MODIS_SWATH_Type_L1B' '/' 'Data Fields' '/' 'Band_1KM_Emissive'
    'MODIS_SWATH_Type_L1B' '/' 'Data Fields' '/' 'Band_250M'

In [0]:
band_info=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['Band_1KM_Emissive'][:]
print('List of MODIS L1B Channels\n=====================================\n{}'.format(band_info))
band_info=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['Band_250M'][:]
print('\nList of MODIS L1B Channels\n=====================================\n{}'.format(band_info))

Then we can chose the channel we want:

In [0]:
C31=L1B_emiss[10, :, :]
C1=L1B_ref[1, :, :]

Scale factor and offset value

Simply read raw data with channel is not enough, in order to maximum the precision, MODIS L1B data formed as

$Data = (RawData - offset) \times scale$

this information is included in the attributes of each variable.

Here we use obj[...].attrs.items() to get all the attributes and see if we can find something interesting.

In [0]:
print('Channel-31 info\n===============================================================')
hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs.items()

Here radiance_scales and radiance_offsets are what we want. Number of channels can also be seen through band_names.

We can use a_list=obj[...].attrs.values() to get these info.

In [0]:
a_list=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs.values()

radiance_scales and radiance_offsets are the 7th and 8th group of a_list, and channel-31 is the 11th element of the group.

In [0]:

We do the same thing for channel-1 data, but now we use reflectance_scales

In [0]:
C1_scale=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'].attrs.values()[9][0]
C1_offset=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'].attrs.values()[10][0]
#corrected_counts_scales=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_250_Aggr1km_RefSB'].attrs.values()[12][0]

Finally, we correct the data, numpy.ones is the same as ones.m in MATLAB.

In [0]:
C31=(C31 - C31_offset * np.ones(C31.shape))*C31_scale
C1=(C1 - C1_offset * np.ones(C1.shape))*C1_scale

Print the maximum data and see if it is reasonable.

In [0]:

Geolocation field

For the geolocation field, we do not using 'Geolocation Fields' in MYD021KM*.h5 file because 'EV_1KM_Emissive' and 'EV_250_Aggr1km_RefSB'as a resolution of 250m, but the resolution of latitude and longitude in 'Geolocation Fields' is 1km.

We use the 'Geolocation Field' in MYD03*.h5, they are in the following place:

    'MODIS_Swath_Type_GEO' '/' 'Geolocation Fields' '/' 'Longitude'
    'MODIS_Swath_Type_GEO' '/' 'Geolocation Fields' '/' 'Latitude'

In [0]:
C_x=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][:]
C_y=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][:]
print('Size of Longitude: {}'.format(C_x.shape))
print('Longitude Interval: {} ~ {}'.format(np.min(C_x), np.max(C_x)))
print('Size of Latitude: {}'.format(C_y.shape))
print('Latitude Interval: {} ~ {}'.format(np.min(C_y), np.max(C_y)))

Reproject MODIS L1B data to a regular grid

Here we define function reproj_L1B to reproject data to a regular grid.

One can also use MODIS Reprojection Tool to do the work.

Define function reproj_L1B

We now using a function to represent the reprojection process as what have done on satellite3.ipynb. One thing a bit different is that we use

numpy.searchsorted(b, a, 'right')

to replace

numpy.digitize(a, b)

Because `numpy.searchsorted` use binary search but `numpy.digitize` based on linear search, the later could be slow when we have lots of data.

In [0]:
def reproj_L1B(raw_data, raw_x, raw_y, xlim, ylim, res):
    Reproject MODIS L1B file to a regular grid
    d_array, x_array, y_array, bin_count = reproj_L1B(raw_data, raw_x, raw_y, xlim, ylim, res)
            raw_data: L1B data, N*M 2-D array.
            raw_x: longitude info. N*M 2-D array.
            raw_y: latitude info. N*M 2-D array.
            xlim: range of longitude, a list.
            ylim: range of latitude, a list.
            res: resolution, single value.
            d_array: L1B reprojected data.
            x_array: reprojected longitude.
            y_array: reprojected latitude.
            bin_count: how many raw data point included in a reprojected grid.
            function do not performs well if "res" is larger than the resolution of input data.
            size of "raw_data", "raw_x", "raw_y" must agree.
    import numpy as np
    x_bins=np.arange(xlim[0], xlim[1], res)
    y_bins=np.arange(ylim[0], ylim[1], res)
#    x_indices=np.digitize(raw_x.flat, x_bins)
#    y_indices=np.digitize(raw_y.flat, y_bins)
    x_indices=np.searchsorted(x_bins, raw_x.flat, 'right')
    y_indices=np.searchsorted(y_bins, raw_y.flat, 'right')
    y_array=np.zeros([len(y_bins), len(x_bins)], dtype=np.float)
    x_array=np.zeros([len(y_bins), len(x_bins)], dtype=np.float)
    d_array=np.zeros([len(y_bins), len(x_bins)], dtype=np.float)
    bin_count=np.zeros([len(y_bins), len(x_bins)], dtype=np.int)
    for n in range(len(y_indices)): #indices
        bin_row=y_indices[n]-1 # '-1' is because we call 'right' in np.searchsorted.
        bin_count[bin_row, bin_col] += 1
        x_array[bin_row, bin_col] += raw_x.flat[n]
        y_array[bin_row, bin_col] += raw_y.flat[n]
        d_array[bin_row, bin_col] += raw_data.flat[n]
    for i in range(x_array.shape[0]):
        for j in range(x_array.shape[1]):
            if bin_count[i, j] > 0:
                x_array[i, j]=x_array[i, j]/bin_count[i, j]
                y_array[i, j]=y_array[i, j]/bin_count[i, j]
                d_array[i, j]=d_array[i, j]/bin_count[i, j] 
                d_array[i, j]=np.nan
                x_array[i, j]=np.nan
    return d_array, x_array, y_array, bin_count

Test if reproj_L1B works well

Here we reproject channel-31 data to see if reproj_L1B works well.

In [0]:
xlim=[np.min(C_x), np.max(C_x)]
ylim=[np.min(C_y), np.max(C_y)]
C31_grid, longitude, latitude, bin_count = reproj_L1B(C31, C_x, C_y, xlim, ylim, 0.1)

Mask NaN for plot, also make sure that the data is not too big to plot.

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)

Plot the result

In [0]:
fig=plt.figure(figsize=(10.5, 9.5))
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
image=ax.pcolormesh(longitude, latitude, C31_grid)

Convert channel-31 and channel-1 data

Be careful to chose res here, high resolution makes your computer slow. Also, it affects the result on 2-D histgram part.

In [0]:
xlim=[np.min(C_x), np.max(C_x)]
ylim=[np.min(C_y), np.max(C_y)]
C31_grid, longitude, latitude, bin_count = reproj_L1B(C31, C_x, C_y, xlim, ylim, res)
C1_grid, longitude, latitude, bin_count = reproj_L1B(C1, C_x, C_y, xlim, ylim, res)

Identify cloud and ocean via 2-D histogram

Histogram is some kind of a basic image segmentation technique. Here we apply 2-D histogram to distinguish clouds from ocean.

How to use numpy.histogram2d

**`numpy.histogram2d`** is the main function we used here to create a 2-D histgram, it partitions two 1-D array into two 1-D bin array, and returns 2-D counts in a combination of 2 bins, as well as the 2-D bin edges.

The I/O format of numpy.histogram2d is not very clear, based on my understanding, a proper way is:

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`** what we used before. numpy.digitize do not returns counts in each bin and we have to do this in a for loop (as what we did in our function reproj_L1B).

There is a counter part of numpy.histogram2d for 1-D histogram named **`numpy.histogram2d`**. One can also make histograms through **`pyplot.hist`** and **`pyplot.hist2d`**.

MATLAB users can use function **`histcounts.m`** and **`histogram.m`** for 1-D histgram. There is no 2-D histgram function in MATLAB's official toolboxes, but one can find lots of them in MATLAB's file exchange center (e.g. here).

Create 2-D Histgram for channel-31 and channel-1 data

Create bins for channel-31 and channel-1

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 as what's in above section to get bin_count value x_edges and y_edges. Noted that masked NumPy array has no attribute flat.

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 histgram to see the difference between clouds and ocean, the core idea is:

# 2-D histgram
ax.contourf(X, Y, H/np.max(H)) # use percentage, because H sensitive to resolution 'res' we used before.
# try to distinguish clouds from ocean through linear function
# x is channel-1
axMain.plot(x, x*5.5+6.5*np.ones(x.shape))

The rest are codes for figures and axises.

In [0]:
# make_axes_locatable ---> for axis control
from mpl_toolkits.axes_grid1 import make_axes_locatable
# set axis
left=0.1; width = 0.8; bottom=0.1; height = 0.65
gap=0.02; hist_len=0.2; cbar_len=0.12
# three boxes
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)
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))
# ========== Main ========== #
# axis
axMain.set_xlabel('Channel-1', fontsize=12)
axMain.set_ylabel('Channel-31', fontsize=12)
axMain.set_title('2-D Histgram', fontsize=16, fontweight='bold', x=1.15, y=1.15)
# grid and frame
plt.grid() # grid on
[i.set_linewidth(2) for i in axMain.spines.itervalues()] # a bold frame
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
axMain.plot(x_edges, x_edges*5.5+6.5*np.ones(x_edges.shape), \
            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 ========== #
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 ========== #
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))
[i.set_linewidth(2) for i in axHisty.spines.itervalues()]
# savefig
plt.savefig('_figures/01_MODIS_L1B_histgram.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,
# show

We can see that, there are generally two place where data points are very dense, one is typical ocean, one is typical cloud.

One can set an arbitrary criteria to segment the image. An over strict criteria will wrongly eliminate the information of clouds, but a over loose one will view some of the ocean points as clouds. Here we use:

$data\left.\right|_{Channel-31} < data\left.\right|_{Channel-1} \times 5.5 + 6.5 $

to identify clouds in all the data points.

We repeat the line to all data points to make it clear.

In [0]:
fig=plt.figure(figsize=(8, 8))
ax.set_xlim(xlim_bin[0], xlim_bin[1])
ax.set_ylim(ylim_bin[0], ylim_bin[1])
ax.set_xlabel('Channel-1', fontsize=12)
ax.set_ylabel('Channel-31', fontsize=12)
ax.plot(x, y, color=[0.5, 0.5, 0.5], marker='.', linestyle='None')
ax.plot(x_edges, x_edges*5.5+6.5*np.ones(x_edges.shape), linestyle='--', color='k', linewidth=5)
ax.text(0.4, 6.25, 'Cloud', fontsize=16, fontweight='bold', \
                    ha='center', va='center', color='k')
ax.text(0.10725, 7.75, 'Ocean', fontsize=16, fontweight='bold', \
                    ha='center', va='center', color='k')
plt.savefig('_figures/01_MODIS_L1B_Divide_Cloud_and_Ocean.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,

Eliminate ocean points based on 2-D histgram

We replace the ocean data point to NaN based on the criteria above.

Using for loops and if command is a general way to do that, sometimes there will be more simple ways like:

id=C1*5.5+6.5; C31[C31<id]=np.nan

In [0]:
criteria_k=5.5 # less than
C1_clouds=np.empty((C31_grid.shape[0], C31_grid.shape[1],))
C31_clouds=np.empty((C31_grid.shape[0], C31_grid.shape[1],))
for i in range(C31_grid.shape[0]):
    for j in range(C31_grid.shape[1]):
        if(C31_grid[i, j] < C1_grid[i, j]*0.3+cirteria_b):
#            print(C31_grid[i, j])
            C31_clouds[i, j]=C31_grid[i, j]
            C1_clouds[i, j]=C1_grid[i, j]

Test if the "criteria" works well

Then we mask and plot C31_clouds to see if our criteria works well

In [0]:
C31_clouds_masked=np.ma.masked_where(np.isnan(C31_clouds), C31_clouds)
fig=plt.figure(figsize=(10.5, 9.5))
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
image=ax.pcolormesh(longitude, latitude, C31_clouds_masked)

Calculate brightness temperature on the top of the clouds

We use planckInvert function as what we have done on satellite3.ipynb to get the brightness temperature at the center of channel-31 (11.02 $\mu m$).

In [0]:
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)
    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

    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]:
cloud_Tbright=planckInvert(11.02, C31_clouds)

Print the maximum and see if it is reasonable.

In [0]:

Save the output as *.mat

We use *.mat as the output data format.

In [0]:
import scipy.io
# save as *.mat
scipy.io.savemat('_share/01_MODIS_L1B_TBright', {'longitude': longitude, 'latitude': latitude, 'cloud_Tbright': cloud_Tbright})

Plot cloud_TBright in mpl_toolkits.basemap

Here we use mpl_toolkits.basemap to Visualize the result.

In [0]:
# mask the result
cloud_Tbright_masked=np.ma.masked_where(np.isnan(cloud_Tbright), cloud_Tbright)

In [0]:
from mpl_toolkits.basemap import Basemap

In [0]:
# Colormap
CMap=plt.cm.hot_r#gist_heat # hot, afmhot, gnuplot
proj=Basemap(resolution='l', projection='lcc', \
            lat_1=30, lat_2=60, lat_0=45, lon_0=-140, \
            llcrnrlon=-155, llcrnrlat=30, \
            urcrnrlon=-110, urcrnrlat=56)
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
## parallels and meridians.
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)
# mask ocean/land to distinguish clouds 
proj.drawlsmask(land_color=[0.925, 0.875, 0.375], ocean_color=[0.375, 0.5, 0.75], \
                lakes=False, resolution='l')
# draw coast 
proj.drawcoastlines(linewidth=1.5, linestyle='solid', color=[0.25, 0.25, 0.25])
# compute native x,y coordinates of grid.
x, y=proj(longitude, latitude)
x_van, y_van=proj(vancity_lon, vancity_lat)
x_text, y_text=proj(vancity_lon+4.5, vancity_lat-0.25)
# pcolor plot
CS=proj.pcolor(x, y, cloud_Tbright_masked, cmap=CMap, vmin=210, vmax=275)
# colorbar
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('Brightness Temperature ( K )', fontsize=12, fontweight='bold')
CBar.ax.tick_params(axis='y', length=0)
# Vancouver
proj.plot(x_van, y_van, marker='o', markersize=18, mfc='k', mec='k')
plt.text(x_text, y_text, 'Vancouver', fontsize=16, fontweight='bold',
                    ha='center', va='center', color='k')
# title
ax.set_title('Brightness Temperature\nMYD021KM.A2006303.2220 channel-31 ',\
             fontweight='bold', fontsize=14)
# Save figure
plt.savefig('_figures/01_MODIS_L1B_TBright.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,
# Show

What we can find on MODIS L1B Image Galary :

In [0]: