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
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.ncattrs())
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
In [0]:
size=100
lat_data=the_lat[:size,:size]
lon_data=the_long[:size,:size]
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+')
ylim=(22.4,23.8)
xlim=(146.,150.)
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)
In [0]:
lon_bins.shape
In [0]:
lon_indices
lat_indices
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
In [0]:
import numpy.ma as ma
lat_m= ma.masked_where(np.isnan(lat_array),lat_array)
help(ma.masked_where)
In [0]:
from matplotlib import pyplot as plt
fig=plt.figure(1,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('binned latitude')
image=ax1.pcolormesh(lon_bins,lat_bins,lat_m)
plt.colorbar(image)
In [0]:
lon_m= ma.masked_where(np.isnan(lon_array),lon_array)
fig=plt.figure(2,figsize=(9,9))
plt.clf()
ax1=fig.add_subplot(1,1,1)
ax1.set_title('binned longitude')
image=ax1.pcolormesh(lon_bins,lat_bins,lon_m)
plt.colorbar(image)
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)
In [0]:
longwave=nc_file.variables['EV_1KM_Emissive']
longwave.ncattrs()
In [0]:
band_names=longwave.band_names
In [0]:
band_names=band_names.split(',')
index31=band_names.index('31')
In [0]:
index31
In [0]:
sat_data=longwave[index31,:size,:size]
In [0]:
sat_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]
sat_array[bin_row,bin_col]=sat_array[bin_row,bin_col] + sat_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:
sat_array[i,j]=sat_array[i,j]/bin_count[i,j]
else:
sat_array[i,j]=np.nan
In [0]:
sat_m= ma.masked_where(np.isnan(sat_array),sat_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,sat_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
In [0]:
sat_bright=planckInvert(11.02,sat_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)
In [0]: