<jupyter_start><jupyter_text>open ltm dictionary, plot and analyze variables from each of the three datasets<jupyter_code># data processing
import os
import pandas as pd, numpy as np, dask, json
import ogh
import geopandas as gpd
import ogh_xarray_landlab as oxl
import xarray as xr
import scipy.io

# data migration library
from utilities import hydroshare

# plotting and shape libraries
%matplotlib inline

import warnings
warnings.filterwarnings("ignore")
#hs=hydroshare.hydroshare()
#homedir = hs.getContentPath(os.environ["HS_RES_ID"])
homedir = '/home/jovyan/work/notebooks/GriddedMetDataDownload/projects/TemperatureSensorRegion'
os.chdir(homedir)
print('Data will be saved to:'+homedir)
# initialize ogh_meta
meta_file = dict(ogh.ogh_meta())
sorted(meta_file.keys())<jupyter_output><empty_output><jupyter_text>Load dictionary of data from all three datasets (livneh 2013, salathe 2014, PNNL 2018)<jupyter_code>import pickle 
infile = open('TempSensorArea_2001to2015.pkl','rb') #rb - read binary
ltm = pickle.load(infile)
infile.close()
#check variables saved in dictionary
ltm.keys()<jupyter_output><empty_output><jupyter_text>Load temperature observations from field<jupyter_code>filename = 'NCascades_AT.mat'

mat = scipy.io.loadmat(filename)

#range of dates to create new files for
start_date = '2007-09-01 00:30' 
end_date = '2008-09-30 23:00'

Hdates = [x.strftime('%Y-%m-%d %H:%M') for x in pd.date_range(start=start_date, end=end_date, freq='1.5H')]
Hddates = pd.to_datetime(Hdates)

start_dated = '2007-09-01' 
end_dated = '2008-09-30'
Ddates = [x.strftime('%Y-%m-%d') for x in pd.date_range(start=start_dated, end=end_dated, freq='D')]
Dddates = pd.to_datetime(Ddates)

#create hourly time series
Hdata = mat['AT']
dfh = pd.DataFrame(Hdata, columns = mat['sensor'].tolist(), index = Hddates )

#create daily time series Tmax
Ddata = mat['AT_maxT']
dfd_mx = pd.DataFrame(Ddata, columns = mat['sensor'].tolist(),index = Dddates )

#create daily time series Tmin
Ddata = mat['AT_minT']
dfd_min = pd.DataFrame(Ddata, columns = mat['sensor'].tolist(),index = Dddates )


#create daily time series Tmin
Ddata = mat['AT_meanT']
dfd_min = pd.DataFrame(Ddata, columns = mat['sensor'].tolist(),index = Dddates )

dfh.columns[12]<jupyter_output><empty_output><jupyter_text>Define wateryear summary functions (ogh.aggregate_space_time_sum and average currently use annual year)<jupyter_code>def WaterYear_Accum(df):
    
    """
    #df = S_PNNL_1981to1987['PREC_ACC_NC_hourlywrf_pnnl'] or any other dataframe in dictionary created by grid_clim_dict
    
    """
    #step 1, add water year (WY) column to data
    aa = df.index.to_series() #convert index to series so that lambda function can be applied 
    df['WY'] = aa.apply(lambda x: x.year if x.month < 10 else x.year+1)
    
    #step 2, remove datetime index to column, label column
    df_re = df.rename_axis('Date').reset_index()
    
    #step 3, move WY as index
    df_re.set_index('WY', inplace=True)
    
    #step 4, sum values in each column (grid point) by WY
    WYs_Sta = df_re.groupby(df_re.index).sum()
    
    #step 5, take mean of each grid cell as basin value
    WYs_basin = WYs_Sta.mean(axis=1)
    
    return(WYs_basin)


