/plot_and_correct_humidity095b.py

'''
Created on Dec 2, 2016
script to plot and correct relative humidity data from 095b based on the plots from 
humidity_analysis.py. and time series plot from this script. also updates the netCDF file
now located on bunnyhill. plots up before and after histograms.

@author: pkormos
'''
 
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from pandas import Series, date_range
 
def rhstretch(data,old_low,new_low,old_high):
    slp = (1-new_low)/(old_high-old_low)
    itc = new_low - old_low * slp
    return(data*slp+itc)
 
stn0 = 'rc-095b'
wyoi = np.arange(1984, 2015)  # water years of interest
 
# bring in precipitation data
pnc = '/Volumes/megadozer/CZO/nc_working_files/rcew_ppt_pk.nc' # netcdf file name
pn = nc.Dataset(pnc, 'r')            # create precip ncdf file
pvar = pn.variables["precipitation"]
p_st = pn.variables['stations']    # id for stations
tt_list = ["station_%s" % i for i in range(np.shape(pvar)[1])]
p_list = [p_st.getncattr(tt_list[k]) for k in range(np.shape(pvar)[1])] # sorted rh station list
ppt = Series(np.ma.filled(pvar[:,p_list.index(stn0)],np.nan).ravel(),index=date_range('1962-1-1', '2014-12-31 23:00', freq='H')) # make pandas dataset
pn.close()
 
# # bring in relative humidity data
rnc = '/Volumes/megadozer/CZO/nc_working_files/rcew_rh_pk.nc' # netcdf file name
rn = nc.Dataset(rnc, 'r')            # create precip ncdf file
rvar = rn.variables["relative_humidity"]
rh_st = rn.variables['stations']    # id for stations
tt_list = ["station_%s" % i for i in range(np.shape(rvar)[1])]
r_list = [rh_st.getncattr(tt_list[k]) for k in range(np.shape(rvar)[1])] # sorted rh station list
rh = Series(np.ma.filled(rvar[:,r_list.index(stn0)],np.nan).ravel(),index=date_range('18-Jun-1981 11:00:00', '01-Oct-2014', freq='H')) # make pandas dataset
rn.close()

# # plot up data to check it
# fig0 = plt.figure(num=10, figsize=(17,8.5))
# ax3 = fig0.add_subplot(1,1,1)
# rh.plot()
# ax3.set_xlim('1999/10/1','2014/10/1')
# ax3.hlines(y=1,xmin='1983/10/1',xmax='2014/10/1')
# ax3.hlines(y=.1,xmin='1983/10/1',xmax='2014/10/1')

# plot original and corrected data
rh = rh.clip_upper(1)
t1 = '6/1/2001 0:00'
t2 = '9/30/2014 23:00'
fig2 = plt.figure(num=2, figsize=(14,8.5), dpi=100, facecolor='w')
ax12 = fig2.add_axes([.05, .05, .9, .9]) 
ppt[t1:t2].plot(style= 'g-', ax=ax12)
ax2 = ax12.twinx()
rh[t1:t2].plot(style = 'b-', ax=ax2) # plot time series before modification
l0 = ax2.hlines(y=1,xmin=t1,xmax = t2,color='k')
l1 = ax2.hlines(y=.1,xmin=t1,xmax = t2,color='k')

## correct first section of data
st1 = '04/01/2003 1:00'
nd1 = '10/01/2004 0:00'
strt  = 0.98  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='c-',ax=ax2)

## correct 2nd section of data
st1 = '10/01/2004 0:00'
nd1 = '10/01/2005 0:00'
strt  = 0.95  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='k-',ax=ax2)

## correct 3rd section of data
st1 = '10/01/2005 1:00'
nd1 = '10/01/2006 0:00'
strt  = 0.97  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='m-',ax=ax2)

## correct 4th section of data
st1 = '10/01/2006 1:00'
nd1 = '10/01/2007 0:00'
strt  = 0.93  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='c-',ax=ax2)

## correct 5th section of data
st1 = '10/01/2007 1:00'
nd1 = '10/01/2008 0:00'
strt  = 0.96  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='k-',ax=ax2)

## correct 6th section of data
st1 = '10/01/2011 1:00'
nd1 = '10/01/2012 0:00'
strt  = 0.98  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='m-',ax=ax2)

## correct 7th section of data
st1 = '10/01/2012 1:00'
nd1 = '10/01/2014 0:00'
strt  = 0.99  # what the upper rh does go up to (increase this.)
strb  = 0.1  # what the lower rh should drop to 
strbb = 0.1  # what the lower rh does drop to
l1 = ax2.hlines(y=strt, xmin=st1, xmax=nd1,color='r', linewidth=2) # plot line showing area to be changed
l2 = ax2.vlines(ymin=0,ymax=1,x=st1,color='r', linewidth=2) # plot line showing area to be changed
l3 = ax2.vlines(ymin=0,ymax=1,x=nd1,color='r',linewidth=2)
rh[st1:nd1] = rhstretch(rh[st1:nd1],strbb,strb,strt)
rh[st1:nd1] = rh[st1:nd1].clip_upper(1)
l4 = rh[st1:nd1].plot(style='c-',ax=ax2)

# plot histograms for 144 rh during precip events
year_st = 2003
fig1 = plt.figure(num=1, figsize=(8.5,14), dpi=100, facecolor='w')
for y in np.arange(year_st,year_st+12):    
    rh0  = rh['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y] # get rh time slice for stn0
    ppt0 = ppt['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y] # get rh time slice for stn0]
    pind = ppt0>0      # index to times with precip
    rind = np.isfinite(rh0)
    ind = pind&rind
    b_edg = np.linspace(0, 1, num=51)       # bin edges
    width = 0.02
    wts = ppt0[ind].dropna()
    whst = np.histogram(rh0[ind].dropna(),bins = b_edg,weights=wts)  # get weighted histogram bar heights
    ax1 = fig1.add_subplot(6,2,y-year_st+1)
    rects1 = ax1.bar(whst[1][0:np.size(whst[0])]+.01,whst[0], width, color='green') 
    plt.axis([0,1.01,0,plt.ylim()[1]])
    plt.title('WY%s R.H.'%(y),y=.83, x=0.4)
    plt.grid(True)
    plt.show()

## update the netcdf file with the corrected data
rn = nc.Dataset(rnc, 'r+')            # create precip ncdf file
rvar = rn.variables["relative_humidity"]
rvar[:,r_list.index(stn0)] = np.array(rh)
 
# a = rvar[:,r_list.index(stn0)]
# plt.plot(a)
rn.close()