<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Objectives
The objectives of this notebook is to provide examples about working with AQUA MODIS L1B files, including
Step1:  Content
Convert HDF4 format file to HDF5 format
Convert HDF4 to HDF5 on Windows
Convert HDF4 to HDF5 on Linux
Import modules & tools
Read MODIS L1B data through h5py
Get HDF5 file object
The structure of HDF file
Get attribute and sub-attribute information
Some external tools
Read MODIS L1B channel-31 and channel-1 data
Read raw data
Select the channel
Scale factor and offset value
Geolocation field
Reproject MODIS L1B data to a regular grid
Define function reproj_L1B
Test if reproj_L1B works well
Convert channel-31 and channel-1 data
Identify cloud and ocean via 2-D histogram
How to use numpy.histogram2d
Create 2-D Histgram for channel-31 and channel-1 data
Eliminate ocean points based on 2-D histgram
Test if the "criteria" works well
Calculate brightness temperature on the top of the clouds
Save the brightness temperature result as *.mat file
Plot cloud_TBright in mpl_toolkits.basemap
Convert HDF4 format file to HDF5 format
The hdf4 module in Python, <a href='http
Step2:  We import h5py to read HDF5 files
Step3:  scipy.io for saving data in *.mat format
Step4:  For the map view of data, we need mpl_toolkits.basemap
Step5:  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.
Step6:  Get HDF5 file object
Step7:  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
Step8:  Some external tools
If you have tools like <a href='http
Step9:  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
Step10:  Then we can chose the channel we want
Step11:  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.
Step12:  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.
Step13:  radiance_scales and radiance_offsets are the 7th and 8th group of a_list, and channel-31 is the 11th element of the group.
Step14:  We do the same thing for channel-1 data, but now we use reflectance_scales
Step15:  Finally, we correct the data, numpy.ones is the same as ones.m in MATLAB.
Step16:  Print the maximum data and see if it is reasonable.
Step17:  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
Step18:  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 <a href='http
Step19:  Test if reproj_L1B works well
Here we reproject channel-31 data to see if reproj_L1B works well.
Step20:  Mask NaN for plot, also make sure that the data is not too big to plot.
Step21:  Plot the result
Step22:  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 <a href='#Create-2-D-Histgram-for-channel-31-and-channel-1-data'>2-D histgram</a> part.
Step23:  Identify cloud and ocean via 2-D histogram
Histogram is some kind of a basic <a href='http
Step24:  Here, we define channel-1 data on x-axis and call np.histogram2d as what's in <a href='#How-to-use-numpy.histogram2d'>above section</a> to get bin_count value x_edges and y_edges. Noted that masked NumPy array has no attribute flat.
Step25:  Then we make 2-D histgram to see the difference between clouds and ocean, the core idea is
Step26:  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
Step27:  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
Step28:  Test if the "criteria" works well
Then we mask and plot C31_clouds to see if our criteria works well
Step30:  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$).
Step31:  Print the maximum and see if it is reasonable.
Step32:  Save the output as *.mat
We use *.mat as the output data format.
Step33:  Plot cloud_TBright in mpl_toolkits.basemap
Here we use mpl_toolkits.basemap to Visualize the result.
<ASSISTANT_TASK:>
Python Code:
__author__ = 'ATSC-301 UBC'
Explanation: Objectives
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.A2006303.2220.006.2012078143305.hdf
        MYD03.A2006303.2220.006.2012078135515.hdf
        MYD021KM.A2006303.2220.006.2012078143305.h5
        MYD03.A2006303.2220.006.2012078135515.h5
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 <a href='http://ladsweb.nascom.nasa.gov/data/search.html'>LAADS Web</a>.
*.h5 is the HDF5 file which can be <a href='#Convert-HDF4-format-file-to-HDF5-format'>converted</a> from HDF4 files.
End of explanation
import glob
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
from __future__ import print_function
% matplotlib inline
Explanation: Content
Convert HDF4 format file to HDF5 format
Convert HDF4 to HDF5 on Windows
Convert HDF4 to HDF5 on Linux
Import modules & tools
Read MODIS L1B data through h5py
Get HDF5 file object
The structure of HDF file
Get attribute and sub-attribute information
Some external tools
Read MODIS L1B channel-31 and channel-1 data
Read raw data
Select the channel
Scale factor and offset value
Geolocation field
Reproject MODIS L1B data to a regular grid
Define function reproj_L1B
Test if reproj_L1B works well
Convert channel-31 and channel-1 data
Identify cloud and ocean via 2-D histogram
How to use numpy.histogram2d
Create 2-D Histgram for channel-31 and channel-1 data
Eliminate ocean points based on 2-D histgram
Test if the "criteria" works well
Calculate brightness temperature on the top of the clouds
Save the brightness temperature result as *.mat file
Plot cloud_TBright in mpl_toolkits.basemap
Convert HDF4 format file to HDF5 format
The hdf4 module in Python, <a href='http://pysclint.sourceforge.net/pyhdf'>pyhdf</a> 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 <a href='https://github.com/fhs/python-hdf4'>here</a>. Linux user can try to install python module <a href='https://conda.binstar.org/mutirri'>mutirri</a> , which includes pyhdf in it.
Convert HDF4 to HDF5 on Windows
Windows user can convert from HDF4 to HDF5 using <a href='http://www.hdfgroup.org/h4toh5/download.html'>h4toh5</a> 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:
bash
        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:
