#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from datetime import datetime
from pylab import *
from pyhdf.SD import SD, SDC # To read a HDF file
import numpy as np
import matplotlib.pyplot as plt
import os
import glob
import subprocess

resol = 0.5
coverage_matrix = np.zeros( (180/resol,360/resol) )
mean_obs_matrix = np.zeros( (180/resol,360/resol) )
print coverage_matrix.shape
coverage_matrix[:,:] = - 1.0

year = '2008'
month = '01'
polder_path = '/Volumes/Data/Polder/RB2/'

for i in range(31):
	if i+1 < 10:
		day = '0' + str(i+1)
	else:
		day = str(i+1)
	polder_full_path = polder_path + year + '/' + year + '_' + month + '_' + day + '/' 		
	polder_day_list = glob.glob(polder_full_path + '/*KD' )
	for polder_granule_name in polder_day_list:
		lat = []
		long = []
		cc = []
		command = ['parasolascii', '-g', polder_granule_name, 'cc']
		process = subprocess.Popen(command, stdout=subprocess.PIPE)
		#print polder_granule_name
		while True:
			text = process.stdout.readline()
			if text != '':
				polder_lat, polder_long, polder_cc = text.split()
				lat.append( float(polder_lat) )
				long.append( float(polder_long) )
				cc.append( float(polder_cc) )            
			else:
				break
		#print len(lat)
		for d in range(len(lat)):
			index_i = (lat[d] + 90)/resol 
			index_j = (long[d] + 180)/resol 
			mean_obs_matrix[index_i,index_j] = mean_obs_matrix[index_i,index_j] + 1.0
			coverage_matrix[index_i,index_j] = coverage_matrix[index_i,index_j] + 100.0 * cc[d]
	print "Day",i ,"done"

coverage_matrix[:,:] = coverage_matrix[:,:] / mean_obs_matrix[:,:]
	
m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, \
	llcrnrlon=-180,urcrnrlon=180,resolution='c')
m.drawcoastlines()
m.drawlsmask(land_color='white',ocean_color='white',lakes=True)
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
m.drawmapboundary(fill_color='aqua')
plt.title("POLDER monthly cloud cover \n (defined as the number of cloudy pixels \n divided by the total number of pixels)")
#colors = [(1,1,1)] + [(cm.gist_rainbow(i)) for i in reversed(xrange(1,256))]
colors = [(0,0,0)] + [(cm.jet(i)) for i in xrange(1,256)]
new_map = matplotlib.colors.LinearSegmentedColormap.from_list('new_map', colors, N=256)
#img = m.imshow(A, interpolation='nearest' )
img = m.imshow(coverage_matrix, cmap=new_map, interpolation='bilinear', vmin=-1, vmax=100 )
bounds=[0,10,20,30,40,50,60,70,80,90,100]
cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, ticks=[0,10,20,30,40,50,60,70,80,90,100])
plt.savefig("PlotPolderMonthlyProductsFromL2_CC.png", dpi=100 )
plt.show()

About / Au sujet de

Ben
Created:
22 octobre 2016 03:33:38

Updated:
22 octobre 2016 03:33:38

License / Licence

MIT License

Abstract / Résumé