In [0]:

    




    

In [0]:

    




    

In [0]:

    




    

In [0]:

    




    

In [0]:

    




    

In [0]:

    

import glob
import h5py
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
from __future__ import print_function
% matplotlib inline


    

In [0]:

    

hdf5_filename=glob.glob('_data/AIRS_L2_IR/*.h5')
print("HDF5 file found {}".format(hdf5_filename))


    

In [0]:

    




    

In [0]:

    

print('{}\n--------------------------------------------------------------------'.format(hdf5_filename[0]))
test_obj=h5py.File(hdf5_filename[0], 'r')
test_lon=test_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Longitude'][:]
test_lat=test_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Latitude'][:]
print('Longitude Range: \n{} ~ {}'.format(np.min(test_lon), np.max(test_lon)))
print('Latitude Range: \n{} ~ {}'.format(np.min(test_lat), np.max(test_lat)))
print('====================================================================')
print('{}\n--------------------------------------------------------------------'.format(hdf5_filename[1]))
test_obj=h5py.File(hdf5_filename[1], 'r')
test_lon=test_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Longitude'][:]
test_lat=test_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Latitude'][:]
print('Longitude Range: \n{} ~ {}'.format(np.min(test_lon), np.max(test_lon)))
print('Latitude Range: \n{} ~ {}'.format(np.min(test_lat), np.max(test_lat)))


    

In [0]:

    

lat_min=9999
lat_max=-9999
lon_min=9999
lon_max=-9999
for i in range(len(hdf5_filename)):
    hdf_obj=h5py.File(hdf5_filename[i], 'r')
    temp_lat=hdf_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Latitude'][:]
    temp_lon=hdf_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Longitude'][:]
    if(temp_lat.min() < lat_min):
        lat_min=temp_lat.min()
    if(temp_lat.max() > lat_max):
        lat_max=temp_lat.max()
    if(temp_lon.min() < lon_min):
        lon_min=temp_lon.min()
    if(temp_lon.max() > lon_max):
        lon_max=temp_lon.max()

lat_lim=[lat_min, lat_max]
lon_lim=[lon_min, lon_max]
print('Overall Latitude Range: {} ~ {}'.format(lat_lim[0], lat_lim[1]))
print('Overall Longitude Range: {} ~ {}'.format(lon_lim[0], lon_lim[1]))


    

In [0]:

    




    

In [0]:

    

def muti_proj_init(lon_lim, lat_lim, res):
    lon_bins=np.arange(lon_lim[0], lon_lim[1], res)
    lat_bins=np.arange(lat_lim[0], lat_lim[1], res)

    overall_lon=np.zeros([len(lat_bins), len(lon_bins)], dtype=np.float)
    overall_lat=np.zeros([len(lat_bins), len(lon_bins)], dtype=np.float)
    overall_data=np.zeros([len(lat_bins), len(lon_bins)], dtype=np.float)
    overall_bincount=np.zeros([len(lat_bins), len(lon_bins)], dtype=np.int)
    
    return lon_bins, lat_bins, overall_lon, overall_lat, overall_data, overall_bincount


    

In [0]:

    




    

In [0]:

    




    

In [0]:

    

def multi_proj_process(raw_data, raw_x, raw_y, lon_bins, lat_bins, overall_data, overall_lon, overall_lat, overall_bincount):
    
    '''
    =========================================================================================
    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)
    -----------------------------------------------------------------------------------------
    =========================================================================================
    '''
    import numpy as np

    x_indices=np.searchsorted(lon_bins, raw_x.flat, 'right')
    y_indices=np.searchsorted(lat_bins, raw_y.flat, 'right')
            
    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
        if((bin_row > 0) and (bin_col > 0)):
            overall_bincount[bin_row, bin_col] += 1
            overall_lon[bin_row, bin_col] += raw_x.flat[n]
            overall_lat[bin_row, bin_col] += raw_y.flat[n]
            overall_data[bin_row, bin_col] += raw_data.flat[n]

                
    return overall_lon, overall_lat, overall_data, overall_bincount


    

In [0]:

    




    

In [0]:

    

def multi_proj_end(overall_lon, overall_lat, overall_data, overall_bincount):                   
    for i in range(overall_lon.shape[0]):
        for j in range(overall_lon.shape[1]):
            if overall_bincount[i, j] > 0:
                overall_lon[i, j]=overall_lon[i, j]/overall_bincount[i, j]
                overall_lat[i, j]=overall_lat[i, j]/overall_bincount[i, j]
                overall_data[i, j]=overall_data[i, j]/overall_bincount[i, j]
            else:
                overall_lon[i, j]=np.nan
                overall_lat[i, j]=np.nan
                overall_data[i, j]=np.nan
                
    return overall_lon, overall_lat, overall_data


    

