Reading MODIS level 1b data in netcdf format

Here is a walkthrough of how to read metadata and data from a level1b file. I've used the the converter I downloaded from this link to produce the Aqua granule netcdf file in this download directory: http://clouds.eos.ubc.ca/~phil/Downloads/a301/MYD021KM.A2005188.0405.005.2009232180906.nc (right click to save to your local drive)

I'm going to use http://unidata.github.io/netcdf4-python/ and two new modules (modismeta.py and netcdflib.py) to read the granule



In [0]:

    
from __future__ import print_function
from __future__ import division
import os,site
currdir=os.getcwd()
head,tail=os.path.split(currdir)
libdir=os.path.join(head,'lib')
site.addsitedir(libdir)
from modismeta import parseMeta
from netcdflib import ncdump
import glob
from netCDF4 import Dataset
from matplotlib import pyplot as plt
import numpy.ma as ma
import numpy as np

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



In [0]:

    
nc_filename=glob.glob('*.nc')
print("found {}".format(nc_filename))



In [0]:

    
nc_file=Dataset(nc_filename[0])

netcdf files have attributes



In [0]:

    
print(nc_file.Earth_Sun_Distance)



In [0]:

    
print(nc_file.HDFEOSVersion)

netcdf files have variables -- stored in a dictionary



In [0]:

    
#print(nc_file.variables.keys())



In [0]:

    
the_long=nc_file.variables['longitude']
the_lat=nc_file.variables['latitude']

the following cell shows how to write code out to a file

%%file snippet.py size=100 lat_data=the_lat[:10,:10] long_data=the_long[:10,:10]

the following cell shows how to read code into a cell

%load snippet.py

Get 500 x 500 slice of lats and lons



In [0]:

    
size=500
lat_data=the_lat[:size,:size]
lon_data=the_long[:size,:size]

Plot some lat/lon pairs



In [0]:

    
%matplotlib inline



In [0]:

    
fig=plt.figure(1,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
the_points=ax1.plot(lon_data,lat_data,'r+')
xlim=(138,150)
ylim=(21,28)
ax1.set_xlim(xlim)
ax1.set_ylim(ylim)



In [0]:

    
lon_bins=np.arange(xlim[0],xlim[1],0.05)
lat_bins=np.arange(ylim[0],ylim[1],0.05)
lon_indices=np.digitize(lon_data.flat,lon_bins)
lat_indices=np.digitize(lat_data.flat,lat_bins)
#print(help(np.arange))
#print(help(lon_data.flat))



In [0]:

    
lon_bins.shape

Working with a partner, make a copy of this notebook and add comments that explain what this code is doing and why. Upload you notebook with this commented cell to connect



In [0]:

    
lat_array=np.zeros([len(lat_bins),len(lon_bins)],dtype=np.float)
bin_count=np.zeros([len(lat_bins),len(lon_bins)],dtype=np.int)
lon_array=np.zeros([len(lat_bins),len(lon_bins)],dtype=np.float)

for point_num in range(len(lat_indices)):
    bin_row=lat_indices[point_num]
    bin_col=lon_indices[point_num]   
    lat_array[bin_row,bin_col]=lat_array[bin_row,bin_col] + lat_data.flat[point_num]
    lon_array[bin_row,bin_col]+=lon_data.flat[point_num]
    bin_count[bin_row,bin_col]+=1
    
for i in range(lon_array.shape[0]):
    for j in range(lon_array.shape[1]):
        if bin_count[i,j] > 0:
            lat_array[i,j]=lat_array[i,j]/bin_count[i,j]
            lon_array[i,j]=lon_array[i,j]/bin_count[i,j]
        else:
            lat_array[i,j]=np.nan
            lon_array[i,j]=np.nan

Plot these arrays to see if they make sense using pcolormesh



In [0]:

    
fig=plt.figure(3,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('counts in each grid cell')
image=ax1.pcolormesh(lon_bins,lat_bins,bin_count)
plt.colorbar(image)

Now get data from Channel 31 (see modis level1b manual page 22.)



In [0]:

    
longwave=nc_file.variables['EV_1KM_Emissive']
longwave.ncattrs()



In [0]:

    
band_names=longwave.band_names
print(band_names)

We want channel 31, which is centered at 11.02 $\mu m$ as describe in this list of MODIS wavelengths



In [0]:

    
band_names=band_names.split(',')
index31=band_names.index('31')



In [0]:

    
index31

Note that the netCDF4 module automatically finds the scale and offset attributes for the longwave netcdf variable and applies them to the 16 bit integer values, converting them to radiances as described in the modis level 1b manual on page 12. See the documentation for set_auto_maskandscale



In [0]:

    
c31_data=longwave[index31,:size,:size]

For Wednesday, use planck_inverse to calculate brightness temperatures a 100 x 100 pixel section of your Channel 31 image and bin the brightness temperatures into a regular 0.05 x 0.05 degree bin



In [0]:

    
c31_array=np.zeros([len(lat_bins),len(lon_bins)],dtype=np.float)

for point_num in range(len(lat_indices)):
    bin_row=lat_indices[point_num]
    bin_col=lon_indices[point_num]   
    c31_array[bin_row,bin_col]=c31_array[bin_row,bin_col] + c31_data.flat[point_num]

for i in range(lon_array.shape[0]):
    for j in range(lon_array.shape[1]):
        if bin_count[i,j] > 0:
            c31_array[i,j]=c31_array[i,j]/bin_count[i,j]
        else:
            c31_array[i,j]=np.nan



In [0]:

    
c31_m= ma.masked_where(np.isnan(c31_array),c31_array)
fig=plt.figure(4,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('channel 31 radiance (W/m^2/micron/sr)')
image=ax1.pcolormesh(lon_bins,lat_bins,c31_m)
plt.colorbar(image)



In [0]:

    
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

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)
    """
    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

Now use our planckInvert function to get the brightness temperature at the center of Channel 31 (11.02 $\mu m$)



In [0]:



In [0]:

    
sat_bright=planckInvert(11.02,c31_array)



In [0]:

    
bright_m= ma.masked_where(np.isnan(sat_bright),sat_bright)
fig=plt.figure(5,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('channel 31 brightness temperatures (K)')
image=ax1.pcolormesh(lon_bins,lat_bins,bright_m)
plt.colorbar(image)

For Friday, make a gridded image of your scene's Channel 1 reflectance and compare it to the thumbnail image that is on the Laadsweb site. Make your grid large enough so you can identify or land features (say 500 km x 500 km). Channel 1 is stored in the EV_250_Aggr1km_RefSB variable.



In [0]:

    
#nc_file.variables.keys()



In [0]:

    
shortwave=nc_file.variables['EV_250_Aggr1km_RefSB']



In [0]:

    
#shortwave.shape
refl_data=shortwave[0,:size,:size]



In [0]:

    
refl_array=np.zeros([len(lat_bins),len(lon_bins)],dtype=np.float)

for point_num in range(len(lat_indices)):
    bin_row=lat_indices[point_num]
    bin_col=lon_indices[point_num]   
    #print(point_num,bin_row)
    refl_array[bin_row,bin_col]=refl_array[bin_row,bin_col] + refl_data.flat[point_num]

for i in range(lon_array.shape[0]):
    for j in range(lon_array.shape[1]):
        if bin_count[i,j] > 0:
            refl_array[i,j]=refl_array[i,j]/bin_count[i,j]
        else:
            refl_array[i,j]=np.nan



In [0]:

    
sat_m= ma.masked_where(np.isnan(refl_array),refl_array)
fig=plt.figure(5,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('channel 1 reflectance')
image=ax1.pcolor(lon_bins,lat_bins,sat_m)
plt.colorbar(image)



In [0]: