from pylab import *
import netCDF4 as NC
import glob
import scipy.ndimage.filters as F
import datetime
import matplotlib.dates as mdates
import matplotlib.pyplot as pyplot

files = sorted(glob.glob('*.nc'))
D = NC.Dataset(files[0])
hhl = D.variables['range'][:]/1e3
D.close()

for f in files:
    D = NC.Dataset(f)
    beta_raw = D.variables['beta_raw']
    stddev = D.variables['stddev']
    rangebin = D.variables['range']

    beginningoftime = datetime.datetime(1904, 1, 1, 0, 0)
    time = [beginningoftime + datetime.timedelta(seconds=t) for t in D.variables['time'][:]]

    signal = (beta_raw[:].T * stddev[:]).T
    rcsignal = signal * rangebin[:]**2
    if 'rcs' in locals():
        rcs = vstack([rcs, rcsignal])
        times = hstack([times, time])
    else:
        rcs = rcsignal
        times = time
    D.close() 

width = 5000
frames = 3000
start_indices = np.arange(0, len(rcs)-width, len(rcs)/frames)

figure(figsize=(16,4),dpi=2000)
for i, istart in enumerate(start_indices):
    print 'Plotting', i, istart
    iend = istart+width

    tstart = times[istart]
    tend = times[iend]
    x_lims = mdates.date2num([tstart, tend])
    extent = [x_lims[0], x_lims[1],  hhl[10], hhl[-1]]

    clf()
    
    data = rcs[istart:iend, 20:]
    data = F.uniform_filter(data, (10, 10))
    data = np.log10(data)

    plt.xticks(rotation=00)

    imshow(data.T, origin='lower', cmap='nipy_spectral', vmin=3, vmax=7, extent=extent, aspect='auto')
    colorbar()
    xlabel('local time, Munich [hrs]')
    ylabel('height [km]')
    title('YALIS Ceilometer backscatter signal')

    gca().xaxis_date()
    gca().get_xaxis().set_major_formatter(mdates.DateFormatter('%H'))
    locator = mdates.AutoDateLocator(maxticks=6, interval_multiples=True)
    locator.intervald[HOURLY] = [3]
    gca().get_xaxis().set_major_locator(locator)

    savefig('lidar{0:04d}.png'.format(i), bbox_inches='tight')