def WaterYear_Mean(df):
    
    """
    #df = S_PNNL_1981to1987['PREC_ACC_NC_hourlywrf_pnnl'] or any other dataframe in dictionary created by grid_clim_dict
    
    """
    #step 1, add water year (WY) column to data
    aa = df.index.to_series() #convert index to series so that lambda function can be applied 
    df['WY'] = aa.apply(lambda x: x.year if x.month < 10 else x.year+1)
    
    #step 2, remove datetime index to column, label column
    df_re = df.rename_axis('Date').reset_index()
    
    #step 3, move WY as index
    df_re.set_index('WY', inplace=True)
    
    #step 4, average values in each column (grid point) by WY
    WYm_Sta = df_re.groupby(df_re.index).mean()
    
    #step 5, take mean of each grid cell as basin value
    WYm_basin = WYm_Sta.mean(axis=1)
    
    return(WYm_basin)<jupyter_output><empty_output><jupyter_text>Analyze and plot variables<jupyter_code>#create precip, temperature and solar radiation variables
PNNL_hp = ltm['PREC_ACC_NC_hourlywrf_pnnl']
#water year accumulated precip
WYp = WaterYear_Accum(PNNL_hp)
#water year mean temp
PNNL_hT = ltm['T2_hourlywrf_pnnl']
WYt = WaterYear_Mean(PNNL_hT)
#water year mean solar radiation W/m2
PNNL_hRs = ltm['SWDOWN_hourlywrf_pnnl']
WYs = WaterYear_Mean(PNNL_hRs)<jupyter_output><empty_output><jupyter_text>Note: pnnl data is missing multi column index, add index back to pnnl variables<jupyter_code>#Tmx_p_d = ltm['T2_hourlywrf_pnnl'].resample('D').max()
ltm['T2_hourlywrf_pnnl'].head()
#PNNL_hp.head()
# check pnnl index
#reload column index
#RELOAD 1988 to 1994 because columns of orginal dataset also reset to index

infile = open('TempSensorArea_2001to2015_SepartePereods.pkl','rb') #rb - read binary
Ldict = pickle.load(infile)
infile.close()
os.chdir(homedir)

PNNL_COL = Ldict[0]['PREC_ACC_NC_hourlywrf_pnnl'].columns
#Ldict[0]['PREC_ACC_NC_hourlywrf_pnnl'].head(50).to_csv('PNNL_MultiIndex')
Ldict[1]['T2_hourlywrf_pnnl'].iloc[:,26].head(100)-273.15
# check liv and sal index
ltm['PRECIP_dailywrf_salathe2014'].head(50).to_csv('VIC_MultiIndex')
#Ldict[0]['PREC_ACC_NC_hourlywrf_pnnl'].head(50).to_csv('PNNL_MultiIndex')<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>## Daily Tmax

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_mx.iloc[:,0],'c.',linewidth=2.0,alpha=.2,label = 'b1_Tmax')
ax.plot(dfd_mx.iloc[:,6],'c.',linewidth=2.0,alpha=.5,label = 'b7_Tmax')
ax.plot(dfd_mx.iloc[:,8],'c.',linewidth=2.0,alpha=.8,label = 'b9_Tmax')
ax.plot(dfd_mx.iloc[:,10],'c.',linewidth=2.0,label = 'b11_Tmax')
ax.plot(ltm['TMAX_dailywrf_salathe2014'].iloc[:,12],'b.',linewidth=2.0,label = 'S_Tmax')
ax.plot(ltm['TMAX_dailymet_livneh2013'].iloc[:,12],'g.',linewidth=2.0,label = 'L_Tmax')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,26].resample('D').max()-273.15,'r.',linewidth=2.0,label = 'P_Tmax')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('tmax [degC]')
plt.title('Tmax, daily')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('LE_Tmax_d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
## 3 Day Tmax

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_mx.iloc[:,0].resample('3D').mean(),'c-',linewidth=2.0,alpha=.2,label = 'b1_Tmax')
ax.plot(dfd_mx.iloc[:,6].resample('3D').mean(),'c-',linewidth=2.0,alpha=.8,label = 'b7_Tmax')
ax.plot(dfd_mx.iloc[:,8].resample('3D').mean(),'c-',linewidth=2.0,alpha=.5,label = 'b9_Tmax')
ax.plot(dfd_mx.iloc[:,10].resample('3D').mean(),'c-',linewidth=2.0,label = 'b11_Tmax')
ax.plot(ltm['TMAX_dailywrf_salathe2014'].iloc[:,12].resample('3D').mean(),'b-',linewidth=2.0,label = 's_Tmax')
ax.plot(ltm['TMAX_dailymet_livneh2013'].iloc[:,12].resample('3D').mean(),'g-',linewidth=2.0,label = 'l_Tmax')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,26].resample('3D').max()-273.15,'r-',linewidth=2.0,label = 'P_Tmax')
ax.plot(Ldict[1]['T2_hourlywrf_pnnl'].iloc[:,26].resample('3D').max()-273.15,'y-',linewidth=2.0,label = 'P_Tmax')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('tmax [degC]')
plt.title('Tmax, 3 day average')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('LE_Tmax_3d.jpg', dpi = 300, bbox_inches='tight')

