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 = "C:/Users/Mine/Documents/3.5/DATA/bedekkingsgraad_nubiscope-2011/" # set data directory
#direcInput_b = "C:/Users/Mine/Documents/3.5/DATA/CESAR-surfdata-2011/" # set data directory

#direcOutput = "C:/Users/Mine/Documents/3.5/Results/" # set directory containing plots of results

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"][:]
RH = (100*(M.e**((17.625*Td_2)/(243.04+Td_2))/M.e**((17.625*T_2)/(243.04+T_2))))

#------------------------------------------NEW
dRH = 2
binmax = 100/dRH
number = N.zeros((binmax), dtype='d')
x_axis = N.zeros((binmax), dtype='d')
RHmax=0

for nt in range(1,nt_len_a): 
  if RH[nt] > 100:
   RH[nt] = 100.0

for bin in range(binmax):
 RHmax = RHmax+dRH
 x_axis[bin] = RHmax-(0.5*dRH)
 number[bin] = 0  
 for nt in range(1,nt_len_a): 
  if (RH[nt]>RHmax-dRH and RH[nt]<=RHmax):
   number[bin] = number[bin]+1

RHmax=0
print 'bin  RH   frequency'
for bin in range(binmax):
 RHmax = RHmax+dRH
 print bin, x_axis[bin], number[bin]/nt_len_a 

plt.figure(figsize=(12,10))
plt.plot(x_axis, 100*number/nt_len_a ,label='frequency', color='black')
plt.axis([0,100,0,100]) 
plt.legend()   
plt.xlabel('Relative humidity [%]')
plt.ylabel('frequency [%] (bin= 2%)')
plt.grid(True)
plt.savefig(direcOutput+"Test_frequency"+".png") # save plot as png-file 
plt.show()

#------------------------------------------NEW

print min(RH)


######################################
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
 
plt.figure(figsize=(12,10))
plt.plot(day_b, RH,label='Humidity', color='blue')
plt.axis([0,31,55,110]) 
plt.legend()   
plt.xlabel('Day')
plt.ylabel('Relative Humidity [%]')
plt.grid(True)
plt.savefig(direcOutput+"Test_RH"+".png") # save plot as png-file 
plt.show()

#ik heb dus 4464 metingen van de relative humidity en 4442 metingen van de cloud cover.
#nu dus de gemiddelde cloudcover berekenen bij een gegeven relative humidity interval van (2)%.
#hoe?