import numpy as N
import netCDF4 as S
import matplotlib.pyplot as plt
import mpl_toolkits.basemap as bm
import math as M

direcInput_a = "/Users/Shared/DATA/bedekkingsgraad_nubiscope-2011/" # set data directory
direcInput_b = "/Users/Shared/DATA/CESAR-surfdata-2011/" # set data directory

direcOutput = "/Users/Shared/DATA/CloudCoverResearch/" # set directory containing plots of results


g = 9.81
pi = M.pi
nt_len_a = 4442
nt_len_b = 4464
day_len = 31

#open file
a = S.Dataset(direcInput_a+"cesar_nubiscope_cloudcover_la1_t10_v1.0_201101.nc", mode='r') # Nubiscope data (Cabauw)
b = S.Dataset(direcInput_b+"cesar_surface_meteo_lc1_t10_v1.0_201101.nc", mode='r') # Near surface measurements at Cabauw


t_a = a.variables["time"][:]   #hours
t_b = b.variables["time"][:]   #hours

cldcover_total = a.variables["cldcover_total"][:]
cldcover_low = a.variables["cldcover_low"][:]
T_2 = b.variables["TA002"][:]
Td_2 = b.variables["TD002"][:]

day_a = N.zeros((nt_len_a), dtype='d')
day_b = N.zeros((nt_len_b), dtype='d')
for nt in range(nt_len_a):
 day_a[nt] = t_a[nt] / 24
for nt in range(nt_len_b):
 day_b[nt] = t_b[nt] / 24

 #print t[nt], day[nt], TotalCloudCover[nt]
 
plt.figure(figsize=(12,10))
plt.plot(day_a, cldcover_total,label='Total cloud cover', color='black')
plt.axis([0,31,-10,110]) 
plt.legend()   
plt.xlabel('Day')
plt.ylabel('Cloud cover fraction [%]')
plt.grid(True)
plt.savefig(direcOutput+"TotalCloudCover_Jan2011"+".png") # save plot as png-file 
plt.show()

plt.figure(figsize=(12,10))
plt.plot(day_a, cldcover_low,label='Low cloud cover', color='blue')
plt.axis([0,31,-10,110]) 
plt.legend()   
plt.xlabel('Day')
plt.ylabel('Cloud cover fraction [%]')
plt.grid(True)
plt.savefig(direcOutput+"LowCloudCover_Jan2011"+".png") # save plot as png-file 
plt.show()

plt.figure(figsize=(12,10))
plt.plot(day_b, T_2,label='temperature', color='red')
plt.plot(day_b, Td_2,label='dewpoint temperature', color='blue')
plt.axis([0,31,-10,30]) 
plt.legend()   
plt.xlabel('Day')
plt.ylabel('Temperature (red) and dewpoint (blue) [C]')
plt.grid(True)
plt.savefig(direcOutput+"TempTd_Jan2011"+".png") # save plot as png-file 
plt.show()