Objectives

  • The objectives of this notebook is to apply the steps we did in 01_MODIS_L1B.ipynb to MODIS L2 (06) Cloud Product.
  • The files required for this notebook including:
      MYD06_L2.A2006303.2220.051.2009072115006.hdf
      MYD03.A2006303.2220.006.2012078135515.hdf
      MYD06_L2.A2006303.2220.051.2009072115006.h5
      MYD03.A2006303.2220.006.2012078135515.h5
  • MYD06* are MODIS AQUA L2 Cloud Product, observed on Oct 30 2006. The temporal and spatial info. is the same as what's in 01_MODIS_L1B.ipynb.
  • MYD03* is the geolocation file.
  • *.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'

The following cells actually repeats the steps in 01_MODIS_L1B.ipynb, so there is not extra documentation for them.

Here we use (and save) 5 different variables in MYD06 product:

  • Cloud Fraction, resolution is 1km×1km
  • Cirrus Reflectance, 250m×250m
  • Cloud Top Pressure (hPa), 1km×1km
  • Cloud Top Temperature (K), 1km×1km
  • Cloud Optical Depth, 250m×250m

Geographical information for 250m×250m datasets is in MYD03 product, so we still need that file in the directory of _data/MODIS_L1B.


In [0]:
import h5py
import glob
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
from __future__ import print_function
from mpl_toolkits.basemap import Basemap
% matplotlib inline

In [0]:
hdf5_06=glob.glob('_data/MODIS_L2/MYD06*.h5')
print("MODIS MYD06 file found {}".format(hdf5_06))
hdf5_07=glob.glob('_data/MODIS_L2/MOD07*.h5')
print("MODIS MOD07 file found {}".format(hdf5_07))
hdf5_Geo=glob.glob('_data/MODIS_L1B/MYD03*.h5')
print("MODIS Geolocation file found {}".format(hdf5_Geo))

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

In [0]:
lon1km= mod06_obj['mod06']['Geolocation Fields']['Longitude'][:]
lat1km= mod06_obj['mod06']['Geolocation Fields']['Latitude'][:]
lonqkm=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][:]
latqkm=geo_obj['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][:]

In [0]:
CFrac = mod06_obj['mod06']['Data Fields']['Cloud_Fraction'][:].astype(float)
Ci_ref = mod06_obj['mod06']['Data Fields']['Cirrus_Reflectance'][:].astype(float)
CTP = mod06_obj['mod06']['Data Fields']['Cloud_Top_Pressure'][:].astype(float)
CTT = mod06_obj['mod06']['Data Fields']['Cloud_Top_Temperature'][:].astype(float)
COT = mod06_obj['mod06']['Data Fields']['Cloud_Optical_Thickness'][:].astype(float)

In [0]:
FillVal = mod06_obj['mod06']['Data Fields']['Cloud_Fraction'].attrs.values()[2]
scale = mod06_obj['mod06']['Data Fields']['Cloud_Fraction'].attrs.values()[5]
offset = mod06_obj['mod06']['Data Fields']['Cloud_Fraction'].attrs.values()[6]
CFrac[CFrac == FillVal] = np.nan
CFrac = (CFrac - offset*np.ones(CFrac.shape))*scale

FillVal = mod06_obj['mod06']['Data Fields']['Cirrus_Reflectance'].attrs.values()[3]
scale = mod06_obj['mod06']['Data Fields']['Cirrus_Reflectance'].attrs.values()[6]
offset = mod06_obj['mod06']['Data Fields']['Cirrus_Reflectance'].attrs.values()[7]
Ci_ref[Ci_ref == FillVal] = np.nan
Ci_ref = (Ci_ref - offset*np.ones(Ci_ref.shape))*scale

FillVal = mod06_obj['mod06']['Data Fields']['Cloud_Top_Pressure'].attrs.values()[2]
scale = mod06_obj['mod06']['Data Fields']['Cloud_Top_Pressure'].attrs.values()[5]
offset = mod06_obj['mod06']['Data Fields']['Cloud_Top_Pressure'].attrs.values()[6]
CTP[CTP == FillVal] = np.nan
CTP = (CTP - offset*np.ones(CTP.shape))*scale

FillVal = mod06_obj['mod06']['Data Fields']['Cloud_Top_Temperature'].attrs.values()[2]
scale = mod06_obj['mod06']['Data Fields']['Cloud_Top_Temperature'].attrs.values()[5]
offset = mod06_obj['mod06']['Data Fields']['Cloud_Top_Temperature'].attrs.values()[6]
CTT[CTT == FillVal] = np.nan
CTT = (CTT - offset*np.ones(CTT.shape))*scale

FillVal = mod06_obj['mod06']['Data Fields']['Cloud_Optical_Thickness'].attrs.values()[2]
scale = mod06_obj['mod06']['Data Fields']['Cloud_Optical_Thickness'].attrs.values()[5]
offset = mod06_obj['mod06']['Data Fields']['Cloud_Optical_Thickness'].attrs.values()[6]
COT[COT == FillVal] = np.nan
COT = (COT - offset*np.ones(COT.shape))*scale

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)
    -----------------------------------------------------------------------------------------
    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, [resolution_x, resolution_y].
    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[0])
    y_bins=np.arange(ylim[0], ylim[1], res[1])
#    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

In [0]:
xlim=[np.min(lon1km), np.max(lon1km)]
ylim=[np.min(lat1km), np.max(lat1km)]
CFrac_grid, longitude1km, latitude1km, _ = reproj_L1B(CFrac, lon1km, lat1km, xlim, ylim, [0.33, 0.33])
CTP_grid, _, _, _ = reproj_L1B(CTP, lon1km, lat1km, xlim, ylim, [0.33, 0.33])
CTT_grid, _, _, _ = reproj_L1B(CTT, lon1km, lat1km, xlim, ylim, [0.33, 0.33])
Ci_ref_grid, longitudeqkm, latitudeqkm, _ = reproj_L1B(Ci_ref, lonqkm, latqkm, xlim, ylim, [0.05, 0.05])
COT_grid, _, _, _ = reproj_L1B(COT, lonqkm, latqkm, xlim, ylim, [0.05, 0.05])

In [0]:
longitude1km=np.ma.masked_where(np.isnan(longitude1km), longitude1km)
latitude1km=np.ma.masked_where(np.isnan(latitude1km), latitude1km)
CFrac_mask=np.ma.masked_where(np.isnan(CFrac_grid), CFrac_grid)
Ci_ref_mask=np.ma.masked_where(np.isnan(Ci_ref_grid), Ci_ref_grid)
CTP_mask=np.ma.masked_where(np.isnan(CTP_grid), CTP_grid)
CTT_mask=np.ma.masked_where(np.isnan(CTT_grid), CTT_grid)
COT_mask=np.ma.masked_where(np.isnan(COT_grid), COT_grid)

In [0]:
def plot_MODIS(lon, lat, z, zlim, CMap, name_var, name_title, name_to_save):
    vancity_lat=49.25
    vancity_lon=-123.1
    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)
    fig=plt.figure(figsize=(12, 12))
    ax=plt.gca()
    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)
#    proj.drawlsmask(land_color=[0.925, 0.875, 0.375], ocean_color=[0.375, 0.5, 0.75], \
#                    lakes=False, resolution='l')
    proj.drawcoastlines(linewidth=1.5, linestyle='solid', color=[0.25, 0.25, 0.25])
    x, y=proj(lon, lat)
    x_van, y_van=proj(vancity_lon, vancity_lat)
    x_text, y_text=proj(vancity_lon+4.5, vancity_lat-0.25)
    CS=proj.pcolor(x, y, z, cmap=CMap, vmin=zlim[0], vmax=zlim[1])
    CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
    CBar.set_label(name_var, fontsize=12, fontweight='bold')
    CBar.ax.tick_params(axis='y', length=0)
    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')
    ax.set_title(name_title, fontweight='bold', fontsize=14)
    path='_figures/'
    plt.savefig(path+name_to_save, dpi=450, facecolor='w', edgecolor='w',
                orientation='portrait', papertype='a4', format='png',
                transparent=True, bbox_inches='tight', pad_inches=0,
                frameon=None)
    plt.show()

In [0]:
#hdf5_06[0][-43:-21]
plot_MODIS(longitude1km, latitude1km, CFrac_mask*100, [0, 100], plt.cm.jet_r, \
           'Cloud Fraction ( % )', 'Cloud Fraction\n'+hdf5_06[0][-43:-21], '01_MODIS_L2_CFrac.png')

In [0]:
plot_MODIS(longitude1km, latitude1km, CTP_mask, [50, 1050], plt.cm.jet, \
           'Cloud Top Pressure ( hPa )', 'Cloud Top Pressure\n'+hdf5_06[0][-43:-21], '01_MODIS_L2_CTP.png')

In [0]:
plot_MODIS(longitude1km, latitude1km, CTT_mask, [210, 290], plt.cm.hot_r, \
           'Cloud Top Temperature ( K )', 'Cloud Top Temperature\n'+hdf5_06[0][-43:-21], '01_MODIS_L2_CTT.png')

In [0]:
plot_MODIS(longitudeqkm, latitudeqkm, Ci_ref_mask, [0, 0.5], plt.cm.gray, \
           'Cirrus Reflectance', 'Cirrus Reflectance\n'+hdf5_06[0][-43:-21], '01_MODIS_L2_CiRef.png')

In [0]:
plot_MODIS(longitudeqkm, latitudeqkm, COT_mask, [0, 100], plt.cm.YlGnBu_r, \
           'Cloud Optical Thickness', 'Cloud Optical Thickness\n'+hdf5_06[0][-43:-21], '01_MODIS_L2_COT.png')

In [0]:
scipy.io.savemat('_share/01_MODIS_L2_CFrac', {'longitude': longitude1km, 'latitude': latitude1km, 'CFrac': CFrac_grid})
scipy.io.savemat('_share/01_MODIS_L2_CTP', {'longitude': longitude1km, 'latitude': latitude1km, 'CTP': CTP_grid})
scipy.io.savemat('_share/01_MODIS_L2_CTT', {'longitude': longitude1km, 'latitude': latitude1km, 'CTT': CTT_grid})
scipy.io.savemat('_share/01_MODIS_L2_CiRef', {'longitude': longitudeqkm, 'latitude': latitudeqkm, 'CiRef': Ci_ref_grid})
scipy.io.savemat('_share/01_MODIS_L2_COT', {'longitude': longitudeqkm, 'latitude': latitudeqkm, 'COT': COT_grid})