<jupyter_start><jupyter_text>NREL WIND Toolkit - NSRDB Comparison This notebook compares the National Renewable Energy Laboratory (NREL) Wind Integration National Dataset (WIND) Toolkit and National Solar Radiation Database (NSRDB) data. The data is provided from Amazon Web Services using the HDF Group's Highly Scalable Data Service (HSDS).For this to work you must first install h5pyd:```pip install --user h5pyd```Next you'll need to configure HSDS:```hsconfigure```and enter at the prompt:```hs_endpoint = https://developer.nrel.gov/api/hsdshs_username = Nonehs_password = <PASSWORD>hs_api_key = <KEY>```*The example API key here is for demonstation and is rate-limited per IP. To get your own API key, visit https://developer.nrel.gov/signup/*You can also add the above contents to a configuration file at ~/.hscfg<jupyter_code>from collections import OrderedDict
import h5pyd
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import wilcoxon, ttest_rel
from scipy.spatial import cKDTree
from pyproj import Proj

#custom plotting properties
sns.set_style("white")
sns.set_style("ticks")
meanprops = dict(marker='o', markeredgecolor='black', markerfacecolor="None",
                 markersize=5)


def WTK_idx(wtk, lat_lon):
    """
    Function to find the nearest x/y WTK indices for a given lat/lon using Proj4 projection library
    
    Parameters
    ----------
    wtk : 'h5pyd.File'
        h5pyd File instance for the WTK
    lat_lon : tuple | list
        (lat, lon) coordinates of interest
        
    Results
    -------
    ij : 'tuple'
        x/y coordinate in the database of the closest pixel to coordinate of interest
    """
    dset_coords = wtk['coordinates']
    projstring = """+proj=lcc +lat_1=30 +lat_2=60 
                    +lat_0=38.47240422490422 +lon_0=-96.0 
                    +x_0=0 +y_0=0 +ellps=sphere 
                    +units=m +no_defs """
    projectLcc = Proj(projstring)
    origin_ll = reversed(dset_coords[0][0])  # Grab origin directly from database
    origin = projectLcc(*origin_ll)
    
    lon_lat = reversed(lat_lon)
    coords = projectLcc(*lon_lat)
    delta = np.subtract(coords, origin)
    ij = [int(round(x / 2000)) for x in delta]
    return tuple(reversed(ij))


def NSRDB_idx(nsrdb, lat_lon):
    """
    Function to find the NSRDB site index for a given lat/lon using a KDTree
    
    Parameters
    ----------
    nsrdb : 'h5pyd.File'
        h5pyd File instance for the NSRDB
    lat_lon : tuple | list
        (lat, lon) coordinates of interest
        
    Results
    -------
    ij : 'tuple'
        x/y coordinate in the database of the closest pixel to coordinate of interest
    """
    dset_coords = nsrdb['coordinates'][...]
    tree = cKDTree(dset_coords)
    dist, pos = tree.query(np.array(lat_lon))
    return pos<jupyter_output><empty_output><jupyter_text>Compare NSRDB to WIND Toolkit<jupyter_code># Specify a latitude and longitude
lat_lon = (32.785450, -96.800734) # Dallas, TX<jupyter_output><empty_output><jupyter_text>NSRDB Site<jupyter_code># Find nearest NSRDB site to desired lat_lon
with h5pyd.File('/nrel/nsrdb/nsrdb_2017.h5', 'r') as nsrdb:
    idx = NSRDB_idx(nsrdb, lat_lon)
    site_meta = pd.DataFrame(nsrdb['meta'][[idx]])
    nsrdb_lat_lon = nsrdb['coordinates'][idx]
    print('NSRDB nearest site coordinates = ', nsrdb_lat_lon)
    display(site_meta)
# Extract hourly irradiance, temperature, wind speed, and solar zenith angle from 2007-2013 NSRDB files
nsrdb_data = []
for year in range(2007, 2014):
    file = '/nrel/nsrdb/v3/nsrdb_{}.h5'.format(year)
    df = pd.DataFrame()
    with h5pyd.File(file, 'r') as nsrdb:
        for dset in ['ghi', 'dni', 'dhi', 'wind_speed', 'air_temperature', 'solar_zenith_angle']:
            ds = nsrdb[dset]
            df[dset] = ds[::2, idx] / ds.attrs['psm_scale_factor']
        
        df.index = pd.to_datetime(nsrdb['time_index'][::2].astype(str))
        df.index.name = 'datetime'
    
    nsrdb_data.append(df)
    
nsrdb_data = pd.concat(nsrdb_data)
nsrdb_data.head()<jupyter_output><empty_output><jupyter_text>WIND Toolkit Site<jupyter_code># Find nearest WTK site to desired lat_lon
with h5pyd.File("/nrel/wtk-us.h5", 'r') as wtk:
    wtk_ij = WTK_idx(wtk, lat_lon)
    wtk_lat_lon = wtk['coordinates'][wtk_ij[0], wtk_ij[1]]
    print('WTK nearest site coordinates = ', wtk_lat_lon)
# Extract hourly irradiance, temperature, and wind speed from WTK
with h5pyd.File("/nrel/wtk-us.h5", 'r') as wtk:
    dt = pd.to_datetime(wtk['datetime'][:].astype(str))
    wtk_data = OrderedDict()
    for var in ['DIF', 'DNI', 'GHI', 'temperature_2m', 'windspeed_10m']:
        wtk_data[var] = wtk[var][:, wtk_ij[0], wtk_ij[1]]

    wtk_data = pd.DataFrame(wtk_data, index=dt)
    wtk_data.head()<jupyter_output><empty_output><jupyter_text>Comparison Combine NSRDB and WTK datasets to allow for comparison<jupyter_code>nsrdb_df = nsrdb_data.copy()
nsrdb_df = nsrdb_df.rename(columns={'dhi': 'DHI', 'dni': 'DNI', 'ghi': 'GHI',
                                        'air_temperature': 'Temperature', 'wind_speed': 'Wind Speed'})
nsrdb_df['dataset'] = 'NSRDB'

wtk_df = wtk_data.copy()
wtk_df = wtk_df.rename(columns={'DIF': 'DHI', 'temperature_2m': 'Temperature', 'windspeed_10m': 'Wind Speed'})
wtk_df['dataset'] = 'WTK'
wtk_df['Temperature'] = wtk_df['Temperature'] - 273.15 # Temperature from C to K

comp_ds = pd.concat([nsrdb_df, wtk_df])
# Extract year, month and hour
comp_ds['year'] = comp_ds.index.year
comp_ds['month'] = comp_ds.index.month
# Shift hour to local time
comp_ds['hour'] = (comp_ds.index + pd.to_timedelta('{:}h'.format(site_meta['timezone'][0]))).hour

comp_ds.head()<jupyter_output>/Users/mrossol/miniconda3/envs/hsds/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  # This is added back by InteractiveShellApp.init_path()<jupyter_text>[Irradiance](Irradiance-Variables) | [Meteorological](Meterological-Variables)  Irradiance VariablesFor irradiance variables non-daylight hours need to be removed<jupyter_code># Daylight hours have a zenith angle <= 89 degrees
daylight = nsrdb_df['solar_zenith_angle'] <= 89
daylight = nsrdb_df.loc[daylight].index
cols = ['DNI', 'GHI', 'DHI', 'dataset', 'year', 'month', 'hour']
solar_comp = comp_ds.loc[daylight, cols].copy()
solar_comp.head()<jupyter_output><empty_output><jupyter_text>[GHI](Global-Horizontal-Irradiance) | [DHI](Diffuse-Horizontal-Irradiance) | [DNI](Direct-Normal-Irradiance)  Global Horizontal Irradiance<jupyter_code>variable = 'GHI'
label = 'Global Horizontal Irradiance (W m-2)'

nsrdb = solar_comp.loc[solar_comp['dataset'] == 'NSRDB', variable]
wtk = solar_comp.loc[solar_comp['dataset'] == 'WTK', variable]
fig, ax = plt.subplots(figsize=((6,5)), dpi=100)
sns.distplot(nsrdb, bins=100, ax=ax, kde=False, norm_hist=True)
sns.distplot(wtk, bins=100, ax=ax, kde=False, norm_hist=True)
sns.despine(offset=10, trim=False)
plt.ylabel('Probability')
plt.xlabel(label)
plt.legend(['NSRDB', 'WTK'])
plt.show()

p_val = wilcoxon(nsrdb, wtk)[1]
print('p-value = {:.2f}'.format(p_val))

fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Annual Statistics', fontsize=16)
sns.boxplot(x='year', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Year', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('year')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Monthly Statistics', fontsize=16)
sns.boxplot(x='month', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Month', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('month')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Diurnal Statistics', fontsize=16)
sns.boxplot(x='hour', y=variable, hue='dataset', data=solar_comp,
           showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Hour', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('hour')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))<jupyter_output><empty_output><jupyter_text>Take aways:- The WTK reports great GHI than the NSRDB.- The variability of GHI is also larger in the WTK than the NSRDB.- The diurnal pattern of GHI is consistent between the NSRDB and WTK.  Diffuse Horizontal Irradiance<jupyter_code>variable = 'DHI'
label = 'Diffuse Horizontal Irradiance (W m-2)'

nsrdb = solar_comp.loc[solar_comp['dataset'] == 'NSRDB', variable]
wtk = solar_comp.loc[solar_comp['dataset'] == 'WTK', variable]
fig, ax = plt.subplots(figsize=((6,5)), dpi=100)
sns.distplot(nsrdb, bins=100, ax=ax, kde=False, norm_hist=True)
sns.distplot(wtk, bins=100, ax=ax, kde=False, norm_hist=True)
sns.despine(offset=10, trim=False)
plt.ylabel('Probability')
plt.xlabel(label)
plt.legend(['NSRDB', 'WTK'])
plt.show()

p_val = wilcoxon(nsrdb, wtk)[1]
print('p-value = {:.2f}'.format(p_val))

fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Annual Statistics', fontsize=16)
sns.boxplot(x='year', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Year', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('year')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Monthly Statistics', fontsize=16)
sns.boxplot(x='month', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Month', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('month')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Diurnal Statistics', fontsize=16)
sns.boxplot(x='hour', y=variable, hue='dataset', data=solar_comp,
           showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Hour', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('hour')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))<jupyter_output><empty_output><jupyter_text>Take aways:- The NSRDB reports greather DHI than the WTK.- The variability of DHI is also larger in the NSRDB than the WTK.- The monthly pattern of DHI is not well correlated between the NSRDB and WTK.  Direct Normal Irradiance<jupyter_code>variable = 'DNI'
label = 'Direct Normal Irradiance (W m-2)'

nsrdb = solar_comp.loc[solar_comp['dataset'] == 'NSRDB', variable]
wtk = solar_comp.loc[solar_comp['dataset'] == 'WTK', variable]
fig, ax = plt.subplots(figsize=((6,5)), dpi=100)
sns.distplot(nsrdb, bins=100, ax=ax, kde=False, norm_hist=True)
sns.distplot(wtk, bins=100, ax=ax, kde=False, norm_hist=True)
sns.despine(offset=10, trim=False)
plt.ylabel('Probability')
plt.xlabel(label)
plt.legend(['NSRDB', 'WTK'])
plt.show()

p_val = wilcoxon(nsrdb, wtk)[1]
print('p-value = {:.2f}'.format(p_val))

fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Annual Statistics', fontsize=16)
sns.boxplot(x='year', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Year', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('year')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Monthly Statistics', fontsize=16)
sns.boxplot(x='month', y=variable, hue='dataset', data=solar_comp,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Month', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('month')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Diurnal Statistics', fontsize=16)
sns.boxplot(x='hour', y=variable, hue='dataset', data=solar_comp,
           showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Hour', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = solar_comp.groupby('hour')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))<jupyter_output><empty_output><jupyter_text>Take aways:- The average DNI is slightly higher in the WTK on an annual basis.- The median DNI is substantially higher in the WTK on an annual and monthly basis.- The variability of DNI is comparable between the NSRDB and WTK on an annuaul and monthly basis. - The montly pattern of DNI does not match up well between NSRDB and WTK.- Of significant note is the fact that DNI varies little in the WTK during the middle of the day,    when the variability is greatest in the NSRDB.  Conclusions for Irradiance VariablesThe discrepencies in irradiance values between NSRDB and WTK can be attributed to poor modeling of  clouds in the underlying WRF model used to derive the WTK.- WRF appears to underpredict the total cloud over.    - A known issue with WRF is the ability to represent thin clouds.- This is consistent with the higher GHI (and DNI), but lower DHI values- This also explains the lack of mid-day variability in DHI in the WTK data  Meterological Variables [Temperature](Surface-Temperature) | [Wind Speed](Surface-Wind-Speed)  Surface Temperature<jupyter_code>variable = 'Temperature'
label = 'Temperature (C)'

nsrdb = comp_ds.loc[comp_ds['dataset'] == 'NSRDB', variable]
wtk = comp_ds.loc[comp_ds['dataset'] == 'WTK', variable]
fig, ax = plt.subplots(figsize=((6,5)), dpi=100)
sns.distplot(nsrdb, bins=100, ax=ax, kde=False, norm_hist=True)
sns.distplot(wtk, bins=100, ax=ax, kde=False, norm_hist=True)
sns.despine(offset=10, trim=False)
plt.ylabel('Probability')
plt.xlabel(label)
plt.legend(['NSRDB', 'WTK'])
plt.show()

p_val = wilcoxon(nsrdb, wtk)[1]
print('p-value = {:.2f}'.format(p_val))

fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Annual Statistics', fontsize=16)
sns.boxplot(x='year', y=variable, hue='dataset', data=comp_ds,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Year', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('year')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Monthly Statistics', fontsize=16)
sns.boxplot(x='month', y=variable, hue='dataset', data=comp_ds,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Month', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('month')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Diurnal Statistics', fontsize=16)
sns.boxplot(x='hour', y=variable, hue='dataset', data=comp_ds,
           showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Hour', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('hour')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
with pd.option_context('display.max_rows', None, 'display.max_columns', 24):
    display(p_vals.round(2))<jupyter_output><empty_output><jupyter_text>Take aways:- The average Temperature is slightly higher in the WTK- The variability of Temperature is comparable between the NSRDB and the WTK,    though the distribution of tempertures in the WTK is slightly offest from the NSRDB.- The WTK suggests that it is slightly warmer in the mornings and evenings but shows the same    diurnal pattern as the NSRDB.  Surface Wind Speed<jupyter_code>variable = 'Wind Speed'
label = 'Wind Speed (m/s)'

nsrdb = comp_ds.loc[comp_ds['dataset'] == 'NSRDB', variable]
wtk = comp_ds.loc[comp_ds['dataset'] == 'WTK', variable]
fig, ax = plt.subplots(figsize=((6,5)), dpi=100)
sns.distplot(nsrdb, bins=100, ax=ax, kde=False, norm_hist=True)
sns.distplot(wtk, bins=100, ax=ax, kde=False, norm_hist=True)
sns.despine(offset=10, trim=False)
plt.ylabel('Probability')
plt.xlabel(label)
plt.legend(['NSRDB', 'WTK'])
plt.show()

p_val = wilcoxon(nsrdb, wtk)[1]
print('p-value = {:.2f}'.format(p_val))

fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Annual Statistics', fontsize=16)
sns.boxplot(x='year', y=variable, hue='dataset', data=comp_ds,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Year', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('year')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Monthly Statistics', fontsize=16)
sns.boxplot(x='month', y=variable, hue='dataset', data=comp_ds,
            showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Month', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('month')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
display(p_vals.round(2))


fig = plt.figure(figsize=(6,5), dpi=100)
axis = fig.add_subplot(111)
axis.set_title('Diurnal Statistics', fontsize=16)
sns.boxplot(x='hour', y=variable, hue='dataset', data=comp_ds,
           showfliers=False, ax=axis, showmeans=True, meanprops=meanprops)
sns.despine(offset=10, trim=False)
axis.set_xlabel('Hour', fontsize=14)
axis.set_ylabel(label, fontsize=14)
axis.tick_params(axis='x', labelsize=12, width=1, length=6)
axis.tick_params(axis='y', labelsize=12, width=1, length=6)
fig.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
plt.close()

p_vals = pd.DataFrame()
groups = comp_ds.groupby('hour')
for group, df in groups:
    nsrdb = df.loc[df['dataset'] == 'NSRDB', variable].values
    wtk = df.loc[df['dataset'] == 'WTK', variable].values
    p_vals.loc['p-value', group] = wilcoxon(nsrdb, wtk)[1]
    
with pd.option_context('display.max_rows', None, 'display.max_columns', 24):
    display(p_vals.round(2))<jupyter_output><empty_output><jupyter_text>Take aways:- The average Wind Speed is slightly higher in the WTK.- But, the WTK and NSRDB indicate the same variability in the wind speed.- The diurnal patters are poorly correlated between WTK and NSRDB.    The NSRDB indicated it is windiest during midday, while in the WTK    it is windiest in the evening.  Conclusions for Meteorological VariablesThe NSRDB's Temperature and Wind Speed is pulled from spatially coarse re-analysis data (MERRA).  The WTK is also derived from coarse re-analysis data (ERA-interim),  but is then resolved to a much high spatial resolution using the Weather Research and Forecasting Model (WRF).As temperature is not expected to vary significantly spatially it is expected that the NSRDB and WTK temperatures  should match up closely, which is indeed the case.  - The slight differences are likley due to differences in the underlying re-analysis data.Conversely wind speed does vary substantially over small spatial length scales and therefore it is unsurpising  that the NSRDB and WTK measurements for a given site differ.- Smoothing over the coarser grid would account for the lower average annual and monthly values in the NSRDB.- Differences in the dirunal pattern is likely due to spatial discrepencies in the source data,    i.e. the WTK pixel is not coincident with the MERRA pixel used by the NSRDB.  Statistical SignificanceTwo sided p-values<jupyter_code>annual_stats = pd.DataFrame()
groups = solar_comp.groupby('year')
for group, df in groups:
    for var in ['GHI', 'DHI', 'DNI']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        annual_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]

annual_means = solar_comp.groupby(['dataset', 'year']).mean()
for var in ['GHI', 'DHI', 'DNI']:
    nsrdb = annual_means.loc['NSRDB'][var].values
    wtk = annual_means.loc['WTK'][var].values
    annual_stats.loc['Annual Means', var] = ttest_rel(nsrdb, wtk)[1]


groups = comp_ds.groupby('year')
for group, df in groups:
    for var in ['Temperature', 'Wind Speed']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        annual_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]
        
annual_means = comp_ds.groupby(['dataset', 'year']).mean()
for var in ['Temperature', 'Wind Speed']:
    nsrdb = annual_means.loc['NSRDB'][var].values
    wtk = annual_means.loc['WTK'][var].values
    annual_stats.loc['Annual Means', var] = ttest_rel(nsrdb, wtk)[1]
    
annual_stats.index.name = 'year'
annual_stats.round(2)
monthly_stats = pd.DataFrame()
groups = solar_comp.groupby('month')
for group, df in groups:
    for var in ['GHI', 'DHI', 'DNI']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        monthly_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]

monthly_means = solar_comp.groupby(['dataset', 'month']).mean()
for var in ['GHI', 'DHI', 'DNI']:
    nsrdb = monthly_means.loc['NSRDB'][var].values
    wtk = monthly_means.loc['WTK'][var].values
    monthly_stats.loc['Monthly Means', var] = ttest_rel(nsrdb, wtk)[1]


groups = comp_ds.groupby('month')
for group, df in groups:
    for var in ['Temperature', 'Wind Speed']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        monthly_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]
        
monthly_means = comp_ds.groupby(['dataset', 'month']).mean()
for var in ['Temperature', 'Wind Speed']:
    nsrdb = monthly_means.loc['NSRDB'][var].values
    wtk = monthly_means.loc['WTK'][var].values
    monthly_stats.loc['Monthly Means', var] = ttest_rel(nsrdb, wtk)[1]
    
monthly_stats.index.name = 'month'
monthly_stats.round(2)
diurnal_stats = pd.DataFrame()
groups = solar_comp.groupby('hour')
for group, df in groups:
    for var in ['GHI', 'DHI', 'DNI']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        diurnal_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]

groups = comp_ds.groupby('hour')
for group, df in groups:
    for var in ['Temperature', 'Wind Speed']:
        nsrdb = df.loc[df['dataset'] == 'NSRDB', var].values
        wtk = df.loc[df['dataset'] == 'WTK', var].values
        diurnal_stats.loc[group, var] = wilcoxon(nsrdb, wtk)[1]
        
diurnal_stats = diurnal_stats.sort_index()
        
diurnal_means = solar_comp.groupby(['dataset', 'hour']).mean()
for var in ['GHI', 'DHI', 'DNI']:
    nsrdb = diurnal_means.loc['NSRDB'][var].values
    wtk = diurnal_means.loc['WTK'][var].values
    diurnal_stats.loc['Diurnal Means', var] = ttest_rel(nsrdb, wtk)[1]
        
diurnal_means = comp_ds.groupby(['dataset', 'hour']).mean()
for var in ['Temperature', 'Wind Speed']:
    nsrdb = diurnal_means.loc['NSRDB'][var].values
    wtk = diurnal_means.loc['WTK'][var].values
    diurnal_stats.loc['Diurnal Means', var] = ttest_rel(nsrdb, wtk)[1]
    
diurnal_stats.index.name = 'hour'
diurnal_stats.round(2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15747563191805716, 'fraction_numerical': 0.022963819593589955, 'mean_word_length': 5.162899613113419, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 50, 'https://': 2, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '321397', 'license': 'FSFAP', 'n_tokens_mistral': 11956, 'n_tokens_neox': 11250, 'n_words': 2488, 'path': 'notebooks/04_WTK-NSRDB_comparison.ipynb', 'repo_name': 'lkilcher/hsds-examples', 'size': '30294'}