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)


"""
Explanation: Reading the Lidar LayerTops variable
1) First, grab the lidar file
End of explanation
"""


the_file=glob.glob('2006*LIDAR*h5')[0]
print(the_file)

"""
Explanation: 2) use glob.glob wildcards to read the filename
from disk without having to get the name exactly right
End of explanation
"""


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))

"""
Explanation: 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
End of explanation
"""


lat,lon,date_times,prof_times,dem_elevation=get_geo(the_file)

time_minutes = prof_times/60.

"""
Explanation: 4) get the time values (in seconds) for the orbit and convert to decimal minutes for plotting
End of explanation
"""


hit = layer_tops < 10    
#http://cswww.cira.colostate.edu/dataSpecs.php    
layer_tops = (layer_tops - offset)/factor
layer_tops[hit] = np.nan



"""
Explanation: 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.
End of explanation
"""


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))

"""
Explanation: 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
End of explanation
"""


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))

#%matplotlib inline
from IPython.display import display
plt.close('all')
plt.ioff()


def plot_layers(layer_tops,ax):
    ntimes,nlayers=layer_tops.shape
    for i in range(nlayers):
        label='layer {}'.format(i)
        ax.plot(layer_tops[:,i],label=label)
    ax.legend()
    return ax

fig, ax = plt.subplots(1,1)
ax=plot_layers(layer_tops,ax)
#display(fig);
#help(display)
print(plt.isinteractive())


ax.set(title = 'new title')
display(fig)

ax.set(title='new title')
display(fig)

"""
Explanation: 7) What fraction of the time is there a layer 2 cloud above layer 1 cloud?
End of explanation
"""


