#!/usr/bin/env python

from pyhdf.SD import SD, SDC # To read a HDF file (SD data)
from pyhdf.HDF import *
from pyhdf.VS import *

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

from mpl_toolkits.axes_grid1 import make_axes_locatable
import pylab

#---------- bits_stripping Function ----------#

def bits_stripping(bit_start,bit_count,value):
	bitmask=pow(2,bit_start+bit_count)-1
	return np.right_shift(np.bitwise_and(value,bitmask),bit_start)

#---------- Read HDF Files (VD data) ----------#

f = HDF('2008001002019_08922_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf') 
vs = f.vstart() 
Latitude = vs.attach('Latitude')
Longitude = vs.attach('Longitude')

a = Latitude[:]

Latitude.detach() # "close" the vdata
Longitude.detach() # "close" the vdata
vs.end() # terminate the vdata interface
f.close() 

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

cpr_2b_geoprof = SD('2008001002019_08922_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf', SDC.READ)

cpr_2b_geoprof_height = cpr_2b_geoprof.select('Height')
cpr_2b_geoprof_cloud_mask = cpr_2b_geoprof.select('CPR_Cloud_mask')
cpr_2b_geoprof_radar_reflectivity = cpr_2b_geoprof.select('Radar_Reflectivity')

cpr_2b_geoprof_height_data = cpr_2b_geoprof_height.get()
cpr_2b_geoprof_cloud_mask_data = cpr_2b_geoprof_cloud_mask.get()
cpr_2b_geoprof_radar_reflectivity_data = cpr_2b_geoprof_radar_reflectivity.get()

fillvalue = 15360
missing = -8888

img = np.zeros((800,125))
img.fill(-30)

for i in np.arange(800):
for j in np.arange(125):
	k = int( 125 * (cpr_2b_geoprof_height_data[i,j] + 5000) / 35000 )
	if cpr_2b_geoprof_radar_reflectivity_data[i,j] > -3000 and \
	cpr_2b_geoprof_radar_reflectivity_data[i,j] < 2100: 		
		img[i,k] = cpr_2b_geoprof_radar_reflectivity_data[i,j] / 100

ax = plt.subplot(111)
im = ax.imshow(img.T, interpolation='bilinear', cmap=cm.jet, origin='lower', extent=[0,200,-10,60])

plt.title( 'CLOUDSAT Radar Reflectivity' )

plt.ylabel('Height (km)',fontsize=14)
pylab.yticks([-10,0,10,20,30,40,50,60], [-5,0,5,10,15,20,25,30])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.10)
plt.colorbar(im, cax=cax)

plt.savefig("cloudsat_radar_reflectivity.png")
plt.show()

About / Au sujet de

Ben
Created:
22 octobre 2016 03:09:51

Updated:
22 octobre 2016 03:09:51

License / Licence

MIT License

Abstract / Résumé