scale_factor and offset mpl_toolkits.basemap. MYD021KM.A2006303.2220.006.2012078143305.hdf
MYD03.A2006303.2220.006.2012078135515.hdf
MYD021KM.A2006303.2220.006.2012078143305.h5
MYD03.A2006303.2220.006.2012078135515.h5MYD021KM* 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
In [0]:
__author__ = 'ATSC-301 UBC'
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.
Windows user can convert from HDF4 to HDF5 using h4toh5 which provided by hdf group.
One need to do the following steps:
h4h5tools-2.2.2-win32.zip or h4h5tools-2.2.2-win64.zip on local machine.bin folder (e.g. cd F:\HDF_Group\H4TOH5\2.2.2\bin)h4toh5convert input.hdf output.h5
h4toh5convert -na input.hdf output.h5
h4toh5convert -nc4 input.hdf output.h5`
PATH, or set PATH=F:\HDF_Group\H4TOH5\2.2.2\bin and then you can use h4toh5convert everywhere. h4toh5convert in cmd:h4toh5 has a Linux version, so we can use it by repeating the steps above and edit PATH in .bashrc file.
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
Searching all *.h5 files in _data/MODIS_L1B/ directory, Make sure that you have the required HDF5 files exist in it.
In [0]:
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))
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.
Something abou using h5py to discover the structure of HDF files:
a_obj.keys() to get all the attributes in a file object.a_obj['an_attribute'].keys() to see the sub-attribute. dir(a_obj['an_attribute']) to see if the option keys still available.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]:
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))
If you have tools like HDF Explorer or HDF Viewer. You may have a more direct view of HDF's file structure.
Based on chapter 3.B, we read channe-31 and channel-1 data.
'EV_1KM_Emissive' (from channel-20 to 36, 16 different channels)
'MODIS_SWATH_Type_L1B' '/' 'Data Fields' '/' 'EV_1KM_Emissive''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.
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, :, :]
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()
print(a_list)
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]:
C31_scale=a_list[6][10]
C31_offset=a_list[7][10]
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]:
np.max(C1)
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('===================================================')
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)))
Here we define function reproj_L1B to reproject data to a regular grid.
One can also use MODIS Reprojection Tool to do the work.
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)
-----------------------------------------------------------------------------------------
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
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)
longitude.shape
Plot the result
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)
#plt.colorbar(image)
plt.show
Be careful to chose res here, high resolution makes your computer slow. Also, it affects the result on 2-D histgram part.
In [0]:
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)
Histogram is some kind of a basic image segmentation technique. Here we apply 2-D histogram to distinguish clouds from ocean.
**`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 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=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()
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=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()
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
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]
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=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
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
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
In [0]:
cloud_Tbright=planckInvert(11.02, C31_clouds)
Print the maximum and see if it is reasonable.
In [0]:
np.nanmax(cloud_Tbright)
np.nanmin(cloud_Tbright)
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})
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]:
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()
What we can find on MODIS L1B Image Galary :
In [0]: