#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm 
from pyhdf.SD import SD, SDC 
from scipy import interpolate
from scipy.interpolate import griddata

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

from matplotlib import colors

import matplotlib.cm as cm

# fix nearest interpolation create mask points too far
THRESHOLD = 5000.0

from scipy.interpolate.interpnd import _ndim_coords_from_arrays
from scipy.spatial import cKDTree

#----------------------------------------------------------------------------------------#
# Inputs

myd03_file_name = './inputs/myd06_granules_viz/MYD03.A2008015.1435.006.2012066174347.hdf'
myd06_file_name = './inputs/myd06_granules_viz/MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files

myd03 = SD(myd03_file_name, SDC.READ)
myd06 = SD(myd06_file_name, SDC.READ)

myd03_Latitude = myd03.select('Latitude')
myd03_Longitude = myd03.select('Longitude')
data_selected_id = myd06.select('Cloud_Phase_Optical_Properties') 

myd03_Latitude_data = myd03_Latitude.get()
myd03_Longitude_data = myd03_Longitude.get()

#----------------------------------------------------------------------------------------#
# Find lat_0 and lon_0 of the MODIS granule

#print myd03_Latitude_data.shape
#print myd03_Longitude_data.shape

ArrayShape = myd03_Latitude_data.shape

#print ArrayShape[0], ArrayShape[1]

as0 = ArrayShape[0] - 1
as1 = ArrayShape[1] - 1

latmin = myd03_Latitude_data[0,0]
latmax = myd03_Latitude_data[as0,as1]

lat_0 = latmin + (latmax - latmin) / 2.

#print "lat_0:", lat_0
	
#----------------------------------------------------------------------------------------#

tmp_01 = myd03_Longitude_data[0,0]
tmp_02 = myd03_Longitude_data[as0,as1]

lonmin = min(myd03_Longitude_data[0,0],myd03_Longitude_data[as0,as1])
lonmax = max(myd03_Longitude_data[0,0],myd03_Longitude_data[as0,as1])

lon_0 = lonmin + (lonmax - lonmin) / 2.

if lon_0 > 180:
	lon_0 = - ( 360 - lon_0 )

#print "lon_0:", lon_0

myd03_Longitude_data[0,0] = tmp_01
myd03_Longitude_data[as0,as1] = tmp_02

#print latmin, latmax
#print lonmin, lonmax
#print "lat_0 and lon_0 done"

#----------------------------------------------------------------------------------------#
# Orthographic Map

fig = plt.figure(figsize=(6, 6), dpi=100)

ax = fig.add_subplot(111)
#ax.patch.set_facecolor('black')
ax.patch.set_facecolor((0.75,0.75,0.75))

m1 = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution=None)

xpt0, ypt0 = m1(lon_0,lat_0) 

#print 'xpt0, ypt0',xpt0, ypt0

xpt1, ypt1 = m1(myd03_Longitude_data[0,0],myd03_Latitude_data[0,0]) 
xpt2, ypt2 = m1(myd03_Longitude_data[0,as1],myd03_Latitude_data[0,as1]) 
xpt3, ypt3 = m1(myd03_Longitude_data[as0,as1],myd03_Latitude_data[as0,as1])
xpt4, ypt4 = m1(myd03_Longitude_data[as0,0],myd03_Latitude_data[as0,0])

llx = min(xpt1,xpt2,xpt3,xpt4) - xpt0  # lower left
lly = min(ypt1,ypt2,ypt3,ypt4) - ypt0

urx = max(xpt1,xpt2,xpt3,xpt4) - xpt0  # upper right
ury = max(ypt1,ypt2,ypt3,ypt4) - ypt0

m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution='l',\
    llcrnrx=llx,llcrnry=lly,urcrnrx=urx,urcrnry=ury)

#print "Orthographic map done"

#----------------------------------------------------------------------------------------#
# Plot MODIS data: Using Scipy Griddata

long = myd03_Longitude.get()
lat = myd03_Latitude.get()

x_igrid, y_igrid = m(long,lat)

x_igrid = x_igrid - xpt0
y_igrid = y_igrid - ypt0

z_igrid = data_selected_id.get()

attributes = data_selected_id.attributes()

z_igrid = z_igrid * attributes['scale_factor']

z_igrid_shape = z_igrid.shape

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()
z1_igrid = z_igrid.ravel()

#print "z1_igrid.min(),z1_igrid.max()", z1_igrid.min(),z1_igrid.max()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T

xi, yi = np.mgrid[llx:urx:1000j, lly:ury:1000j]

z = griddata(xy1_igrid, z1_igrid, (xi, yi), method='nearest')

#----------------------------------------------------------------------------------------#
# Construct kd-tree, functionality copied from scipy.interpolate

tree = cKDTree(xy1_igrid)
arr_x = _ndim_coords_from_arrays((xi, yi))
dists, indexes = tree.query(arr_x)

z = z[:]
z[dists > THRESHOLD] = np.nan

cmap =  [(0.0,0.0,0.0)] + [(32./255., 107./255., 164./255.)] + \
       [(187./255., 217./255., 238./255.)] + [(255./255, 51./255, 51./255.)] 
cmap = colors.ListedColormap(cmap)

bounds = [-0.5, 1.5, 2.5, 3.5, 4.5]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(z.T, origin='lower', cmap=cmap, norm=norm )

cbar_bounds = [1.5, 2.5, 3.5, 4.5]
cbar_ticks = [2.0,3.0,4.0]  
cbar_labels = ['Liquid','Ice','Und.']  
		   
cbar = m.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

m.drawcoastlines()

#m.fillcontinents(color='coral',lake_color='white')
#m.drawmapboundary(fill_color='white')  

m.drawparallels(np.arange(-90.,120.,5.), color='k', labels=[True,False,False,False])
m.drawmeridians(np.arange(0.,420.,10.), color='k', labels=[False,False,False,True])

ax.set_xlabel("", fontsize=10)
ax.set_ylabel("", fontsize=10)       

plt.title('Cloud Phase Optical Properties \n MYD06 C6 (2008-01-15; 14h35)', fontsize=10)

plt.savefig("myd06_c6_cloud_phase_optical_properties_orthographic_projection.png", dpi=200)
plt.show()
plt.close()

About / Au sujet de

Ben
Created:
3 avril 2017 17:33:43

Updated:
3 avril 2017 17:33:43

License / Licence

MIT License

Abstract / Résumé