01_MODIS_L1B.ipynb to MODIS L2 (06) Cloud Product. 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.h5MYD06* 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.
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:
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})