<img src="_figures/01_h4toh5convert.PNG" width=550 align="left"/>
Convert HDF4 to HDF5 on Linux
<a href='http://www.hdfgroup.org/h4toh5/download.html'>h4toh5</a> 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:
bash
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 <a href='http://docs.h5py.org/en/latest/index.html'>documentation</a> and their <a href='http://www.h5py.org/'>website</a>. Also, one can find more information about HDF5's file format from <a href='http://www.hdfgroup.org/HDF5/doc/H5.format.html'>here</a>.
End of explanation
import h5py
Explanation: We import h5py to read HDF5 files:
End of explanation
import scipy.io
Explanation: scipy.io for saving data in *.mat format
End of explanation
from mpl_toolkits.basemap import Basemap
Explanation: For the map view of data, we need mpl_toolkits.basemap
End of explanation
hdf5_L1B=glob.glob('_data/MODIS_L1B/MYD021*.h5')
print("MODIS L1B file found {}".format(hdf5_L1B))
hdf5_Geo=glob.glob('_data/MODIS_L1B/MYD03*.h5')
print("MODIS Geolocation file found {}".format(hdf5_Geo))
Explanation: 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.
End of explanation
hdf5_obj=h5py.File(hdf5_L1B[0], 'r')
geo_obj=h5py.File(hdf5_Geo[0], 'r')
Explanation: Get HDF5 file object
End of explanation
hdf5_attr=hdf5_obj.keys()
print('Attributes in {}'.format(hdf5_L1B))
print('=============================================================')
print('hdf5_attr=hdf5_obj.keys() \n\n{}'.format(hdf5_attr))
sub_attr=hdf5_obj['MODIS_SWATH_Type_L1B'].keys()
print("\n\n\tSub-attributes in 'MODIS_SWATH_Type_L1B'")
print('\t=============================================================')
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\t=============================================================')
print("\t\tsubsub_attr=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields'].keys() \n\n\t\t{}".format(subsub_attr))
Explanation: 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:
End of explanation
# 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))
Explanation: Some external tools
If you have tools like <a href='http://hdfeos.org/software/hdfexplorer.php'>HDF Explorer</a> or <a href='http://www.hdfgroup.org/products/java/hdfview/'>HDF Viewer</a>. You may have a more direct view of HDF's file structure.
<img src="_figures/01_HDF_file_structure.PNG" width=850 align="left"/>
Read MODIS L1B channel-31 and channel-1 data
Read raw data
Based on chapter <a href='#Structure-of-HDF-file'>3.B</a>, 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'
End of explanation
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))
Explanation: 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'
End of explanation
C31=L1B_emiss[10, :, :]
C1=L1B_ref[1, :, :]
Explanation: Then we can chose the channel we want:
End of explanation
print('Channel-31 info\n===============================================================')
hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs.items()
Explanation: 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.
End of explanation
a_list=hdf5_obj['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'].attrs.values()
print(a_list)
Explanation: 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.
End of explanation
C31_scale=a_list[6][10]
C31_offset=a_list[7][10]
Explanation: radiance_scales and radiance_offsets are the 7th and 8th group of a_list, and channel-31 is the 11th element of the group.
End of explanation
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]
Explanation: We do the same thing for channel-1 data, but now we use reflectance_scales
End of explanation
C31=(C31 - C31_offset * np.ones(C31.shape))*C31_scale
C1=(C1 - C1_offset * np.ones(C1.shape))*C1_scale
Explanation: Finally, we correct the data, numpy.ones is the same as ones.m in MATLAB.
End of explanation
np.max(C1)
Explanation: Print the maximum data and see if it is reasonable.
End of explanation
C_x=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][:]
C_y=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][:]
print('===================================================')
print('Size of Longitude: {}'.format(C_x.shape))
print('Longitude Interval: {} ~ {}'.format(np.min(C_x), np.max(C_x)))
print('===================================================')
print('Size of Latitude: {}'.format(C_y.shape))
print('Latitude Interval: {} ~ {}'.format(np.min(C_y), np.max(C_y)))
Explanation: 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'
End of explanation
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)
    -----------------------------------------------------------------------------------------
    Input:
            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.
    Output:
            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.
    Note:
            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_col=x_indices[n]-1
        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] 
            else:
                d_array[i, j]=np.nan
                x_array[i, j]=np.nan
                y_array[i,j]=np.nan
                
    return d_array, x_array, y_array, bin_count
