#!/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()

Abstract / Résumé

Announcement

Bienvenue sur "The Open Science Notebook Project" ! Le site est en développement, vous pouvez nous contacter par le biais du formulaire de contact si vous avez des questions ou si vous avez simplement découvert une erreur !

Date de la prochaine mise à jour le: 27 janvier 2018

About

Created: April 3, 2017, 3:33 p.m.

Last Edit: April 3, 2017, 3:33 p.m.

MIT License