In [0]:

    

lat_lim=[-90, 90]
lon_lim=[-180, 180]


    

In [0]:

    

lon_bins, lat_bins, overall_lon, overall_lat, overall_data, overall_bincount = muti_proj_init(lon_lim, lat_lim, 1)
for i in range(len(hdf5_filename)):
    hdf_obj=h5py.File(hdf5_filename[i], 'r')
    raw_y=hdf_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Latitude'][:]
    raw_x=hdf_obj['L2_Standard_atmospheric&surface_product']['Geolocation Fields']['Longitude'][:]
    raw_CldFrc=hdf_obj['L2_Standard_atmospheric&surface_product']['Data Fields']['CldFrcTot'][:]
    raw_CldFrc[raw_CldFrc == -9999]=np.nan
    overall_lon, overall_lat, overall_data, overall_bincount = \
        multi_proj_process(raw_CldFrc, raw_x, raw_y, lon_bins, lat_bins, overall_data, overall_lon, overall_lat, overall_bincount)
    
longitude, latitude, CldFrc = multi_proj_end(overall_lon, overall_lat, overall_data, overall_bincount)


    

In [0]:

    

CldFrc_masked=np.ma.masked_where(np.isnan(CldFrc), CldFrc)
longitude=np.ma.masked_where(np.isnan(longitude), longitude)
latitude=np.ma.masked_where(np.isnan(latitude), latitude)
longitude.shape


    

In [0]:

    

fig=plt.figure(figsize=(12, 8))
#clev=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
ax=plt.gca()
ax.set_xlim(lon_lim[0], lon_lim[1])
ax.set_ylim(lat_lim[0], lat_lim[1])
image=ax.pcolormesh(longitude, latitude, CldFrc_masked, cmap=plt.cm.gray_r)
plt.colorbar(image)
plt.show


    

In [0]:

    

from mpl_toolkits.basemap import Basemap

lonlim=lon_lim
latlim=lat_lim
vancity_lat=49.25
vancity_lon=-123.1

proj = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,\
            llcrnrlon=-180,urcrnrlon=180,resolution='c')
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
ax=plt.gca()
## parallels and meridians.
parallels=np.arange(-90, 90, 30)
meridians=np.arange(0, 360, 60)
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
proj.drawcoastlines(linewidth=1.0, linestyle='solid', color='k')
# 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+15, vancity_lat)
# pcolor plot
CS=proj.pcolor(x, y, CldFrc_masked, cmap=plt.cm.gray_r, vmin=0, vmax=1)
# colorbar
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('CldFrc', fontsize=10)
CBar.ax.tick_params(axis='y', length=0)
#CBar.ax.invert_yaxis() 
# Vancouver
proj.plot(x_van, y_van, marker='o', markersize=12, mfc=[0.3, 0.6, 0.8], mec='k')
#plt.text(x_text, y_text, 'Vancouver', fontsize=16, fontweight='bold',
#                    ha='center', va='center', color=[0.15, 0.3, 0.4])
# title
ax.set_title('AIRS L2 CldFrc 240-file-combined',\
             fontweight='bold', fontsize=14)
# Save figure
plt.savefig('test.png', dpi=400, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,
            frameon=None)
# Print
plt.show()


    

In [0]:

    




    

In [0]:

    

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)
# draw coast & fill continents
#map.fillcontinents(color=[0.25, 0.25, 0.25], lake_color=None) # coral
proj.drawcoastlines(linewidth=1.0, linestyle='solid', color='k')
# 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, CldFrc_masked, cmap=plt.cm.gray_r, vmin=0, vmax=1)
#CS=proj.contourf(x, y, CldFrc_masked, cmap=plt.cm.gray_r)
# colorbar
CBar=proj.colorbar(CS, 'right', size='5%', pad='5%')
CBar.set_label('CldFrc', fontsize=10)
CBar.ax.tick_params(axis='y', length=0)
#CBar.ax.invert_yaxis() 
# Vancouver
proj.plot(x_van, y_van, marker='o', markersize=18, mfc=[0.3, 0.6, 0.8], mec='k')
plt.text(x_text, y_text, 'Vancouver', fontsize=16, fontweight='bold',
                    ha='center', va='center', color=[0.15, 0.3, 0.4])
# title
ax.set_title('AIRS L2 CldFrc Vancouver',\
             fontweight='bold', fontsize=14)
# Save figure
#plt.savefig('test2.png', dpi=400, facecolor='w', edgecolor='w',
#            orientation='portrait', papertype='a4', format='png',
#            transparent=True, bbox_inches='tight', pad_inches=0,
#            frameon=None)
# Print
plt.show()


    

In [0]:

    




    

In [0]:

    




    

In [0]:

    




    

In [0]: