/plot_and_correct_humidity076.py
Python | 239 lines | 91 code | 30 blank | 118 comment | 0 complexity | 40f864a0b58c63f968ead714ab3db2af MD5 | raw file
- '''
- Created on Oct 31, 2016
- script to plot and correct relative humidity data based on the plots from
- humidity_analysis.py.
- @author: pkormos
- '''
- import netCDF4 as nc
- import numpy as np
- import matplotlib.pyplot as plt
- from pandas import DataFrame, 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)
- wyoi = np.arange(1984, 2015) # water years of interest
- # bring in precipitation data
- p_list = np.array(['rc-047', 'rc-055', 'rc-059', 'rc-074', 'rc-075', 'rc-076', 'rc-078', 'rc-083', # manually bring in column headers
- 'rc-088', 'rc-095b', 'rc-097', 'rc-106', 'rc-114', 'rc-119', 'rc-128', 'rc.fl-057',
- 'rc.lsc-127', 'rc.mc-043041', 'rc.mc-053', 'rc.mc-054', 'rc.mc-061', 'rc.mc-072',
- 'rc.ng-098c', 'rc.ng-108', 'rc.sc-012', 'rc.sc-015', 'rc.sc-023', 'rc.sc-024',
- 'rc.sc-031','rc.sc-045', 'rc.sc.mp-033', 'rc.sum-049', 'rc.tg-116c','rc.tg-126',
- 'rc.tg-145', 'rc.tg-147', 'rc.tg-155', 'rc.tg-156', 'rc.tg-165', 'rc.tg-167',
- 'rc.tg-174_ppt', 'rc.tg.dc-144', 'rc.tg.dc-154', 'rc.tg.dc-163', 'rc.tg.dc-163',
- 'rc.tg.dc.jd-124', 'rc.tg.dc.jd-124b', 'rc.tg.dc.jd-125', 'rc.tg.rme-166b',
- 'rc.tg.rme-176', 'rc.tg.rme-rmsp', 'rc.tg.rmw-166x94','rc.usc-138031',
- 'rc.usc-138044', 'rc.usc-d03','rc.usc-j10', 'rc.usc-l21'])
- pnc = '/Volumes/bunnyhill/HUMIDITY_WORK_2016/rcew_ppt_pk.nc' # netcdf file name
- pn = nc.Dataset(pnc, 'r') # create precip ncdf file
- pvar = pn.variables["precipitation"]
- ppt = DataFrame(pvar[:], index=date_range('1962-1-1', '2014-12-31 23:00', freq='H'), columns = p_list) # make pandas dataset
- pn.close()
- # bring in relative humidity data
- r_list = np.array(['rc-076','rc-095b','rc-128','rc.lsc-127','rc.sc-012','rc.sc-031',
- 'rc.tg-145','rc.tg-167','rc.tg-167b','rc.tg-174','rc.tg.dc-144',
- 'rc.tg.dc-163',
- 'rc.tg.dc.jd-124',
- 'rc.tg.dc.jd-124b',
- 'rc.tg.dc.jd-125',
- 'rc.tg.rme-166b',
- 'rc.tg.rme-176',
- 'rc.tg.rme-rmsp',
- 'rc.usc-d03',
- 'rc.usc-j10',
- 'rc.usc-l21',
- 'rc-076_tf',
- 'rc-128_tf',
- 'rc.lsc-127_tf',
- 'rc.tg.rme-176_tf'])
- rnc = '/Volumes/bunnyhill/HUMIDITY_WORK_2016/rcew_rh_pk.nc' # netcdf file name
- rn = nc.Dataset(rnc, 'r') # create precip ncdf file
- rvar = rn.variables["relative_humidity"]
- rh = DataFrame(rvar[:], index=date_range('18-Jun-1981 11:00:00', '01-Oct-2014', freq='H'), columns = r_list) # make pandas dataset
- rn.close()
- # code to correct 076 from 9/22/1985 17:00 to 10/26/1995 0:00
- stn0 = 'rc-076'
- # # plot histograms for rh during precip events
- # fig1 = plt.figure(num=1, figsize=(8.5,14), dpi=100, facecolor='w')
- # for y in np.nditer(wyoi[0:14]):
- # rh0 = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
- # ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
- # pind = ppt0>0 # index to times with precip
- # b_edg = np.linspace(0, 1, num=51) # bin edges
- # width = 0.02
- # wts = ppt0[pind].dropna()
- # whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts) # get weighted histogram bar heights
- # ax1 = fig1.add_subplot(7,2,y-1983)
- # 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()
- # fig1.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_before.tif', dpi=100)
- # plot original and corrected data
- stretcher1=0.84
- stretcher2=0.88
- # fig2 = plt.figure(num=2, figsize=(14,8.5), dpi=100, facecolor='w')
- # ax2 = fig2.add_axes([.05, .05, .9, .9])
- # rh.ix['10/1/1982 0:00':'9/30/1999 23:00',[stn0]].plot(style = 'b-', ax=ax2) # plot time series before modification
- # ax2.hlines(y=stretcher1,xmin='9/22/1985 17:00', xmax='6/1/1992 12:00',color='r', linewidth=2) # plot line showing area to be changed
- # ax2.hlines(y=stretcher2,xmin='6/1/1992 13:00', xmax='10/26/1995 0:00',color='r', linewidth=2) # plot line showing area to be changed
- # xmin, xmax = ax2.get_xlim()
- # ax2.hlines(y=.13, xmin=xmin, xmax=xmax-1, color='r', linewidth=2)
- # apply stretch correction
- rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]] = rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]]/stretcher1
- rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]] = rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]].clip_upper(1)
- rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]] = rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]]/stretcher2
- rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]] = rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]].clip_upper(1)
- # rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]].plot(style = 'g-', ax=ax2) # plot time series after mod. green
- # rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]].plot(style = 'k-', ax=ax2) # plot time series after mod. green
- # # this method made the minimums too high compared to wy 1984-5 and wy96-99
- # rh.ix['9/22/1985 17:00':'9/30/1995 23:00',[stn0]] = rh.ix['9/22/1985 17:00':'9/30/1995 23:00',[stn0]]+0.15
- # rh.ix['9/22/1985 17:00':'9/30/1995 23:00',[stn0]] = rh.ix['9/22/1985 17:00':'9/30/1995 23:00',[stn0]].clip_upper(1)
- # ax12 = ax2.twinx()
- # ppt.ix['10/1/1982 0:00':'9/30/1999 23:00',[stn0]].plot(style= 'g-', ax=ax12)
- # fig2.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_076_before_after.tif', dpi=100)
- # ## replot histograms after correction
- # fig3 = plt.figure(num=3, figsize=(8.5,14), dpi=100, facecolor='w')
- # for y in np.nditer(wyoi[0:14]):
- # rh0 = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
- # ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
- # pind = ppt0>0 # index to times with precip
- # b_edg = np.linspace(0, 1, num=51) # bin edges
- # width = 0.02
- # wts = ppt0[pind].dropna()
- # whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts) # get weighted histogram bar heights
- # ax3 = fig3.add_subplot(7,2,y-1983)
- # rects1 = ax3.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)
- # fig3.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_after.tif', dpi=100)
- # # plot histograms for rh 2000 -> 2014 during precip events
- # fig4 = plt.figure(num=4, figsize=(8.5,14), dpi=100, facecolor='w')
- # for y in np.nditer(wyoi[17:31]):
- # rh0 = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
- # ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
- # pind = ppt0>0 # index to times with precip
- # b_edg = np.linspace(0, 1, num=51) # bin edges
- # width = 0.02
- # wts = ppt0[pind].dropna()
- # whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts) # get weighted histogram bar heights
- # ax1 = fig4.add_subplot(7,2,y-2000)
- # 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()
- # fig4.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_before_2001_2014.tif', dpi=100)
- ################################################################################################################
- ## plot original and corrected data 2003-2014
- ################################################################################################################
- fig5 = plt.figure(num=5, figsize=(14,8.5), dpi=100, facecolor='w')
- ax15 = fig5.add_axes([.05, .05, .9, .9])
- ppt.ix['10/1/1984 0:00':'9/30/2014 23:00',[stn0]].plot(style= 'g-', ax=ax15)
- ax5 = ax15.twinx()
- rh.ix['10/1/1984 0:00':'9/30/2014 23:00',[stn0]].plot(style = 'b-', ax=ax5) # plot time series before modification
- plt.legend(loc='center', bbox_to_anchor=(1, 0.5))
- ## correct first section of data 05-14
- stt='2005-10-01 01:00'
- ndt='2007-10-01 00:00'
- strt = 0.95 # what the upper rh does go up to
- strb = 0.00 # what the lower rh should drop to
- strbb = 0.00 # what the lower rh does drop to
- # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
- # correct first section of data and plot it 05-14
- rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
- rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
- # l2 = rh.ix[stt:ndt,[stn0]].plot(style='c-',ax=ax5)
- ## correct 2nd section of data 05-14
- stt='2007-10-01 01:00'
- ndt='2009-10-01 00:00'
- strt = 0.97 # what the upper rh does go up to
- strb = 0.00 # what the lower rh should drop to
- strbb = 0.00 # what the lower rh does drop to
- # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
- # correct 2nd section of data and plot it 05-14
- rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
- rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
- # l2 = rh.ix[stt:ndt,[stn0]].plot(style='y-',ax=ax5)
- ## correct 3rd section of data 05-14
- stt='2011-10-01 01:00'
- ndt='2012-10-01 00:00'
- strt = 0.90 # what the upper rh does go up to
- strb = 0.00 # what the lower rh should drop to
- strbb = 0.00 # what the lower rh does drop to
- # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
- # correct 3rd section of data and plot it 05-14
- rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
- rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
- # l2 = rh.ix[stt:ndt,[stn0]].plot(style='m-',ax=ax5)
- ## correct 4th section of data 05-14
- stt='2012-10-01 01:00'
- ndt='2014-10-01 00:00'
- strt = 0.94 # what the upper rh does go up to
- strb = 0.00 # what the lower rh should drop to
- strbb = 0.00 # what the lower rh does drop to
- # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
- # correct 4th section of data and plot it 05-14
- rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
- rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
- # l2 = rh.ix[stt:ndt,[stn0]].plot(style='m-',ax=ax5)
- #
- # fig5.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_076_before_after2015.tif', dpi=100)
- # fig6 = plt.figure(num=6, figsize=(8.5,14), dpi=100, facecolor='w')
- # for y in np.nditer(wyoi[17:31]):
- # rh0 = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
- # ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
- # pind = ppt0>0 # index to times with precip
- # b_edg = np.linspace(0, 1, num=51) # bin edges
- # width = 0.02
- # wts = ppt0[pind].dropna()
- # whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts) # get weighted histogram bar heights
- # ax1 = fig6.add_subplot(7,2,y-2000)
- # 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()
- # fig6.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_after_2001_2014.tif', dpi=100)
- rh[stn0] = rh[stn0].clip_upper(1) # clip all of the values greater than 1
- # rh[stn0].plot()
- ## 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==stn0] = np.array(rh[stn0])
- # a = rvar[:,r_list.index(stn0)]
- # plt.plot(a)
- rn.close()