## Plotting brightness temperatures in a Lambert Conformal Conic map projection

In this notebook we're going to continue working with http://clouds.eos.ubc.ca/~phil/Downloads/a301/MYD021KM.A2005188.0405.005.2009232180906.h5. We will also need to go to Laadsweb and do a wildcard search on the filename:

MYD03.A2005188.0405.*

which will produce the geometry file which can be downloaded and converted to hdf5 to get

The MYD03 file contains the lat/lon for every pixel and is described at http://modaps.nascom.nasa.gov/services/about/products/MYD03.html



In [0]:

from __future__ import print_function
import os,site
import glob
import h5py
from IPython.display import Image
import numpy as np
from matplotlib import pyplot as plt
#
# add the lib folder to the path assuming it is on the same
# level as the notebooks folder
#
libdir=os.path.abspath('../lib')
import h5dump



the glob function finds a file using a wildcard to save typing (google: python glob wildcard)



In [0]:

h5_filename=glob.glob('../data/MYD02*.h5')
print("found {}".format(h5_filename))




In [0]:

h5_file=h5py.File(h5_filename[0])



note that channel 31 occurs at index value 10



In [0]:

index31=10



the data is stored as unsigned, 2 byte integers which can hold values from 0 to $2^{16}$ - 1 = 65,535



In [0]:

chan31=h5_file['MODIS_SWATH_Type_L1B']['Data Fields']['EV_1KM_Emissive'][index31,:,:]
print(chan31.shape,chan31.dtype)




In [0]:

chan31[:3,:3]



we need to apply a scale and offset to convert to radiance (the netcdf module did this for us automatically

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

this information is included in the attributes of each variable.

(see page 36 of the Modis users guide )

here is the scale for all 16 channels



In [0]:

print(scale)



and here is the offset for 16 channels



In [0]:

print(offset)




In [0]:

chan31=(chan31 - offset[index31])*scale[index31]



now convert this to brightness temperature



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]:

chan31_Tbright=planckInvert(11.02, chan31)




In [0]:

%matplotlib inline



histogram the calibrated radiances and show that they lie between 0-10 $W\,m^{-2}\,\mu m^{-1}\,sr^{-1}$



In [0]:

import matplotlib.pyplot as plt
out=plt.hist(chan31_Tbright.flat)



note that the longitude and latitude arrays are (406,271) while the actual data are (2030,1354). These lat/lon arrays show only every fifth row and column. We need to get the full lat/lon arrays from the MYD03 file



In [0]:

geom_filename=glob.glob('../data/MYD03*.h5')
print("found {}".format(h5_filename))




In [0]:

geom_h5=h5py.File(geom_filename[0])




In [0]:

h5dump.dumph5(geom_h5)




In [0]:

the_long=geom_h5['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][...]
the_lat=geom_h5['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][...]




In [0]:

print(the_long.shape,the_lat.shape)




In [0]:

print('===================================================')
print('Size of Longitude: {}'.format(the_long.shape))
print('Longitude Range: {} ~ {}'.format(np.min(the_long), np.max(the_long)))
print('===================================================')
print('Size of Latitude: {}'.format(the_lat.shape))
print('Latitude Range: {} ~ {}'.format(np.min(the_lat), np.max(the_lat)))




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.
=========================================================================================
'''

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



now regrid the radiances and brightness temperatures on a 0.1 x 0.1 degree regular lat/lon grid



In [0]:

xlim=[np.min(the_long), np.max(the_long)]
ylim=[np.min(the_lat), np.max(the_lat)]
chan31_grid, longitude, latitude, bin_count = reproj_L1B(chan31, the_long, the_lat, xlim, ylim, 0.1)
tbright_grid,longitude,latitude,bin_count=reproj_L1B(chan31_Tbright, the_long, the_lat, xlim, ylim, 0.1)




In [0]:

longitude.shape



Plot this gridded data without a map projections



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, chan31_grid)



Now replot using an lcc (Lambert conformal conic) projection from basemap at http://matplotlib.org/basemap/users/examples.html



In [0]:

from mpl_toolkits.basemap import Basemap




In [0]:

lcc_values=dict(resolution='l',projection='lcc',
lat_1=20,lat_2=40,lat_0=30,lon_0=135,
llcrnrlon=120,llcrnrlat=20,
urcrnrlon=150,urcrnrlat=42)




In [0]:

proj=Basemap(**lcc_values)




In [0]:

fig=plt.figure(figsize=(12, 12))
## define parallels and meridians to draw.
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)
# draw coast & fill continents
#map.fillcontinents(color=[0.25, 0.25, 0.25], lake_color=None) # coral
out=proj.drawcoastlines(linewidth=1.5, linestyle='solid', color='k')
x, y=proj(longitude, latitude)
# contourf the bathmetry
CS=proj.pcolor(x, y, chan31_grid, cmap=plt.cm.hot)
# colorbar
CBar.set_label('Channel 31 radiance ($W\,m^{-2}\,\mu m\,sr^{-1})$', fontsize=10)
CBar.ax.tick_params(axis='y', length=0)



repeat for brightness temperature



In [0]:

fig=plt.figure(figsize=(12, 12))
## define parallels and meridians to draw.
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)
# draw coast & fill continents
#map.fillcontinents(color=[0.25, 0.25, 0.25], lake_color=None) # coral
out=proj.drawcoastlines(linewidth=1.5, linestyle='solid', color='k')
x, y=proj(longitude, latitude)
# contourf the bathmetry
CS=proj.pcolor(x, y, tbright_grid, cmap=plt.cm.hot)
# colorbar
CBar.set_label('Channel 31 Brightness temperature (K)', fontsize=10)
CBar.ax.tick_params(axis='y', length=0)




In [0]: