#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

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

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

file_name = 'MYD04_L2.A2015190.0010.006.2015190191844.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

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

file = SD(file_name, SDC.READ)

longitude_id = file.select('Longitude')
latitude_id = file.select('Latitude')
data_selected_id = file.select('Optical_Depth_Land_And_Ocean')

longitude_attributes = longitude_id.attributes()
latitude_attributes = latitude_id.attributes()
attributes = data_selected_id.attributes()

longitude = longitude_id.get()
latitude = latitude_id.get()
data = data_selected_id.get()

_FillValue = attributes['_FillValue']

#----------------------------------------------------------------------------------------#
# Find longitude and latitude indexes

long_Halemaumau = -155.28
lat_Halemaumau = 19.41

X = np.sqrt( np.square( longitude - long_Halemaumau ) + \
             np.square( latitude - lat_Halemaumau ) )
LongLat_idx = np.where( X == X.min() )

#----------------------------------------------------------------------------------------#
# Plot data

data[ data == _FillValue ] = 0.0

data = data * attributes['scale_factor']

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0,0.0001,0.1,0.2,0.3,0.4,0.5]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')                  

plt.scatter(135-LongLat_idx[1],LongLat_idx[0],color='r',s=100)
plt.annotate('Halemaumau Crater',xy=(135-LongLat_idx[1]-75,LongLat_idx[0]+10),\
color='red',bbox=dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.80))
                         
plt.xticks([0,25,50,75,100,135], 
           ['135','100','50','75','25','0'], rotation=0, fontsize=8 )
plt.yticks([0,50,100,150,203], 
           ['0','50','100','150','203'], 
           rotation=0, fontsize=8 )

title = 'Optical Depth Land And Ocean \n MYD04 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

cbar_bounds = [0.0,0.1,0.2,0.3,0.4,0.5]
cbar_ticks = [0.0,0.1,0.2,0.3,0.4,0.5] 

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=8)

plt.xlim(0.0,135)
plt.ylim(0.0,203)

plt.savefig("myd04_c6_Optical_Depth_Land_And_Ocean.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2017-04-27

Edited / Modifié: 2017-04-27

License / Licence: MIT License

Abstract / Résumé