Explanation: 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 <a href='http://www.whrc.org/education/rwanda/pdf/Baccini_MRT.pdf'>MODIS Reprojection Tool</a> 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
python
numpy.searchsorted(b, a, 'right')
to replace
python
numpy.digitize(a, b)
Because <a href='http://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html'>numpy.searchsorted</a> use binary search but <a href='http://docs.scipy.org/doc/numpy/reference/generated/numpy.digitize.html'>numpy.digitize</a> based on linear search, the later could be slow when we have lots of data.
End of explanation
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)
Explanation: Test if reproj_L1B works well
Here we reproject channel-31 data to see if reproj_L1B works well.
End of explanation
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
Explanation: Mask NaN for plot, also make sure that the data is not too big to plot.
End of explanation
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)
#plt.colorbar(image)
plt.show
Explanation: Plot the result
End of explanation
res=0.05;
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)
Explanation: 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 <a href='#Create-2-D-Histgram-for-channel-31-and-channel-1-data'>2-D histgram</a> part.
End of explanation
# 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)
Explanation: Identify cloud and ocean via 2-D histogram
Histogram is some kind of a basic <a href='http://en.wikipedia.org/wiki/Image_segmentation#Histogram-based_methods'>image segmentation</a> technique. Here we apply 2-D histogram to distinguish clouds from ocean.
How to use numpy.histogram2d
<a href='http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram2d.html'>numpy.histogram2d</a> 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:
python
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 
<a href='http://docs.scipy.org/doc/numpy/reference/generated/numpy.digitize.html'>numpy.digitize</a> 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 <a href='#Define-function-reproj_L1B'>reproj_L1B</a>).
There is a counter part of numpy.histogram2d for 1-D histogram named 
<a href='http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html'>numpy.histogram2d</a>. One can also make histograms through <a href='http://matplotlib.org/api/pyplot_api.html'>pyplot.hist</a> and <a href='http://matplotlib.org/api/pyplot_api.html'>pyplot.hist2d</a>.
MATLAB users can use function <a href='http://www.mathworks.com/help/matlab/ref/histcounts.html'>histcounts.m</a> and 
<a href='http://www.mathworks.com/help/matlab/ref/histogram.html'>histogram.m</a> 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. 
<a href='http://www.mathworks.com/matlabcentral/fileexchange/45325-efficient-2d-histogram--no-toolboxes-needed'>here</a>).
Create 2-D Histgram for channel-31 and channel-1 data
Create bins for channel-31 and channel-1
End of explanation
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])
Explanation: Here, we define channel-1 data on x-axis and call np.histogram2d as what's in <a href='#How-to-use-numpy.histogram2d'>above section</a> to get bin_count value x_edges and y_edges. Noted that masked NumPy array has no attribute flat.
End of explanation
# 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=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))
# ========== Main ========== #
# 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', 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)
divider=make_axes_locatable(axMain)
# 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=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 ========== #
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))
[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,
            frameon=None)
# show
plt.show()
Explanation: Then we make 2-D histgram to see the difference between clouds and ocean, the core idea is:
```python
 # 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, x5.5+6.5np.ones(x.shape))
``` 
The rest are codes for figures and axises.
End of explanation
fig=plt.figure(figsize=(8, 8))
ax=plt.gca()
plt.grid()
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,
            frameon=None)
plt.show()
Explanation: 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.
End of explanation
criteria_k=5.5 # less than
cirteria_b=6.5
C1_clouds=np.empty((C31_grid.shape[0], C31_grid.shape[1],))
C1_clouds[:]=np.nan
C31_clouds=np.empty((C31_grid.shape[0], C31_grid.shape[1],))
C31_clouds[:]=np.nan
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]
Explanation: 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:
python
id=C1*5.5+6.5; C31[C31&lt;id]=np.nan
End of explanation
C31_clouds_masked=np.ma.masked_where(np.isnan(C31_clouds), C31_clouds)
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_clouds_masked)
#plt.colorbar(image)
plt.show
Explanation: Test if the "criteria" works well
Then we mask and plot C31_clouds to see if our criteria works well
End of explanation
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
    c1=2.*h*c**2.
    c2=h*c/kb
    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
cloud_Tbright=planckInvert(11.02, C31_clouds)
Explanation: 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$).
End of explanation
np.nanmax(cloud_Tbright)
np.nanmin(cloud_Tbright)
Explanation: Print the maximum and see if it is reasonable.
End of explanation
import scipy.io
# save as *.mat
scipy.io.savemat('_share/01_MODIS_L1B_TBright', {'longitude': longitude, 'latitude': latitude, 'cloud_Tbright': cloud_Tbright})
Explanation: Save the output as *.mat
We use *.mat as the output data format.
End of explanation
# mask the result
cloud_Tbright_masked=np.ma.masked_where(np.isnan(cloud_Tbright), cloud_Tbright)
from mpl_toolkits.basemap import Basemap
lonlim=xlim
latlim=ylim
vancity_lat=49.25
vancity_lon=-123.1
# 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))
ax=plt.gca()
## 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)
#CBar.ax.invert_yaxis() 
# 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,
            frameon=None)
# Show
plt.show()
Explanation: Plot cloud_TBright in mpl_toolkits.basemap
Here we use mpl_toolkits.basemap to Visualize the result.
End of explanation
<END_TASK>