plt.show()

## Daily Tmin

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_min.iloc[:,0],'c.',linewidth=2.0,alpha=.2,label = 'b1_Tmin')
ax.plot(dfd_min.iloc[:,6],'c.',linewidth=2.0,alpha=.5,label = 'b7_Tmin')
ax.plot(dfd_min.iloc[:,8],'c.',linewidth=2.0,alpha=.8,label = 'b9_Tmin')
ax.plot(dfd_min.iloc[:,10],'c.',linewidth=2.0,label = 'b11_Tmin')
ax.plot(ltm['TMIN_dailywrf_salathe2014'].iloc[:,12],'b.',linewidth=2.0,label = 'S_Tmin')
ax.plot(ltm['TMIN_dailymet_livneh2013'].iloc[:,12],'g.',linewidth=2.0,label = 'L_Tmin')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,26].resample('D').min()-273.15,'r.',linewidth=2.0,label = 'P_Tmin')
#ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,26]-273.15,'r.',linewidth=2.0,label = 'P_Tmin')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmin [degC]')
plt.title('Tmin, daily')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('LE_Tmin_d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
## 3 Day Tmin

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_min.iloc[:,0].resample('3D').mean(),'c-',linewidth=2.0,alpha=.2,label = 'b1_Tmin')
ax.plot(dfd_min.iloc[:,6].resample('3D').mean(),'c-',linewidth=2.0,alpha=.8,label = 'b7_Tmin')
ax.plot(dfd_min.iloc[:,8].resample('3D').mean(),'c-',linewidth=2.0,alpha=.5,label = 'b9_Tmin')
ax.plot(dfd_min.iloc[:,10].resample('3D').mean(),'c-',linewidth=2.0,label = 'b11_Tmin')
ax.plot(ltm['TMIN_dailywrf_salathe2014'].iloc[:,12].resample('3D').mean(),'b-',linewidth=2.0,label = 's_Tmin')
ax.plot(ltm['TMIN_dailymet_livneh2013'].iloc[:,12].resample('3D').mean(),'g-',linewidth=2.0,label = 'l_Tmin')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,26].resample('3D').min()-273.15,'r-',linewidth=2.0,label = 'P_Tmin')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmin [degC]')
plt.title('Tmin, 3 day average')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('LE_Tmin_3d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
## High El, Daily Tmax

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_mx.iloc[:,12],'c.',linewidth=2.0,alpha=.8,label = 'b13_Tmax')
ax.plot(dfd_mx.iloc[:,14],'c.',linewidth=2.0,label = 'b17_Tmax')
ax.plot(ltm['TMAX_dailywrf_salathe2014'].iloc[:,29],'b.',linewidth=2.0,label = 'S_Tmax')
ax.plot(ltm['TMAX_dailymet_livneh2013'].iloc[:,29],'g.',linewidth=2.0,label = 'L_Tmax')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,5].resample('D').max()-273.15,'r.',linewidth=2.0,label = 'P_Tmax')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmax [degC]')
plt.title('Tmax, daily')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('UE_Tmax_d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
## 3 Day Tmax

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_mx.iloc[:,12].resample('3D').mean(),'c-',linewidth=2.0,alpha=.5,label = 'b13_Tmax')
ax.plot(dfd_mx.iloc[:,14].resample('3D').mean(),'c-',linewidth=2.0,label = 'b17_Tmax')
ax.plot(ltm['TMAX_dailywrf_salathe2014'].iloc[:,29].resample('3D').mean(),'b-',linewidth=2.0,label = 's_Tmax')
ax.plot(ltm['TMAX_dailymet_livneh2013'].iloc[:,29].resample('3D').mean(),'g-',linewidth=2.0,label = 'l_Tmax')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,5].resample('3D').max()-273.15,'r-',linewidth=2.0,label = 'P_Tmax')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmax [degC]')
plt.title('Tmax, 3 day average')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('UE_Tmax_3d.jpg', dpi = 300, bbox_inches='tight')

plt.show()

## Daily Tmin

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_min.iloc[:,12],'c.',linewidth=2.0,alpha=.8,label = 'b13_Tmin')
ax.plot(dfd_min.iloc[:,14],'c.',linewidth=2.0,label = 'b17_Tmin')
ax.plot(ltm['TMIN_dailywrf_salathe2014'].iloc[:,29],'b.',linewidth=2.0,label = 'S_Tmin')
ax.plot(ltm['TMIN_dailymet_livneh2013'].iloc[:,29],'g.',linewidth=2.0,label = 'L_Tmin')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,5].resample('D').min()-273.15,'r.',linewidth=2.0,label = 'P_Tmin')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmin [degC]')
plt.title('Tmin, daily')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('UE_Tmin_d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
## 3 Day Tmin

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(dfd_min.iloc[:,12].resample('3D').mean(),'c-',linewidth=2.0,alpha=.5,label = 'b13_Tmin')
ax.plot(dfd_min.iloc[:,14].resample('3D').mean(),'c-',linewidth=2.0,label = 'b17_Tmin')
ax.plot(ltm['TMIN_dailywrf_salathe2014'].iloc[:,29].resample('3D').mean(),'b-',linewidth=2.0,label = 's_Tmin')
ax.plot(ltm['TMIN_dailymet_livneh2013'].iloc[:,29].resample('3D').mean(),'g-',linewidth=2.0,label = 'l_Tmin')
ax.plot(ltm['T2_hourlywrf_pnnl'].iloc[:,5].resample('3D').min()-273.15,'r-',linewidth=2.0,label = 'P_Tmin')
# =============================================================================
# ax.plot(Tmx_S_d.index,Tmx_S_d,'g.',linewidth=2.0,label = 'Salathe')
# ax.plot(PNNL_hT_v_d_mx.index, PNNL_hT_v_d_mx,'r.',alpha = .2,linewidth=2.0,label = 'PNNL')
# 
# =============================================================================
plt.xlim(['2007-10-01','2008-09-30'])
# plt.ylim([0,1])
# plt.xticks(np.arange(0, 1.5, .25))
#plt.xlabel('date')
plt.ylabel('Tmin [degC]')
plt.title('Tmin, 3 day average')
#plt.grid()

legend = ax.legend(loc = 'lower right')
plt.savefig('UE_Tmin_3d.jpg', dpi = 300, bbox_inches='tight')

plt.show()
#TO DO
# ADD Elevaation of grid point, add more grid points
# Add elevation to tidbit
# Add more grid points, are all PNNL temps different from obs?
# Double check way pnnl netcdf fiels are appended. could it be that appended out of order?<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.267477584801555, 'fraction_numerical': 0.05097498275753966, 'mean_word_length': 6.776694295465626, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '468462', 'license': 'MIT', 'n_tokens_mistral': 7128, 'n_tokens_neox': 6612, 'n_words': 920, 'path': 'Ws_MetDn_Step3_Assess.ipynb', 'repo_name': 'keckje/Observatory', 'size': '15949'}