In [1]:
import glob
import h5py
import numpy as np
import glob
from a301lib.cloudsat import get_geo
from a301utils.a301_readfile import download
from matplotlib import pyplot as plt
lidar_name='2006303212128_02702_CS_2B-GEOPROF-LIDAR_GRANULE_P2_R04_E02.h5'
download(lidar_name)
2) use glob.glob wildcards to read the filename from disk without having to get the name exactly right
In [2]:
the_file=glob.glob('2006*LIDAR*h5')[0]
print(the_file)
3) Use hdfview to figure out the path to the LayerTop variable, and to get the factor and offset needed to turn the 16 bit interger data into science values as described on the cloudsat web page
In [3]:
with h5py.File(the_file,'r') as in_file:
layer_tops=in_file['2B-GEOPROF-LIDAR']['Data Fields']['LayerTop'][...]
factor=in_file['2B-GEOPROF-LIDAR']['Data Fields']['LayerTop'].attrs['factor']
offset=in_file['2B-GEOPROF-LIDAR']['Data Fields']['LayerTop'].attrs['offset']
units=in_file['2B-GEOPROF-LIDAR']['Data Fields']['LayerTop'].attrs['units']
missing = in_file['2B-GEOPROF-LIDAR']['Data Fields']['LayerTop'].attrs['missing']
#
# the next line turns the numpy bytes (b'm') object returned by h5py into a unicode string
# for printing
#
units=units.decode('utf-8')
print('missing value, factor, offset and units: {} {} {} {}'
.format(missing,factor,offset,units))
4) get the time values (in seconds) for the orbit and convert to decimal minutes for plotting
In [4]:
lat,lon,date_times,prof_times,dem_elevation=get_geo(the_file)
time_minutes = prof_times/60.
5) turn the missing values (-99) into np.nan ("not a number") so they will be dropped from our plot. Count any cloud height below 10 meters as noise and assign it as missing.
In [5]:
hit = layer_tops < 10
#http://cswww.cira.colostate.edu/dataSpecs.php
layer_tops = (layer_tops - offset)/factor
layer_tops[hit] = np.nan
6) go through each of the layers and find the mean height (not counting the nan values) and the number of timesteps where there was cloud detected
In [6]:
num_times,num_layers = layer_tops.shape
text="""
layer number: {0:}
cloud fraction is {1:4.2f}%
mean height is {2:6.1f} meters
"""
for the_layer in range(num_layers):
missing = np.isnan(layer_tops[:,the_layer])
present = np.logical_not(missing)
num_present = np.sum(present)
percent_present=num_present/num_times*100.
mean_height = np.nanmean(layer_tops[:,the_layer])
print(text.format(the_layer,percent_present,mean_height))
7) What fraction of the time is there a layer 2 cloud above layer 1 cloud?
In [7]:
layer1= layer_tops[:,0]
layer2= layer_tops[:,1]
layer1_count=0
overlap_count=0
for index,layer1_height in enumerate(layer1):
if not np.isnan(layer1_height):
layer1_count += 1
if not np.isnan(layer2[index]):
overlap_count += 1
overlap_freq=100.*overlap_count/layer1_count
print(('when there was cloud in layer 1, there was also cloud in layer2 {:6.2f} '
'percent of the time')
.format(overlap_freq))
In [29]:
%matplotlib inline
from IPython.display import display
plt.close('all')
meters2km=1.e3
seconds2mins=60.
def plot_layers(time_secs,layer_tops,ax):
ntimes,nlayers=layer_tops.shape
time_mins=time_secs/seconds2mins
for i in range(nlayers):
label='layer {}'.format(i)
ax.plot(time_mins,layer_tops[:,i]/meters2km,label=label)
ax.legend()
return ax
fig, ax = plt.subplots(1,1,figsize=(12,4))
ax=plot_layers(prof_times,layer_tops,ax)
ax.set(title='Cloudsat Orbit -- lidar/radar cloud tops',
xlabel='time (minutes)',ylabel='height (km)');
#
# expand to view the 60-70 minute time inteval
#
hit=np.logical_and(prof_times > 60*60,prof_times < 70*60)
fig, ax = plt.subplots(1,1,figsize=(12,3))
ax=plot_layers(prof_times[hit],layer_tops[hit,:],ax)
ax.set(title='Cloudsat Orbit -- lidar/radar cloud tops -- zoomed',
xlabel='time (minutes)',ylabel='height (km)');
In [ ]: