#!/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
Data_matrix = np.zeros( (180/resol,360/resol) )
mean_obs_matrix = np.zeros( (180/resol,360/resol) )
print Data_matrix.shape
Data_matrix[:,:] = - 1.0

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

PolderPressureMethod = 'po2' # 'po2' (oxygen) or 'pray' (Rayleigh)

if PolderPressureMethod == 'pray':
	PolerLabel = 'CTP (Rayleigh)'
	PolderFileLabel = 'Rayleigh'
else:
	PolerLabel = 'CTP (Oxygen)'
	PolderFileLabel = 'Oxygen'

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 = []
		pressure = []
		command = ['parasolascii', '-g', polder_granule_name, PolderPressureMethod]
		process = subprocess.Popen(command, stdout=subprocess.PIPE)
		print polder_granule_name
		while True:
			text = process.stdout.readline()
			if text != '':
				polder_lat, polder_long, polder_pressure = text.split()
				lat.append( float(polder_lat) )
				long.append( float(polder_long) )
				pressure.append( float(polder_pressure) )            
			else:
				break
		#print len(lat)
		for d in range(len(lat)):
			index_i = (lat[d] + 90)/resol 
			index_j = (long[d] + 180)/resol 
			if pressure[d] > 0.0 and pressure[d] < 1100:
				mean_obs_matrix[index_i,index_j] = mean_obs_matrix[index_i,index_j] + 1.0
				Data_matrix[index_i,index_j] = Data_matrix[index_i,index_j] + pressure[d]
	print "Day",i ,"done"

MatrixShape = Data_matrix.shape
for i in range(MatrixShape[0]):
	for j in range(MatrixShape[1]):
		if mean_obs_matrix[i,j] != 0:
			Data_matrix[i,j] = Data_matrix[i,j] / mean_obs_matrix[i,j]

np.savetxt('PolderMeanCTP_'+PolderFileLabel+'.txt', Data_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 mean" + PolerLabel + "\n January 2008")
#colors = [(1,1,1)] + [(cm.gist_rainbow(i)) for i in reversed(xrange(1,256))]
colors = [(1,1,1)] + [(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(Data_matrix, cmap=new_map, \
	interpolation='bilinear', vmin=-1, vmax=1100.0 )
bounds=[0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100]
cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, \
	ticks=[0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100])
plt.savefig("PlotPolderMonthlyProductsFromL2_CTP_"+PolderFileLabel+".png", dpi=100 )
plt.show()

About / Au sujet de

Ben
Created:
22 octobre 2016 03:43:38

Updated:
22 octobre 2016 03:43:38

License / Licence

MIT License

Abstract / Résumé