PageRenderTime 46ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 1ms

/plot_and_correct_humidity076.py

https://gitlab.com/pkormos/phase_python
Python | 239 lines | 91 code | 30 blank | 118 comment | 0 complexity | 40f864a0b58c63f968ead714ab3db2af MD5 | raw file
    

  1. '''
    2. 

  2. Created on Oct 31, 2016
    3. 

  3. 

  4. script to plot and correct relative humidity data based on the plots from 
    5. 

  5. humidity_analysis.py. 
    6. 

  6. 

  7. @author: pkormos
    8. 

  8. '''
    9. 

    10. 

  10. import netCDF4 as nc
    11. 

  11. import numpy as np
    12. 

  12. import matplotlib.pyplot as plt
    13. 

  13. from pandas import DataFrame, date_range
    14. 

  14. 

  15. def rhstretch(data,old_low,new_low,old_high):
    16. 

  16.     slp = (1-new_low)/(old_high-old_low)
    17. 

  17.     itc = new_low - old_low * slp
    18. 

  18.     return(data*slp+itc)
    19. 

  19. 

  20. wyoi = np.arange(1984, 2015)  # water years of interest
    21. 

  21. 

  22. # bring in precipitation data
    23. 

  23. 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
    24. 

  24.     'rc-088', 'rc-095b', 'rc-097', 'rc-106', 'rc-114', 'rc-119', 'rc-128', 'rc.fl-057',
    25. 

  25.     'rc.lsc-127', 'rc.mc-043041', 'rc.mc-053', 'rc.mc-054', 'rc.mc-061', 'rc.mc-072',
    26. 

  26.     'rc.ng-098c', 'rc.ng-108', 'rc.sc-012', 'rc.sc-015', 'rc.sc-023', 'rc.sc-024',
    27. 

  27.     'rc.sc-031','rc.sc-045', 'rc.sc.mp-033', 'rc.sum-049', 'rc.tg-116c','rc.tg-126',
    28. 

  28.     'rc.tg-145', 'rc.tg-147', 'rc.tg-155', 'rc.tg-156', 'rc.tg-165', 'rc.tg-167',
    29. 

  29.     'rc.tg-174_ppt', 'rc.tg.dc-144', 'rc.tg.dc-154', 'rc.tg.dc-163', 'rc.tg.dc-163',
    30. 

  30.     'rc.tg.dc.jd-124', 'rc.tg.dc.jd-124b', 'rc.tg.dc.jd-125', 'rc.tg.rme-166b',
    31. 

  31.     'rc.tg.rme-176', 'rc.tg.rme-rmsp', 'rc.tg.rmw-166x94','rc.usc-138031', 
    32. 

  32.     'rc.usc-138044', 'rc.usc-d03','rc.usc-j10', 'rc.usc-l21'])
    33. 

  33. pnc = '/Volumes/bunnyhill/HUMIDITY_WORK_2016/rcew_ppt_pk.nc' # netcdf file name
    34. 

  34. pn = nc.Dataset(pnc, 'r')            # create precip ncdf file
    35. 

  35. pvar = pn.variables["precipitation"]
    36. 

  36. ppt = DataFrame(pvar[:], index=date_range('1962-1-1', '2014-12-31 23:00', freq='H'), columns = p_list) # make pandas dataset
    37. 

  37. pn.close()
    38. 

  38. 

  39. # bring in relative humidity data
    40. 

  40. r_list = np.array(['rc-076','rc-095b','rc-128','rc.lsc-127','rc.sc-012','rc.sc-031',
    41. 

  41.     'rc.tg-145','rc.tg-167','rc.tg-167b','rc.tg-174','rc.tg.dc-144',
    42. 

  42.     'rc.tg.dc-163',
    43. 

  43.     'rc.tg.dc.jd-124',
    44. 

  44.     'rc.tg.dc.jd-124b',
    45. 

  45.     'rc.tg.dc.jd-125',
    46. 

  46.     'rc.tg.rme-166b',
    47. 

  47.     'rc.tg.rme-176',
    48. 

  48.     'rc.tg.rme-rmsp',
    49. 

  49.     'rc.usc-d03',
    50. 

  50.     'rc.usc-j10',
    51. 

  51.     'rc.usc-l21',
    52. 

  52.     'rc-076_tf',
    53. 

  53.     'rc-128_tf',
    54. 

  54.     'rc.lsc-127_tf',
    55. 

  55.     'rc.tg.rme-176_tf'])
    56. 

  56. rnc = '/Volumes/bunnyhill/HUMIDITY_WORK_2016/rcew_rh_pk.nc' # netcdf file name
    57. 

  57. rn = nc.Dataset(rnc, 'r')            # create precip ncdf file
    58. 

  58. rvar = rn.variables["relative_humidity"]
    59. 

  59. rh = DataFrame(rvar[:], index=date_range('18-Jun-1981 11:00:00', '01-Oct-2014', freq='H'), columns = r_list) # make pandas dataset
    60. 

  60. rn.close()
    61. 

  61. 

  62. # code to correct 076 from 9/22/1985 17:00 to 10/26/1995 0:00
    63. 

  63. stn0 = 'rc-076'
    64. 

  64. 

  65. # # plot histograms for rh during precip events
    66. 

  66. # fig1 = plt.figure(num=1, figsize=(8.5,14), dpi=100, facecolor='w')
    67. 

  67. # for y in np.nditer(wyoi[0:14]):    
    68. 

  68. #     rh0  = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
    69. 

  69. #     ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
    70. 

  70. #     pind = ppt0>0      # index to times with precip
    71. 

  71. #     b_edg = np.linspace(0, 1, num=51)       # bin edges
    72. 

  72. #     width = 0.02
    73. 

  73. #     wts = ppt0[pind].dropna()
    74. 

  74. #     whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts)  # get weighted histogram bar heights
    75. 

  75. #     ax1 = fig1.add_subplot(7,2,y-1983)
    76. 

  76. #     rects1 = ax1.bar(whst[1][0:np.size(whst[0])]+.01,whst[0], width, color='green') 
    77. 

  77. #     plt.axis([0,1.01,0,plt.ylim()[1]])
    78. 

  78. #     plt.title('WY%s R.H.'%(y),y=.83, x=0.4)
    79. 

  79. #     plt.grid(True)
    80. 

  80. #     plt.show()
    81. 

  81. # fig1.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_before.tif', dpi=100) 
    82. 

  82. 

  83. # plot original and corrected data
    84. 

  84. stretcher1=0.84
    85. 

  85. stretcher2=0.88
    86. 

  86. # fig2 = plt.figure(num=2, figsize=(14,8.5), dpi=100, facecolor='w')
    87. 

  87. # ax2 = fig2.add_axes([.05, .05, .9, .9]) 
    88. 

  88. # rh.ix['10/1/1982 0:00':'9/30/1999 23:00',[stn0]].plot(style = 'b-', ax=ax2) # plot time series before modification
    89. 

  89. # 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
    90. 

  90. # 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
    91. 

  91. # xmin, xmax = ax2.get_xlim()
    92. 

  92. # ax2.hlines(y=.13, xmin=xmin, xmax=xmax-1, color='r', linewidth=2)
    93. 

  93. 

  94. # apply stretch correction
    95. 

  95. 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
    96. 

  96. 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)
    97. 

  97. 

  98. 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
    99. 

  99. 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)
    100. 

  100. 

  101. # rh.ix['9/22/1985 17:00':'6/1/1992 12:00',[stn0]].plot(style = 'g-', ax=ax2) # plot time series after mod. green
    102. 

  102. # rh.ix['6/1/1992 13:00':'10/26/1995 0:00',[stn0]].plot(style = 'k-', ax=ax2) # plot time series after mod. green
    103. 

  103. 

  104. # # this method made the minimums too high compared to wy 1984-5 and wy96-99
    105. 

  105. # 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 
    106. 

  106. # 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)
    107. 

  107. 

  108. # ax12 = ax2.twinx()
    109. 

  109. # ppt.ix['10/1/1982 0:00':'9/30/1999 23:00',[stn0]].plot(style= 'g-', ax=ax12)
    110. 

  110. # fig2.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_076_before_after.tif', dpi=100) 
    111. 

  111. 

  112. # ## replot histograms after correction
    113. 

  113. # fig3 = plt.figure(num=3, figsize=(8.5,14), dpi=100, facecolor='w')
    114. 

  114. # for y in np.nditer(wyoi[0:14]):    
    115. 

  115. #     rh0  = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
    116. 

  116. #     ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
    117. 

  117. #     pind = ppt0>0      # index to times with precip
    118. 

  118. #     b_edg = np.linspace(0, 1, num=51)       # bin edges
    119. 

  119. #     width = 0.02
    120. 

  120. #     wts = ppt0[pind].dropna()
    121. 

  121. #     whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts)  # get weighted histogram bar heights
    122. 

  122. #     ax3 = fig3.add_subplot(7,2,y-1983)
    123. 

  123. #     rects1 = ax3.bar(whst[1][0:np.size(whst[0])]+.01,whst[0], width, color='green') 
    124. 

  124. #     plt.axis([0,1.01,0,plt.ylim()[1]])
    125. 

  125. #     plt.title('WY%s R.H.'%(y),y=.83, x=0.4)
    126. 

  126. # fig3.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_after.tif', dpi=100) 
    127. 

  127. 

  128. # # plot histograms for rh 2000 -> 2014 during precip events
    129. 

  129. # fig4 = plt.figure(num=4, figsize=(8.5,14), dpi=100, facecolor='w')
    130. 

  130. # for y in np.nditer(wyoi[17:31]):    
    131. 

  131. #     rh0  = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
    132. 

  132. #     ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
    133. 

  133. #     pind = ppt0>0      # index to times with precip
    134. 

  134. #     b_edg = np.linspace(0, 1, num=51)       # bin edges
    135. 

  135. #     width = 0.02
    136. 

  136. #     wts = ppt0[pind].dropna()
    137. 

  137. #     whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts)  # get weighted histogram bar heights
    138. 

  138. #     ax1 = fig4.add_subplot(7,2,y-2000)
    139. 

  139. #     rects1 = ax1.bar(whst[1][0:np.size(whst[0])]+.01,whst[0], width, color='green') 
    140. 

  140. #     plt.axis([0,1.01,0,plt.ylim()[1]])
    141. 

  141. #     plt.title('WY%s R.H.'%(y),y=.83, x=0.4)
    142. 

  142. #     plt.grid(True)
    143. 

  143. #     plt.show()
    144. 

  144. # fig4.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_before_2001_2014.tif', dpi=100) 
    145. 

  145. 

  146. 

  147. ################################################################################################################
    148. 

  148. ## plot original and corrected data 2003-2014
    149. 

  149. ################################################################################################################
    150. 

  150. fig5 = plt.figure(num=5, figsize=(14,8.5), dpi=100, facecolor='w')
    151. 

  151. ax15 = fig5.add_axes([.05, .05, .9, .9]) 
    152. 

  152. ppt.ix['10/1/1984 0:00':'9/30/2014 23:00',[stn0]].plot(style= 'g-', ax=ax15)
    153. 

  153. ax5 = ax15.twinx()
    154. 

  154. rh.ix['10/1/1984 0:00':'9/30/2014 23:00',[stn0]].plot(style = 'b-', ax=ax5) # plot time series before modification
    155. 

  155. plt.legend(loc='center', bbox_to_anchor=(1, 0.5))
    156. 

  156. 

  157. ## correct first section of data 05-14
    158. 

  158. stt='2005-10-01 01:00'
    159. 

  159. ndt='2007-10-01 00:00'
    160. 

  160. strt  = 0.95 # what the upper rh does go up to
    161. 

  161. strb  = 0.00 # what the lower rh should drop to 
    162. 

  162. strbb = 0.00 # what the lower rh does drop to
    163. 

  163. # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
    164. 

  164. 

  165. # correct first section of data and plot it 05-14
    166. 

  166. rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
    167. 

  167. rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
    168. 

  168. # l2 = rh.ix[stt:ndt,[stn0]].plot(style='c-',ax=ax5)
    169. 

  169. 

  170. ## correct 2nd section of data 05-14
    171. 

  171. stt='2007-10-01 01:00'
    172. 

  172. ndt='2009-10-01 00:00'
    173. 

  173. strt  = 0.97 # what the upper rh does go up to
    174. 

  174. strb  = 0.00 # what the lower rh should drop to 
    175. 

  175. strbb = 0.00 # what the lower rh does drop to
    176. 

  176. # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
    177. 

  177. 

  178. # correct 2nd section of data and plot it 05-14
    179. 

  179. rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
    180. 

  180. rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
    181. 

  181. # l2 = rh.ix[stt:ndt,[stn0]].plot(style='y-',ax=ax5)
    182. 

  182. 

  183. ## correct 3rd section of data 05-14
    184. 

  184. stt='2011-10-01 01:00'
    185. 

  185. ndt='2012-10-01 00:00'
    186. 

  186. strt  = 0.90 # what the upper rh does go up to
    187. 

  187. strb  = 0.00 # what the lower rh should drop to 
    188. 

  188. strbb = 0.00 # what the lower rh does drop to
    189. 

  189. # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
    190. 

  190. 

  191. # correct 3rd section of data and plot it 05-14
    192. 

  192. rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
    193. 

  193. rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
    194. 

  194. # l2 = rh.ix[stt:ndt,[stn0]].plot(style='m-',ax=ax5)
    195. 

  195. 

  196. ## correct 4th section of data 05-14
    197. 

  197. stt='2012-10-01 01:00'
    198. 

  198. ndt='2014-10-01 00:00'
    199. 

  199. strt  = 0.94 # what the upper rh does go up to
    200. 

  200. strb  = 0.00 # what the lower rh should drop to 
    201. 

  201. strbb = 0.00 # what the lower rh does drop to
    202. 

  202. # l1 = ax5.hlines(y=strt,xmin=stt, xmax=ndt,color='r', linewidth=2) # plot line showing area to be changed
    203. 

  203. 

  204. # correct 4th section of data and plot it 05-14
    205. 

  205. rh.ix[stt:ndt,[stn0]] = rhstretch(rh.ix[stt:ndt,[stn0]],strbb,strb,strt)
    206. 

  206. rh.ix[stt:ndt,[stn0]] = rh.ix[stt:ndt,[stn0]].clip_upper(1)
    207. 

  207. # l2 = rh.ix[stt:ndt,[stn0]].plot(style='m-',ax=ax5)
    208. 

  208. #
    209. 

  209. # fig5.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_076_before_after2015.tif', dpi=100) 
    210. 

  210. 

  211. # fig6 = plt.figure(num=6, figsize=(8.5,14), dpi=100, facecolor='w')
    212. 

  212. # for y in np.nditer(wyoi[17:31]):    
    213. 

  213. #     rh0  = rh.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0
    214. 

  214. #     ppt0 = ppt.ix['10/1/%s 0:00'%(y-1):'9/30/%s 23:00'%y,[stn0]] # get rh time slice for stn0]
    215. 

  215. #     pind = ppt0>0      # index to times with precip
    216. 

  216. #     b_edg = np.linspace(0, 1, num=51)       # bin edges
    217. 

  217. #     width = 0.02
    218. 

  218. #     wts = ppt0[pind].dropna()
    219. 

  219. #     whst = np.histogram(rh0[pind].dropna(),bins = b_edg,weights=wts)  # get weighted histogram bar heights
    220. 

  220. #     ax1 = fig6.add_subplot(7,2,y-2000)
    221. 

  221. #     rects1 = ax1.bar(whst[1][0:np.size(whst[0])]+.01,whst[0], width, color='green') 
    222. 

  222. #     plt.axis([0,1.01,0,plt.ylim()[1]])
    223. 

  223. #     plt.title('WY%s R.H.'%(y),y=.83, x=0.4)
    224. 

  224. #     plt.grid(True)
    225. 

  225. #     plt.show()
    226. 

  226. # fig6.savefig('/Users/pkormos/rcew_met_dist/snow_data/phase_analysis/phase_python/rh_hist_076_after_2001_2014.tif', dpi=100) 
    227. 

  227. 

  228. rh[stn0] = rh[stn0].clip_upper(1) # clip all of the values greater than 1
    229. 

  229. # rh[stn0].plot()
    230. 

  230. 

  231. 

  232. ## update the netcdf file with the corrected data
    233. 

  233. rn = nc.Dataset(rnc, 'r+')            # create precip ncdf file
    234. 

  234. rvar = rn.variables["relative_humidity"]
    235. 

  235. rvar[:,r_list==stn0] = np.array(rh[stn0])
    236. 

  236. 

  237. # a = rvar[:,r_list.index(stn0)]
    238. 

  238. # plt.plot(a)
    239. 

  239. rn.close()
    240. 

  240.