#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from scipy.misc import bytescale
from mpl_toolkits.basemap import Basemap, cm 
from scipy import interpolate
from scipy.interpolate import griddata
from matplotlib import colors

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

# References:
# https://earthdata.nasa.gov/sites/default/files/field/document/MODIS_True_Color.pdf
# http://www.idlcoyote.com/ip_tips/brightmodis.html

#----------------------------------------------------------------------------------------#
# inputs

file_name_myd021km = 'MYD021KM.A2013189.1350.006.2013190155358.hdf'
file_name_myd03 = 'MYD03.A2013189.1350.006.2013190152754.hdf'

file_myd021km = SD(file_name_myd021km, SDC.READ)

#----------------------------------------------------------------------------------------#

selected_sds = file_myd021km.select('EV_250_Aggr1km_RefSB')

selected_sds_attributes = selected_sds.attributes()

for key, value in selected_sds_attributes.iteritems():
	#print key, value
	if key == 'reflectance_scales':
		reflectance_scales_250_Aggr1km_RefSB = np.asarray(value)		
	if key == 'reflectance_offsets':
		reflectance_offsets_250_Aggr1km_RefSB = np.asarray(value)	

sds_data_250_Aggr1km_RefSB = selected_sds.get()

data_shape = sds_data_250_Aggr1km_RefSB.shape

along_track = data_shape[1]
cross_trak = data_shape[2]

#----------------------------------------------------------------------------------------#

selected_sds = file_myd021km.select('EV_500_Aggr1km_RefSB')

selected_sds_attributes = selected_sds.attributes()

for key, value in selected_sds_attributes.iteritems():
	if key == 'reflectance_scales':
		reflectance_scales_500_Aggr1km_RefSB = np.asarray(value)	
	if key == 'reflectance_offsets':
		reflectance_offsets_500_Aggr1km_RefSB = np.asarray(value)	

sds_data_500_Aggr1km_RefSB = selected_sds.get()

#----------------------------------------------------------------------------------------#

file_myd03 = SD(file_name_myd03, SDC.READ)

selected_sds = file_myd03.select('Latitude')
myd03_lat = selected_sds.get()

selected_sds = file_myd03.select('Longitude')
myd03_long = selected_sds.get()

#----------------------------------------------------------------------------------------#

data_shape = sds_data_250_Aggr1km_RefSB.shape

along_track = data_shape[1]
cross_trak = data_shape[2]

z = np.zeros((along_track, cross_trak,3))

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,0] = ( sds_data_250_Aggr1km_RefSB[0,i,j] - \
        reflectance_offsets_250_Aggr1km_RefSB[0] ) * \
        reflectance_scales_250_Aggr1km_RefSB[0] 

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,1] = ( sds_data_500_Aggr1km_RefSB[1,i,j] - \
        reflectance_offsets_500_Aggr1km_RefSB[1] ) * \
        reflectance_scales_500_Aggr1km_RefSB[1]  

for i in np.arange(along_track):
    for j in np.arange(cross_trak): 
        z[i,j,2] = ( sds_data_500_Aggr1km_RefSB[0,i,j] - \
        reflectance_offsets_500_Aggr1km_RefSB[0] ) * \
        reflectance_scales_500_Aggr1km_RefSB[0] 

z[ z > 1 ] = 1.0
z[ z < 0 ] = 0.0

#----------------------------------------------------------------------------------------#
# Rough estimation of latitude and longitude at granule center (long_0, lat_0)

lat_min = myd03_lat[0,0]
lat_max = myd03_lat[along_track-1,cross_trak-1]

lat_0 = lat_min + (lat_max - lat_min) / 2.

long_min = min(myd03_long[0,0],myd03_long[along_track-1,cross_trak-1])
long_max = max(myd03_long[0,0],myd03_long[along_track-1,cross_trak-1])

lon_0 = long_min + (long_max - long_min) / 2.

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

fig = plt.figure()

ax = fig.add_subplot(111)

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) 

xpt1, ypt1 = m1(myd03_long[0,0],myd03_lat[0,0]) 
xpt2, ypt2 = m1(myd03_long[0,cross_trak-1],myd03_lat[0,cross_trak-1]) 
xpt3, ypt3 = m1(myd03_long[along_track-1,cross_trak-1], \
                myd03_lat[along_track-1,cross_trak-1])
xpt4, ypt4 = m1(myd03_long[along_track-1,0],myd03_lat[along_track-1,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)

x_igrid, y_igrid = m(myd03_long,myd03_lat)

x_igrid = x_igrid - xpt0
y_igrid = y_igrid - ypt0

z_igrid_01 = np.zeros((along_track, cross_trak))
z_igrid_02 = np.zeros((along_track, cross_trak))
z_igrid_03 = np.zeros((along_track, cross_trak))

for i in np.arange(2030):
    for j in np.arange(1354): 
        z_igrid_01[i,j] = z[i,j,0]
        z_igrid_02[i,j] = z[i,j,1]
        z_igrid_03[i,j] = z[i,j,2]

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()
z_igrid_01 = z_igrid_01.ravel()
z_igrid_02 = z_igrid_02.ravel()
z_igrid_03 = z_igrid_03.ravel()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T
xi, yi = np.mgrid[llx:urx:1000j, lly:ury:1000j]

z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='cubic')
z_02 = griddata(xy1_igrid, z_igrid_02, (xi, yi), method='cubic')
z_03 = griddata(xy1_igrid, z_igrid_03, (xi, yi), method='cubic')

rgb_projected = np.zeros((1000, 1000,3))
for i in np.arange(1000):
    for j in np.arange(1000): 
        rgb_projected[i,j,0] = z_01[i,j]
        rgb_projected[i,j,1] = z_02[i,j]
        rgb_projected[i,j,2] = z_03[i,j]

#rgb_projected[ z > 1 ] = 1.0
#rgb_projected[ z < 0 ] = 0.0
whereAreNaNs = np.isnan(rgb_projected);
rgb_projected[whereAreNaNs] = 0.75;

img = m.imshow(np.rot90(np.fliplr(rgb_projected)), origin='lower')

m.drawcoastlines()

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

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

plt.title('MODIS RGB Orthographic Projection', fontsize=10)

plt.savefig("modis_granule_rgb_orth_projection.png", dpi=100)
plt.show()

#----------------------------------------------------------------------------------------#

About / Au sujet de

Ben
Created:
20 juin 2016 21:50:42

Updated:
20 juin 2016 21:50:42

License / Licence

MIT License

Abstract / Résumé