In the first notebook we calculated the clearsky irradiance for seven SURFRAD stations. Then we compared different models and atmospheric data sets at each station. In this notebook we will now combine all seven SURFRAD stations to compare models accross all regions. Finally we will calculate the grand total MBE and RMS for all SURFRAD stations for each model. The calculated clearsky data is shared in a Public OneDrive folder as hdf5, one file per site, each about 100MB. You can use ViTables or Java HDFView to visualize the data. This notebook assumes the files are all downloaded here, in this directory.
As described in the first notebook, this is a Jupyter notebook. To use it you may need to make sure that Python and the required Python packages are installed.
While process the SURFRAD station data, there was an issue with the Sioux Falls, SD data set from 2006 which contains some data from 2007 from October 8-21. The data appears to be similar to the data from 2007, so it was removed in favor of the the later.
In [54]:
# imports and settings
import os
import h5py
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import pvlib
import seaborn as sns
import statsmodels.api as sm
from pvsc44_clearsky_aod import ecmwf_macc_tools
%matplotlib inline
sns.set_context('notebook', rc={'figure.figsize': (8, 6)})
sns.set(font_scale=1.5)
In [55]:
# get the "metadata" that contains the station id codes for the SURFRAD data that was analyzed
METADATA = pd.read_csv('metadata.csv', index_col=0)
In [56]:
# load calculations for each station
atm_params_3min_clear = {}
for station_id in METADATA.index:
with h5py.File('%s_3min_clear_atm_params.h5' % station_id, 'r') as f:
np_atm_params_3min_clear = pd.DataFrame(np.array(f['data']))
np_atm_params_3min_clear['index'] = pd.DatetimeIndex(np_atm_params_3min_clear['index'])
np_atm_params_3min_clear.set_index('index', inplace=True)
np_atm_params_3min_clear.index.rename('timestamps', inplace=True)
atm_params_3min_clear[station_id] = np_atm_params_3min_clear
In [57]:
# filter out low light
# CONSTANTS
MODELS = {'solis': 'SOLIS', 'lt': 'Linke', 'macc': 'ECMWF-MACC', 'bird': 'Bird'}
CS = ['dni', 'dhi', 'ghi']
LOW_LIGHT = 200 # threshold for low light in W/m^2
is_bright = {}
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
is_bright[station_id] = station_atm_params_3min_clear['ghi'] > LOW_LIGHT
For a given site lets compare all models annual mean bias error (MBE) and root-mean-square error (RMSE). Then lets copmare all models MBE and RMSE for years by month to see if there are any seasonal biases. We can also just pivot the table by months to plot all years by month on the same plot. This should be similar to the combination of the MBE and RMSE on the same plot because the spread in montly MBE across years is representative of the RMSE for all years at that month. Finally we should make some boxplots to compmare models against each other to see if there is a statistical difference.
Absolute mean bias error is given as: $$MBE = \frac{\sum^{N}_{n=1} \left( predicted - measured \right)}{N}$$ Relative MBE is given as $$rMBE = N \frac{MBE}{\sum^{N}_{n=1} measured }$$
Absolute root mean square error is given as: $$RMSE = \sqrt{ \frac{\sum^{N}_{n=1} \left( predicted - measured \right)^2}{N} }$$ Relative RMSE is given as $$rRMSE = N \frac{RMSE}{\sum^{N}_{n=1} measured }$$
When measured and predicted numbers are very small, relative differences can be exagerated relative to a characteristic value. For example, a characteristic GHI might be 1000 W/m2. An error 100 W/m2 at say 800 W/m2 is only 12.5%, but at 200 W/m2 that would suddenly be 50%. There are at least two ways to avoid the curse of small numbers:
Use a characteristice value in the denominator. For example if you use the characteristic value of 1000 W/m2 for both errors you get 10% for 100 W/m2 at 800 W/m2 and at 200 W/m2, so the curse is avoided.
Rollup timeseries up to periods that include both small and big numbers. This doesn't minimizes the exageration of relative differences in small numbers because the demominator is dominated by larger numbers. For example in a typical day there might be 6 kWh/m2, an error of 0.6 kWh/m2 is only 10%, even if most of those errors occured during low light conditions.
Note that in the equation for rMBE, the number of samples N is in the numerator and the denominator so it cancels out allowing you to just sum up the timeseries over the desired period.
In [5]:
# plot annual MBE and RMSE, combine all models on one chart for each station and clear sky component
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
f, ax = plt.subplots(2, 3, figsize=(24, 12), sharex=False)
for n, cs in enumerate(CS):
annual_avg = station_atm_params_3min_clear[cs][is_bright[station_id]].resample('A').mean()
# MBE
for model in MODELS.iterkeys():
((
(station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]]
).resample('A').mean() / annual_avg).plot(ax=ax[0][n])
ax[0][n].legend(MODELS.values())
ax[0][n].set_title('Mean Bias Error (MBE) of %s at %s' % (cs.upper(), METADATA['station name'][station_id]))
ax[0][n].set_ylabel('average relative difference (arb. units)')
# RMSE
for model in MODELS.iterkeys():
(np.sqrt((
((station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]])**2
).resample('A').mean()) / annual_avg).plot(ax=ax[1][n])
ax[1][n].legend(MODELS.values())
ax[1][n].set_title('Root Mean Square Error (RMSE) of %s at %s' % (cs.upper(), METADATA['station name'][station_id]))
ax[1][n].set_ylabel('relative RMSE (arb. units)')
f.tight_layout()
f.savefig('%s_annual_mbe-rmse.png' % station_id)
In [6]:
# plot MBE and RMSE by month, combine all models on one chart for each station and clear sky component
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
f, ax = plt.subplots(2, 3, figsize=(24, 12), sharex=False)
for n, cs in enumerate(CS):
monthly_avg = station_atm_params_3min_clear[cs][is_bright[station_id]].groupby(lambda x: x.month).mean()
for model in MODELS.iterkeys():
((
(station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]]
).groupby(lambda x: x.month).mean() / monthly_avg).plot(ax=ax[0][n])
ax[0][n].legend(MODELS.values())
ax[0][n].set_title('rMBE of %s by month at %s' % (cs.upper(), METADATA['station name'][station_id]))
ax[0][n].set_ylabel('average relative difference (arb. units)')
for model in MODELS.iterkeys():
(np.sqrt((
((station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]])**2
).groupby(lambda x: x.month).mean()) / monthly_avg).plot(ax=ax[1][n])
ax[1][n].legend(MODELS.values())
ax[1][n].set_title('rRMSE of %s by month at %s' % (cs.upper(), METADATA['station name'][station_id]))
ax[1][n].set_ylabel('relative RMSE (arb. units)')
f.tight_layout()
f.savefig('%s_mbe-rmse_by_month.png' % station_id)
In [7]:
# combine MBE all years on one chart for each station, model and clear sky component
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
f, ax = plt.subplots(4, 3, figsize=(24, 24), sharex=False)
for n, cs in enumerate(CS):
for m, model in enumerate(MODELS.iteritems()):
x = station_atm_params_3min_clear.loc[is_bright[station_id]]
x.insert(1, 'mbe', (x['%s_%s' % (model[0], cs)] - x[cs]))
y = x.resample('M').mean()
y.insert(2, 'rmbe', y['mbe'] / y[cs])
y.pivot(index='month', columns='year', values='rmbe').plot(ax=ax[m][n])
ax[m][n].set_title('%s %s rMBE by month for %s'% (model[1], cs.upper(), METADATA['station name'][station_id]))
ax[m][n].set_ylabel('average relative difference (arb. units)')
f.tight_layout()
f.savefig('%s_%s_%s_mbe_by_month_all_years.png' % (station_id, model[0], cs))
In [8]:
# plot annual average error, combine all stations on one chart for each model and clear sky component
for model, name in MODELS.iteritems():
f, ax = plt.subplots(2, 3, figsize=(24, 12), sharex=False)
for n, cs in enumerate(CS):
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
# MBE
annual_avg = station_atm_params_3min_clear[cs][is_bright[station_id]].resample('A').mean()
((
(station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]]
).resample('A').mean() / annual_avg).plot(ax=ax[0][n])
# RMSE
(np.sqrt((
((station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]])**2
).resample('A').mean()) / annual_avg).plot(ax=ax[1][n])
ax[0][n].legend(METADATA.index.tolist())
ax[0][n].set_title('%s %s MBE' % (name, cs.upper()))
ax[0][n].set_ylabel('average relative difference (arb. units)')
ax[1][n].legend(METADATA.index.tolist())
ax[1][n].set_title('%s %s RMSE' % (name, cs.upper()))
ax[1][n].set_ylabel('relative RMSE (arb. units)')
f.tight_layout()
plt.savefig('%s_annual_mbe-rms.png' % model)
In [9]:
# plot monthly average error, combine all stations on one chart for each model and clear sky component
for model, name in MODELS.iteritems():
f, ax = plt.subplots(2, 3, figsize=(24, 12), sharex=False)
for n, cs in enumerate(CS):
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
monthly_avg = station_atm_params_3min_clear[cs][is_bright[station_id]].groupby(lambda x: x.month).mean()
((
(station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]]
).groupby(lambda x: x.month).mean() / monthly_avg).plot(ax=ax[0][n]).plot()
(np.sqrt((
((station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]])**2
).groupby(lambda x: x.month).mean()) / monthly_avg).plot(ax=ax[1][n])
ax[0][n].legend(METADATA.index.tolist())
ax[0][n].set_title('%s %s MBE by month' % (name, cs.upper()))
ax[0][n].set_ylabel('average relative difference (arb. units)')
ax[1][n].legend(METADATA.index.tolist())
ax[1][n].set_title('%s %s RMSE by month' % (name, cs.upper()))
ax[1][n].set_ylabel('relative RMSE (arb. units)')
f.tight_layout()
f.savefig('%s_%s_mbe-rmse_by_month.png' % (model, cs))
In [58]:
# CONSTANTS
MODELS = {0: 'solis', 1: 'lt', 2: 'macc', 3: 'bird'}
# an empty data frame for monthly atmospheric parameters
MONTHLY_ATM_PARAMS_3MIN_CLEAR = None
# loop over stations
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
# filter for bright days
data = station_atm_params_3min_clear[is_bright[station_id]]
# insert errors first
for n, model in MODELS.iteritems():
for m, cs in enumerate(CS):
data.insert(
43 + m + 3*n,
'%s-%s_err' % (model, cs.upper()),
data['%s_%s' % (model, cs)] - data[cs]
)
# average by month
avg_data = data.resample('M').mean()
# divide by montly average
for n, cs in enumerate(CS):
for m, model in MODELS.iteritems():
avg_data.insert(
54 + m + 4*n,
'%s-%s_rel' % (model, cs.upper()),
avg_data['%s-%s_err' % (model, cs.upper())] / avg_data[cs]
)
# insert station id
avg_data.insert(0, 'station', station_id)
# append together
if MONTHLY_ATM_PARAMS_3MIN_CLEAR is None:
MONTHLY_ATM_PARAMS_3MIN_CLEAR = pd.DataFrame(avg_data)
else:
MONTHLY_ATM_PARAMS_3MIN_CLEAR = MONTHLY_ATM_PARAMS_3MIN_CLEAR.append(avg_data)
MONTHLY_ATM_PARAMS_3MIN_CLEAR[['station', 'year', 'month', 'ghi', 'solis-GHI_err', 'solis-GHI_rel']]
Out[58]:
In [60]:
# plot
f, ax = plt.subplots(1, 2, figsize=(16, 6), sharex=False, sharey=False) # 24 x 6 for 3 subplots
for n, cs in enumerate(CS):
if cs == 'dhi': continue
if n == 2: n = 1
data = pd.melt(
MONTHLY_ATM_PARAMS_3MIN_CLEAR,
id_vars=['station', 'year', 'month'],
value_vars=['%s-%s_rel' % (model, cs.upper()) for model in MODELS.itervalues()],
var_name='model-cs',
value_name='error'
)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='red')
sns.boxplot(x='model-cs', y='error', hue='year', data=data, ax=ax[n],
whis=[5,95], showmeans=True, meanline=True, meanprops=meanlineprops)
ax[n].set_title('Interannual variability in %s by model' % cs.upper())
ax[n].set_ylabel('average monthly relative error')
ax[n].set_xlabel('model and clear sky component')
if n == 1:
ax[n].legend(bbox_to_anchor=(-0.1, -0.25), ncol=10, loc='lower center') # x = 0.5 for all 3 components
else:
legend = ax[n].legend()
legend.set_visible(False)
legend.remove()
f.tight_layout()
f.subplots_adjust(bottom=0.2)
f.savefig('y2y-boxplot_by_model-cs_CORRECTED.png')
In [61]:
# plot
f, ax = plt.subplots(1, 2, figsize=(16, 6), sharex=False, sharey=False) # 24 x 6 for 3 subplots
for n, cs in enumerate(CS):
if cs == 'dhi': continue
if n == 2: n = 1
data = pd.melt(
MONTHLY_ATM_PARAMS_3MIN_CLEAR,
id_vars=['station', 'year', 'month'],
value_vars=['%s-%s_rel' % (model, cs.upper()) for model in MODELS.itervalues()],
var_name='model-cs',
value_name='error'
)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='red')
sns.boxplot(x='model-cs', y='error', hue='month', data=data, ax=ax[n],
whis=[5,95], showmeans=True, meanline=True, meanprops=meanlineprops)
ax[n].set_title('Seasonal variation in %s by model' % cs.upper())
ax[n].set_ylabel('average monthly relative error')
ax[n].set_xlabel('model and clear sky component')
if n == 1:
ax[n].legend(bbox_to_anchor=(-0.1, -0.25), ncol=12, loc='lower center') # x = 0.5 for all 3 components
else:
legend = ax[n].legend()
legend.set_visible(False)
legend.remove()
f.tight_layout()
f.subplots_adjust(bottom=0.20)
f.savefig('seasonal-boxplot_by_model-cs_CORRECTED.png')
In [62]:
# plot
f, ax = plt.subplots(1, 2, figsize=(16, 6), sharex=False, sharey=False) # 24 x 6 for 3 subplots
for n, cs in enumerate(CS):
if cs == 'dhi': continue
if n == 2: n = 1
data = pd.melt(
MONTHLY_ATM_PARAMS_3MIN_CLEAR,
id_vars=['station', 'year', 'month'],
value_vars=['%s-%s_rel' % (model, cs.upper()) for model in MODELS.itervalues()],
var_name='model-cs',
value_name='error'
)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='red')
sns.boxplot(x='model-cs', y='error', hue='station', data=data, ax=ax[n],
whis=[5,95], showmeans=True, meanline=True, meanprops=meanlineprops)
ax[n].set_title('Regional variation in %s by model' % cs.upper())
ax[n].set_ylabel('average monthly relative error')
ax[n].set_xlabel('model and clear sky component')
if n == 1:
ax[n].legend(bbox_to_anchor=(-0.1, -0.25), ncol=12, loc='lower center') # x = 0.5 for all 3 components
else:
legend = ax[n].legend()
legend.set_visible(False)
legend.remove()
f.tight_layout()
f.subplots_adjust(bottom=0.20)
f.savefig('boxplot_by_model-cs_and_station_CORRECTED.png')
In [17]:
# plot all
STATION_ORDER = dict([(station_id, n) for n, station_id in enumerate(METADATA.index)])
f, ax = plt.subplots(7, 3, figsize=(24, 42), sharex=False)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
for n, cs in enumerate(CS):
data = pd.concat([
(station_atm_params_3min_clear['%s_%s' % (model, cs)] - station_atm_params_3min_clear[cs])[is_bright[station_id]]
for model in MODELS.itervalues()
], axis=1).resample('M').mean()
data.rename(columns=MODELS, inplace=True)
sns.boxplot(data=data, ax=ax[STATION_ORDER[station_id]][n])
sns.swarmplot(data=data, color=".25", ax=ax[STATION_ORDER[station_id]][n])
ax[STATION_ORDER[station_id]][n].set_title('%s by model at %s' % (cs.upper(), METADATA['station name'][station_id]))
f.tight_layout()
f.savefig('all_boxplot.png')
In [18]:
f, ax = plt.subplots(7, 3, figsize=(24, 42), sharex=False)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
for n, cs in enumerate(CS):
data = pd.concat([
(station_atm_params_3min_clear['%s_%s' % (model, cs)]
- station_atm_params_3min_clear[cs])[is_bright[station_id]]
for model in MODELS.itervalues()
] + [station_atm_params_3min_clear['year']], axis=1).resample('M').mean()
data.rename(columns=MODELS, inplace=True)
data = pd.melt(
data, id_vars=['year'], value_vars=MODELS.values(),
var_name='model', value_name='mbe'
)
sns.boxplot(x='model', y='mbe', data=data, hue='year', ax=ax[STATION_ORDER[station_id]][n])
#ax[STATION_ORDER[station_id]][n].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
ax[STATION_ORDER[station_id]][n].set_title('%s by model and year at %s' % (cs, METADATA['station name'][station_id]))
f.tight_layout()
f.savefig('all_boxplots_by_year.png')
In [19]:
f, ax = plt.subplots(7, 3, figsize=(24, 42), sharex=False)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():
for n, cs in enumerate(CS):
data = pd.concat([
(station_atm_params_3min_clear['%s_%s' % (model, cs)] - station_atm_params_3min_clear[cs])[is_bright[station_id]]
for model in MODELS.itervalues()
] + [station_atm_params_3min_clear['month']], axis=1).resample('M').mean()
data.rename(columns=MODELS, inplace=True)
data = pd.melt(
data, id_vars=['month'], value_vars=MODELS.values(),
var_name='model', value_name='mbe'
)
sns.boxplot(x='model', y='mbe', data=data, hue='month', ax=ax[STATION_ORDER[station_id]][n])
#ax[STATION_ORDER[station_id]][n].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
ax[STATION_ORDER[station_id]][n].set_title('%s by model and month at %s' % (cs.upper(), METADATA['station name'][station_id]))
f.tight_layout()
f.savefig('all_boxplot_by_months.png')
In [51]:
# plot all SURFRAD stations DHI by model
MODELS = {0: 'solis', 1: 'lt', 2: 'macc', 3: 'bird'}
# calculate station average DHI
avg_dhi = pd.concat([
pd.Series(data=station_atm_params_3min_clear['dhi'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# calculate DHI errors for each station and concatenate stations horizontally (axis=1) and average each station
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# calculate relative station DHI error by dividing average error by average station DHI
# then concatenate all models
data = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_dhi' % model] - station_atm_params_3min_clear['dhi'])[is_bright[station_id]],
name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_dhi for model in MODELS.itervalues()
], axis=1)
data.rename(columns=MODELS, inplace=True)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='black')
sns.boxplot(data=data, whis=[5, 95], showmeans=True, meanline=True, meanprops=meanlineprops)
sns.swarmplot(data=data, color=".25")
plt.title('All SURFRAD stations DHI by model')
plt.ylabel('average relative error by station')
plt.savefig('all_stations_all_years_by_model_DHI.png')
data
Out[51]:
In [52]:
# plot all SURFRAD stations DNI by model
MODELS = {0: 'solis', 1: 'lt', 2: 'macc', 3: 'bird'}
# calculate station average DNI
avg_dni = pd.concat([
pd.Series(data=station_atm_params_3min_clear['dni'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# calculate DNI errors for each station and concatenate stations horizontally (axis=1) and average each station
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# calculate relative station DNI error by dividing average error by average station DNI
# then concatenate all models
data = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_dni' % model] - station_atm_params_3min_clear['dni'])[is_bright[station_id]],
name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_dni for model in MODELS.itervalues()
], axis=1)
data.rename(columns=MODELS, inplace=True)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='black')
sns.boxplot(data=data, whis=[5, 95], showmeans=True, meanline=True, meanprops=meanlineprops)
sns.swarmplot(data=data, color=".25")
plt.title('All SURFRAD stations DNI by model')
plt.ylabel('average relative error by station')
plt.savefig('all_stations_all_years_by_model_DNI.png')
data
Out[52]:
In [50]:
# plot all SURFRAD stations GHI by model
MODELS = {0: 'solis', 1: 'lt', 2: 'macc', 3: 'bird'}
# calculate station average GHI
avg_ghi = pd.concat([
pd.Series(data=station_atm_params_3min_clear['ghi'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# sxf 560.271702
# gwn 591.976444
# fpk 543.374845
# tbl 608.619889
# bon 584.017623
# dra 654.731447
# psu 586.754301
# calculate GHI errors for each station and concatenate stations horizontally (axis=1) and average each station
# NOTE: averaging columns using `mean(axis=0)` transposes the dataframe so that stations are now the index
# sxf -10.020557
# gwn -2.605967
# fpk -7.293992
# tbl -24.014558
# bon -8.630023
# dra -29.584633
# psu -6.126202
# calculate relative station GHI error by dividing average error by average station GHI
# then concatenate all models
data = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_ghi' % model] - station_atm_params_3min_clear['ghi'])[is_bright[station_id]],
name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_ghi for model in MODELS.itervalues()
], axis=1)
data.rename(columns=MODELS, inplace=True)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='black')
sns.boxplot(data=data, whis=[5, 95], showmeans=True, meanline=True, meanprops=meanlineprops)
sns.swarmplot(data=data, color=".25")
plt.title('All SURFRAD stations GHI by model')
plt.ylabel('average relative error by station')
plt.savefig('all_stations_all_years_by_model_GHI.png')
data
Out[50]:
In [53]:
MODELS = {0: 'solis', 1: 'lt', 2: 'macc', 3: 'bird'}
# DNI
avg_dni = pd.concat([
pd.Series(data=station_atm_params_3min_clear['dni'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
dni = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_dni' % model] - station_atm_params_3min_clear['dni'])[is_bright[station_id]],
data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAgwAAAFiCAYAAACakIC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVGX7wPHvDDCsLpAi7qamuAuuKOaaqbkrZi6USZG5pJb7RlZvKmqaZUVqJaivGyii5pq45QJobri2uAsIiMCwDr8//DFv4wwyysCg3J/r4rqY5zznzH2GZe55VkVOTk4OQgghhBBPoDR3AEIIIYQo/iRhEEIIIUS+JGEQQgghRL4kYRBCCCFEviRhEEIIIUS+JGEQQgghRL4kYRCF6tixY9StW5dWrVqRkZGhd7xTp04MHz48z8d5ycjIICAggN69e9O0aVPc3d3p378/AQEBpKen69RdtmwZdevW5fjx4wavdfPmTerWrcvUqVP1znn8q2HDhnh6evLRRx/xzz//GHyeJ31FR0cDMHXqVIPHGzVqRMeOHZk+fTpxcXF53v/mzZupW7cuvXr1yrPO8OHD9a7fuHFjOnfuzJw5c7hz506+r/PTuH//Pqmpqc90bnJyMvHx8drHua/lzZs3TRWeyd24cYMFCxbQs2dP3N3dtb+D33//vcHXoW7duk/83c79nfj3Pef+DJ/k+PHj1K1bl+Dg4Ge/mWIqODj4iX+7pj5PPJmluQMQL7Zt27ZhZ2dHYmIi+/fvp1u3bgW+ZlZWFiNHjuT06dP07duXN998k+zsbCIiIli8eDH79+9n9erVqFSqAj/XBx98QM2aNbWP09LSOH36NCEhIURFRbFt2zbKli37xHP+rVKlSjqPp02bhqOjo/ZxcnIyv//+O5s3b+bcuXNs2rTJ4H3kvq6XL1/m7NmzNGrUKM97WLBggfb71NRUrly5wubNm9m5cyfr1q2jVq1aT34RjBAeHs4nn3xCSEgIdnZ2T3XuuXPnGDVqFAsXLqRVq1YAvPbaa1SrVg0nJ6cCx1YYdu/ezZQpU1CpVPTq1YtXXnmFrKwsTpw4wdKlSwkNDSUoKKjYxi/Es5CEQRSajIwMdu/eTZ8+fQgLCyMkJMQkCcPOnTs5ceIEy5Yto2vXrtpyb29vVqxYgb+/P5s2bWLIkCEFfq42bdpo38Ryvfnmm9SsWZOFCxeyceNG3nvvvXzPyUuXLl2oUqWKTtnQoUPx8/Nj3bp17N27lx49eugcj4mJ4fjx44wcOZIVK1YQEhLyxIShT58+emVeXl4MHjyYjz76iNDQUJTKgjU2njlzhqSkpGc69/Lly8TExOiUubq64urqWqCYCsvly5eZOHEirq6u/PjjjzoJ39ChQzl06BC+vr7Mnz+f+fPnmzFSIUxLuiREoQkPDycpKYlWrVrh6enJ4cOHiY2NLfB1T506BUDbtm31jg0ZMgQrKytOnz5d4Od5kn79+gHwxx9/FPn1d+zYgUajoWvXrjRq1Ijt27cb7O55knr16uHr68uVK1f47bffTBJzSfHll18C4O/vr5Ms5GrXrh09e/Zk+/btJCcnF3V4QhQaSRhEodm2bRsKhYIWLVrw2muvkZWVxdatWwt8XXt7ewDWr1+vd8zOzo6oqCidZvjCYGtrC0Bhraz+pOuHhobi4OBA/fr1ee2110hMTGTfvn1P/Ry54x8OHTr0xHq3b99m7NixeHp60qhRI3r06MGPP/6IRqMBHvW9f/PNNwB07txZp59+586dDBs2jGbNmtGwYUM6derEggULtAnOsmXLmDZtGvCohahTp07a8sf78xMSEvDz86Ndu3Y0bNiQ119/nYCAALKzs7V1li1bRqNGjfj777/x9fXFzc2NFi1aMGXKFBISEnTua926dfTq1YsmTZrQqlUrRo8ezZUrV574WsTGxnLs2DFee+01Xn755TzrjRkzhq1bt+Lg4PDE65lScnIy06ZNw93dnZYtWzJ58mRtgq7RaHj11VcZOHCg3nkHDx6kbt26HDhwwOB1c8dIHD16lBkzZtCiRQuaNWvGtGnTSE1NJTw8nD59+tCkSRP69OnD77//rnO+Wq1m0aJFdOrUSfs7sHDhQtRqtU69+/fvM23aNFq3bk2zZs2YPXu2wUQ4PT2dr776Snu9zp07s3Tp0qdOmsXTky4JUSiSk5M5cOAATZs2pVy5crRv3x6VSsWWLVvw8fEp0LV79+7NTz/9xPz58wkODqZLly54eHjg5uaGSqUyydiF/OS+ydarV0/v2MOHD3UG8OUqVaoUVlZWT3X9+vXr65T/9ddfnD9/np49e2JpaUmXLl1YtGgRISEhdO/e/anuoWrVqtja2nLx4sU862RmZuLj40NaWhrvvPMOpUuXJjw8nIULF5Kdnc0HH3zAm2++SXJyMnv27GHatGm88sorAGzcuJGZM2fSqVMnPvnkEzIzM9mzZw8rV67Ezs6OMWPG8NprrxEbG8v69ev54IMP8uxaefDgAYMHD+bWrVsMHjyYl19+mSNHjrBo0SIuXLjAkiVLtHU1Gg3e3t40b96cKVOmcPbsWTZt2kRaWhpLly4FHiVdfn5+9O3bl+HDhxMfH88vv/zC8OHD2bNnD6VKlTIYx8mTJ9FoNLRu3fqJr221atXyfD0N/W4ABX7D++qrr6hSpQrjxo3j3r17BAYGcvbsWUJCQrCxsaFHjx789NNP3Lx5U6cbbPv27ZQtW9Zgi92/TZ06ldq1a/Pxxx9z4sQJgoODuXv3LhcuXGD48OGUKlWKgIAAPvroI/bu3Uvp0qXJyMhgxIgRnD59mv79+9OwYUPOnDnDjz/+SGRkJKtXr8bKyor09HSGDRvGzZs38fb2pnz58oSEhLBjxw6dGLKzs/H19SUqKopBgwZRq1Ytzp07x/fff090dDTfffcdCoWiQK+jyJskDKJQ7Nq1i/T0dO0YAwcHB9q0acOBAwc4c+YMjRs3fuZrv/LKK3zzzTdMnz6dK1eucOXKFb777jvs7Ozo1KkTY8aMeeKnv6fx+Jt/amoqkZGRzJs3DycnJ4YNG6Z3zujRow1ea/Xq1XpjG5KSknSun5yczKFDh/jmm2+oVasWb7zxhk79bdu2AWhf15o1a/LKK69w+PBhYmJicHZ2fqr7K126NImJiXkej46O5tq1ayxdulQ7/sTLywsfHx/++usvANzc3Khbty579uzRGZOxatUq3NzcWL58ufaf+JAhQ+jcuTO7du1izJgxuLq60rRpU9avX//EsR8//vgjf//9N99++y1dunQBHo0X+PTTT1m7di39+vWjffv2wKNBsT169NDOehk8eDD37t1j7969qNVqbG1t2bZtG6+88orOGIN69eqxYMECLl++TLNmzQzGce/ePQBcXFz0jhlKBBwcHHQS2FOnTuHh4ZHXy10gTk5OrF+/XjvotE6dOkydOpWNGzcyfPhwevXqxU8//cTOnTu1424yMjLYu3cvPXv2zDeZdXZ2ZsWKFSiVSgYNGsSJEyc4evQoP/74I6+++irwqIVv5syZnD17lrZt27J582ZOnTrFtGnTeOedd4BHvwO1a9fG39+fDRs2MHToUDZu3Miff/6p8/MdNGgQXl5ePHz4UBvD1q1b+f3331mxYgXt2rXTljdu3JjZs2ezb98+7fnC9CRhEIUiLCwMeDTaPddrr73GgQMHCA4OLlDCANChQwd+++039u3bx4EDBzh69CixsbGEhYWxZ88eVqxYQcuWLZ/qmoY+mRh687eysqJNmzbMnj1bb4YEwJQpUwwO2DNUljtW4d9sbW3p3LkzM2fO1PsnHhYWho2NjfYfNDx6XZcvX87WrVv1BmDmJysr64mfyJydnVEoFPzwww/Y29vTqlUrVCoVK1euzPfaoaGhqNVqnevfv3+f0qVLP/X0y/3791OrVi29N4MPP/yQtWvXsm/fPm3CAOi1ttSrV49Dhw6RmJiIra0tLi4uHDlyhG+++Ya+fftSpUoV2rdvr3MNQ3K7YR7vKkpJSTGYCHz55Zf0799f+/jx6bv/tnLlSg4fPvzE53+SIUOG6MxQ6d27N19++SUHDhxg+PDhNGjQgJo1a+okDOHh4SQnJ9OzZ898r9+5c2ft4FilUknVqlV5+PChzu9ibrKY2xWyf/9+HBwcGDp0qM61vL29+e6779i/fz9Dhw7l4MGDlCtXTufna2dnh5eXl3bMCDyaneLk5ESDBg10ErT27dtjYWHBgQMHJGEoRJIwCJOLiYnh2LFj1KhRA4VCoe2HdnV1RaFQsGPHDqZPn17grgNra2t69OihnUVw/vx5Vq1aRVhYGHPmzGHnzp3aeoBOX/e/5ZYbiif3zT87O5uoqChWrlxJq1atWLBggcFkAaBBgwZGz5Lw9/enXLlyZGZmcujQIdasWUP37t3x8/PTxp3rzJkz/PPPP7Ru3Zr79+9ry3Ob8bds2fJUCUN2djZJSUlPbI1xcXFh0qRJLF68GB8fH+zs7PDw8KBHjx50794dCwuLPM+1srLi5MmThIWF8eeff3L9+nVt3JUrVzY6Tni0Vsa/P1HmKl++PKVLl+bWrVs65Y9PZ8z92eb+rEePHs3p06dZtmwZy5Yto3bt2nTq1AkvL688uxMAbQvO44N3bWxs+Omnn7SPL168aHCGRJkyZWjTpo3Ba4eGhub5vMZ4fCqvhYUFlStX1nltevbsyddff82NGzeoWrUq27dvp2LFijRv3jzf65crV07nsaWlpd7rnJtQ5CZWN2/epGrVqnqJr0qlomrVqtrYbt26RdWqVfWe8/HfzevXrxMfH59nK42p1xYRuiRhECaXO4r/77//pnPnznrHHzx4YHC6oDFSU1P54YcfaNCggc6USnj0Rr1o0SKSkpI4ePAgCQkJODo6Urp0aeDRp0BDHjx4AKCt9/g1c9/8cwfbjR49Gh8fH9asWaP3pv603N3dtZ/K2rdvT/Xq1fn8889JTEzUacqH/72hHDt2zODrevXq1afq7rl69SqZmZn5Tl8cOXIkPXv2ZM+ePYSHh3PkyBH27dvHli1bWLFiRZ7nLVq0iICAAOrXr0/Tpk3p06cPbm5ufPbZZ0/9j/1Jg0s1Go3eG1J+/dguLi5s3bqV48ePs2/fPg4dOkRAQAA//fQTq1atyrN1ys3NDXg0ENDLy0tbbmFhoZMIPCmRKiyG7jknJ0cnll69evH1119rB6MeOHCAt956y6h+f0P3lN95xv7cFAoFaWlp+Z6fnZ1NjRo1mDNnjsFrGvobFqYjCYMwudzZEfPmzdMbJX7x4kWWLVtGSEjIMyUM1tbWrFy5Ejc3N72EIVft2rU5dOgQNjY22scAV65c0ekiyXX58mUA7WC9J8mdBbB69Wr8/f2ZOXPmU9/DkwwfPpzff/+dffv28csvv2j7fbOzs9m5cyd2dnbMnz9fb92EI0eOsHbt2qfq7vn111+195SXxMRELl68iLu7O8OGDWPYsGGkpqYydepUdu3axaVLlwyuRHjr1i0CAgLo06eP3oyVJ61gmZfKlStrx0z8W2xsLMnJyVSsWPGprnfp0iUAPDw8tJ9WIyMjefvttwkMDMwzYahSpQru7u7s27dPb/CguT3eypKZmcnNmzd1Wg+qVatG48aN2b9/PzVr1kStVj9xtdCCqly5MqdPnyYzM1MnqcvIyNCJrUqVKkRERJCVlYWl5f/elm7cuKFzvSpVqnDu3Dlat26t8zeQO6DW0NgSYToyrVKY1F9//cW5c+do2bIlffv2pUuXLjpfvr6+lC9fniNHjmgHkD0NCwsLevTowYkTJwxO0UxMTGTXrl20adNGOzWxcePGlC9fno0bN+pNrcvIyGDdunXY2dnh6elpVAwff/wxVatWZc2aNYWy3sPcuXMpU6YMS5Ys0f7D/P3334mLi6Nbt2507dpV73UdO3YsKpXK6DUZrl69ys8//0yDBg2eOAjvyJEjvP322+zfv19bZmdnR506dYD/ferM/eed+4kwt9UmN1nLFR4ezt9//01WVpa27PFmbEM6duzItWvX2Lt3r055QEAA8GhMy9P46KOPmDx5sk43Vf369bGyssp3Eas5c+aQmZmpnY3wuJiYmCe2vBSW4OBgnfvZuHEjDx8+1OvT79WrF2fOnCE0NJSaNWvqzcQxpU6dOpGcnMyaNWt0yteuXUtKSor259a1a1cePnzIxo0btXUyMzPZsGGD3vUSExNZt26dTvl///tfJkyYoDelU5iWtDAIk8od7Ghovjc86tceMGAA33///TOvyTB16lTOnDnD5MmTCQ0NpV27djg4OHD9+nWCg4PJzMxk9uzZ2voqlQo/Pz/Gjx9Pv379GDhwIBUrVuT+/fts27aNq1ev8tlnnxm9jK+NjQ1+fn6MHDmSmTNnEhISYvR0SWOUK1eOTz75hFmzZuHn58fKlSvzfV2dnJzo2rUrYWFhet09/36dU1NTuXTpElu3bsXW1hZ/f/8nNit37NiRl19+mRkzZnD+/HmqVavGn3/+yZo1a/Dw8NAmBLmv3YoVK3j11Vdp164dlSpV4vvvvyc9PR0XFxfOnDlDSEgI1tbWOt1DueeuW7eOuLg4g594fX192b17N+PHj+ett96iRo0aHDt2jN27d9O1a9d8Bys+Lvdn984779CtWzdycnLYunUr6enp+a4Q6urqyrfffsukSZPo3r073bt3177p/vHHH+zevRu1Ws0bb7xBx44dnyqugrh58yZvv/02PXv25OrVq6xduxY3Nzf69u2rU69Hjx7MmzePXbt2MXbs2EKNycvLi5CQEObNm8fly5dp2LAh586dIzg4mKZNm2q7dfr06cOGDRv47LPPuHbtGjVq1CA0NFRvrEju9T777DPOnz9P48aNuXz5MuvXr6dBgwY6A0yF6UnCIEwqLCyMUqVK5dldAI+mSwUEBBASEvJMz+Hk5ERwcDA///wz+/bt49tvv0WtVuPs7EzXrl354IMP9KYXdunShTVr1vDTTz+xYcMG4uPjKVu2LA0bNmTWrFlPPaPC09OTXr16sW3bNn744QfGjBnzTPeSFy8vL7Zs2cLhw4fZsmULu3fv5uWXX85zuh/AW2+9RVhYGMHBwToJw+TJk7XflylTBhcXFwYMGMB7771HhQoVnhiHnZ0dq1at4uuvv2bbtm3ExcVRvnx5hgwZonPPb7zxBrt37yY4OJgTJ07QuXNnAgICmDdvHqtXryYnJ4dq1aoxffp0srKy+OKLLzh37hwNGzbEw8OD7t2789tvv3Hs2DGDvztly5Zl/fr1LFmyhB07dpCUlETVqlWZPHmyttvmaXh5eWFlZcXq1atZvHgxGo2Ghg0b8uOPPxo1YLV9+/Zs27aNjRs3snfvXnbt2kVmZiYuLi706dOHQYMG0aBBg6eOqyDmzJnD7t27mTdvHjY2NgwePJiJEyfqNPHDo4TUw8ODw4cPGzU7oiBUKhU///wz3377LTt37iQ0NBQXFxd8fX0ZNWqUNtG2sLBg5cqVLF68mJ07d5Kamsqrr77KO++8w4QJEwxeb9euXYSGhuLs7Mxbb73F6NGjta2KonAocgprqTohhBDFko+PDw8ePNDpAhAiPzKGQQghSpB//vmH33//XZrvxVOTFgYhhCgBDh48SEhICCdPngQeLYL0tFuRi5JNWhiEEKIEsLW15dChQ9jb2/P1119LsiCemrQwCCGEECJf0sIghBBCiHzJtMoniI19mH8lIYQQ4gVSvrzh7d2lhUEIIYQQ+SqWCUN2djaLFi3C09MTNzc3xo0b98T158+ePcvgwYNp0qQJXbt2ZcuWLTrH//nnHz788ENatWpF69atGTduHLdv3y7s2xBCCCFeGMUyYcjdnGj+/PkEBQVx9+7dPJcwjY+Px8fHhwYNGhAcHMzw4cOZMWOGdl/51NRURo4ciUaj4ZdffmHlypUkJCTw3nvvGbXmvhBCCCGK4RiGjIwMVq9ezcyZM2nbti0AixcvpnPnzkRFReHu7q5Tf+PGjTg4ODBjxgyUSiW1atXiwoULrFq1Ck9PT44cOcKdO3fYsmWLdufEBQsW0KFDB/744w9atGhR5PcohBBCPG+KXQvDxYsXSUlJ0Vnbv0qVKlSuXJmIiAi9+hEREbRo0UJnh7mWLVsSFRVFTk4OjRs3JiAgQGeb5dy6uTvqCSGEEOLJil0Lw927dwH0NsVxdnbWHnu8/uPbszo7O6NWq0lISKBChQp61woICMDOzk5nn3ghhBBC5K3YJQxqtRqlUqm3XbBKpSI9PV2vflpaGiqVSq8uYHCMwtq1awkKCmLWrFmULVv2ibE4OtphaWnxtLcghBBCvHCKXcJgY2ODRqMhKytLZ1vWjIwMg1uX2tjY6CUGuY8fr//dd9+xZMkSfH19GTZsWL6xJCSkPsstCCGEEM+tvNZhKHYJQ8WKFQGIjY3Vfg8QExOj17UA4OLiQmxsrE5ZTEwMdnZ2lCr16KY1Gg1+fn6sX7+eTz75hPfee68Q70AIIYR48RS7QY+urq7Y29tz4sQJbdnNmze5deuWwRkNzZo1IyIign9viXH8+HHc3d21gxvnzp3Lpk2b+PLLLyVZEEIIIZ5BsUsYVCoVQ4YMYcGCBRw8eJDz588zceJEWrZsSdOmTcnIyCA2Nlbb7TBw4EDi4+OZM2cO165dIzAwkLCwMHx8fAA4cOAA69atY9SoUbRr147Y2Fjtl6ExEUIIIYTQVyx3q8zKymLhwoWEhISQlZVFu3btmD17Nk5OThw/fhxvb29Wr15Nq1atADh9+jSff/45ly5dolKlSowbN4433ngDgI8//piwsDCDz7NgwQL69OmTZxyyl4QQwtQuXrwAgKtr/XxqCmEeeY1hKJYJQ3EhCYMQwpTU6lSmT/+EjIwMRo0aS8OGTcwdkhB6ZPMpIYQws9mzp3Lr1k1iY2OYO3cWx48fNXdIQhhNEgYhhCgCN25c559//tYp27t3l3mCEeIZSMIghBBFwNA6MobKhCiuJGEQQogiUK5ceVq0aKV9bG1tTd++A80YkRBPRwY9PoEMehRCmNqsWZNJT89g5sy5lC5d2tzhCKHnuVnpUQghXmRDh74DIMmCeO5IC8MTSAuDEEKIkkamVQohhBDimUnCIIQQQoh8ScIghBBCiHxJwiCEEEKIfEnCIIQQQoh8ScIghBBCiHxJwiCEEEKIfEnCIIQQQoh8ScIghBBCiHxJwiCEEEKIfEnCIIQQQoh8ScIghBBCiHxJwiCEEEKIfEnCIIQQQoh8WRpb8fr164SHh6NWq9FoNDrHFAoFvr6+Jg9OCCGEEMWDIicnJye/SqGhoUydOlUvUdBeRKEgOjra5MGZW2zsQ3OHIIQQQhSp8uVLGSw3KmHo1q0blStX5vPPP8fFxQWFQmHyAIsjSRiEEEKUNHklDEaNYbh16xY+Pj5UrFixxCQLQgghhPgfoxKGGjVqcPfu3cKORQghhBDFlFEJw4QJE/jmm284efIkWVlZhR0T2dnZLFq0CE9PT9zc3Bg3bhxxcXF51j979iyDBw+mSZMmdO3alS1bthisl5OTg4+PD8uXLy+s0IUQQogXklEJg7+/P/Hx8Xh7e9OoUSMaNmyo92VKy5YtIyQkhPnz5xMUFMTdu3cZO3aswbrx8fH4+PjQoEEDgoODGT58ODNmzODw4cM69TIyMpgxYwaHDh0yaaxCCCFESWDUtMo33nijsOPQysjIYPXq1cycOZO2bdsCsHjxYjp37kxUVBTu7u469Tdu3IiDgwMzZsxAqVRSq1YtLly4wKpVq/D09ATg/PnzzJgxg4cPH1K6dOkiuxchhBDiRWFUwjBmzJjCjkPr4sWLpKSk0LJlS21ZlSpVqFy5MhEREXoJQ0REBC1atECp/F9jScuWLfn000/JyclBoVDw+++/4+HhwejRo+ndu3eR3YsQQgjxojB64ab09HQ2b97MiRMnePjwIY6OjjRv3py+fftiY2NjsoByB1dWqFBBp9zZ2dngwMu7d+9Sv359vbpqtZqEhAScnJzw8fExWXxCCCFESWRUwpCYmIi3tzeXL1+mevXqvPTSS1y/fp2wsDACAwNZu3YtZcqUMUlAarUapVKJlZWVTrlKpSI9PV2vflpaGiqVSq8uPOreKAhHRzssLS0KdA0hhBDiRWBUwrBo0SLi4+PZsGEDjRs31pafOXOGDz/8kK+++go/Pz+TBGRjY4NGoyErKwtLy/+Fl5GRga2trcH6jycGuY8N1X8aCQmpBTpfCCGEeN4UaOGmffv28dFHH+kkCwCNGzfmo48+Yu/evQWP8P9VrFgRgNjYWJ3ymJgYvW4KABcXF4N17ezsKFXK8E0LIYQQ4ukYlTCkpaVp38gf5+LiQlJSkskCcnV1xd7enhMnTmjLbt68ya1bt2jRooVe/WbNmhEREcG/V7g+fvw47u7uOgMhhRBCCPHsjHpHrVOnDtu3bzd4LCwsjNq1a5ssIJVKxZAhQ1iwYAEHDx7k/PnzTJw4kZYtW9K0aVMyMjKIjY3VdjsMHDiQ+Ph45syZw7Vr1wgMDCQsLEwGOgohhBAmZNQYhlGjRuHr60tiYiJvvPEG5cqVIy4ujrCwMMLDw1m0aJFJgxo/fjxZWVlMmjSJrKws2rVrx+zZswE4deoU3t7erF69mlatWlGuXDlWrFjB559/Tt++falUqRLz58/Hw8PDpDEJUZxcvHgBAFfX+vnUFEII0zBqt0qAzZs389VXX+ks0VyuXDnGjx/PwIEDCy1Ac5LdKoW5JSYmsGXLZmJi7tK6dVtefbUjAPPnfwbAlCmzzBmeEOIFlNegR6PXYRgwYAD9+/fnzz//5MGDB5QpU4aaNWvK7pVCFJKcnBy++GIO//zzNwARESfIzMykcuUqXLoUDTxqaZBWBiFEUcgzYbh37x4vvfQSlpaW3Lt3T1vu4OCAg4MD8Gg2Qi5DMxiEEM/uxo1/tMlCrsOHw3UWStu6dbMkDEKIIpFnwtChQwfWr19P48aNad++fb4tCdHR0SYPToiSrHTpMlhYWJCdnf2vstJcuhRNcnIySqWS5ORkM0YohChJ8kwY/vOf/1C1alXt99L1IETRKlvWkYEDB7NhwzpycjQTnnbpAAAgAElEQVQ4OjphaakiPj4eAI1Gw61bN/UWORNCiMJg1KDH27dvU758eb3lmuHRHhPR0dE0bdq0UAI0Jxn0KApiw4Y1nDx5vMDXycrKJCsrC2trG+7evaO3smmlSpWwslLlcfbTa9GiFYMGDTXZ9YQQz5cCrfTYuXPnPLsczpw5w9tvv/3skQkhnsjS0gobG1sUCgUqlbXOMaVSiaWlfiIvhBCmlmcLw/z580lMTAQgJCSEDh064OjoqFcvOjqauLg4Dh8+XLiRmoG0MIjiJiUlmeXLv+bkyWNYWloyY8anNGjQyNxhCSFeIE89rfKVV17h+++/B0ChUHDx4kW9XSGVSiWlS5dm7ty5JgxVCJEXe3sHJk2aziefjEWhUEiyIIQoMnkmDP3796d///4AdOrUieXLl+Pq6lpkgQkh8iaDkIUQRc2oodX79+9/4vGUlBTs7e1NEpAQQgghih+jEoaMjAwCAwM5efIkmZmZ2p0hNRoNarWaS5cucfr06UINVAghhBDmY1TCsHDhQlavXk2dOnWIj4/H2toaJycnLl++TGZmJmPGjCnsOIUQQghhRkZNq9y1axcjRowgNDSUYcOG0bBhQzZu3Mju3bupXLkyGo2msOMUQgghhBkZlTDcv3+fV199FYA6depw9uxZ4NH+Ee+//z47duwovAiFEEIIYXZGJQylSpUiMzMTgOrVq3Pnzh3tGvY1atTgzp07hRehEEIIIczOqIShWbNmBAUFkZaWRvXq1bG1tWXv3r0A/PHHH9rdK4UQQgjxYjIqYRg9ejSRkZG8//77WFpaMmTIEGbPno2XlxdfffUVr7/+emHHKYQQQggzMmqWRL169dixYweXL18G4OOPP8bBwYGoqChGjRqFr69voQYphBBCCPMyKmE4efIk9evXp127dsCjVeY++OADAJKSkti7dy/du3cvvCiFEEIIYVZGdUl4e3tz7do1g8cuXLjAlClTTBqUEEIIIYqXPFsYpkyZop39kJOTg5+fn8HBjX///TflypUrvAiFEEIIYXZ5tjB0794dCwsLLCwsALTf//vLysqKZs2asWTJkiILWAghhBBFL88Whg4dOtChQwcAhg8fjp+fH7Vq1SqquIQQQghRjBg16DEwMFCv7MKFC9y5c4dWrVrJOgxCCCHEC86oQY8xMTG8/fbbLF++HICgoCAGDBjA6NGj6dq1K1evXi3UIIUQQghhXkYlDP7+/ly7do1GjRqh0Wj4/vvvadOmDVu2bKFmzZosXLiwsOMUQgghhBkZlTAcOXKEKVOm0K5dO6KiooiLi8Pb2xtXV1d8fHyIiIgwaVDZ2dksWrQIT09P3NzcGDduHHFxcXnWP3v2LIMHD6ZJkyZ07dqVLVu26BxXq9XMmjWLVq1a0bx5c2bOnElKSopJYxZCCCFeZEYlDCkpKVSsWBGAgwcPolKpaN26NQAqlYqcnByTBrVs2TJCQkKYP38+QUFB3L17l7FjxxqsGx8fj4+PDw0aNCA4OJjhw4czY8YMDh8+rK0ze/ZsIiMj+eGHH/j+++85ceIEs2fPNmnMQgghxIvMqIShRo0anDx5kszMTHbt2kXLli2xtrYGIDQ0lBo1apgsoIyMDFavXs3EiRNp27YtDRo0YPHixURFRREVFaVXf+PGjTg4ODBjxgxq1arF8OHD6d27N6tWrQLg7t27hIWFMWfOHJo2bUrz5s35/PPP2b59O/fu3TNZ3EIIIcSLzKiE4b333uObb77Bw8ODGzduMGLECAC8vLwIDQ3Fx8fHZAFdvHiRlJQUWrZsqS2rUqUKlStXNtj1ERERQYsWLVAq/3crLVu2JCoqipycHKKiolAqlbi7u2uPu7u7Y2FhQWRkpMniFkIIIV5kRk2r7NmzJxUrViQyMpKWLVvStGlTAFq1asWECRNo06aNyQK6e/cuABUqVNApd3Z21h57vH79+vX16qrVahISErh37x5OTk5YWVlpj1taWuLk5KRdyVIIIYQQT2ZUwgDQrFkzmjVrplP2ySefmDwgtVqNUqnUeYOHR2Ml0tPT9eqnpaWhUqn06sKj7g21Wq3tPjHmev/m6GiHpaXF096CEIXOwuJRi1r58qXMHIkQoqQwOmEoKjY2Nmg0GrKysrC0/F94GRkZ2NraGqyfkZGhU5b72NbW1uDx3Dp2dnZPjCUhIfVZbkGIQpedrQEgNvahmSMRQrxo8vogYtQYhqKUOxsjNjZWpzwmJkavmwLAxcXFYF07OztKlSqFi4sL8fHxZGdna49nZWURHx+Ps7NzIdyBEEII8eIpdgmDq6sr9vb2nDhxQlt28+ZNbt26RYsWLfTqN2vWjIiICJ2pncePH8fd3R2lUkmzZs3Iysri1KlT2uORkZFoNBq9LhYhhBBCGFbsEgaVSsWQIUNYsGABBw8e5Pz580ycOFE72DIjI4PY2FhtN8PAgQOJj49nzpw5XLt2jcDAQMLCwrQzNypUqED37t2ZMWMGkZGRREREMGvWLPr06WOwxUIIIYQQ+oxKGLy8vFi3bh1JSUmFHQ8A48ePp1evXkyaNAlvb28qVarE0qVLATh16hSenp7aFoNy5cqxYsUKLly4QN++fQkKCmL+/Pl4eHhor/f555/j7u7O+++/z+jRo2ndujV+fn5Fci9CCCHEi0CRY8QyjRMnTmTfvn0AdOzYkQEDBuDp6YlCoSj0AM1JBpSJ4mrSpHEA+Pt/beZIhBAvmrwGPRo1S2Lx4sUkJyezc+dOtmzZwvvvv0/58uXp3bs3/fr1o1atWiYNVgghhBDFi9FjGBwcHPDy8mLNmjXs3r2b4cOHc/ToUXr27Mmbb77Jpk2b8l3XQAghhBDPp6ce9JiZmUl0dDQXLlzg77//xtbWFkdHR+bNm0eXLl1MvnOlEEIIIczP6IWbIiIiCA0NZdeuXTx48IBmzZoxc+ZMunXrhp2dHSkpKbz77rtMmzaNPXv2FGbMQgghhChiRiUMnTp14s6dOzg7OzN48GD69+9P9erVderY29vTrl07AgMDCyVQIQrLf/7jR0JCvLnDeCq58eYOfnweODo6MX26n7nDEEI8I6MShkaNGjFnzhzatWunsyvk4/r168eAAQNMFpwQRSEhIZ778XEobYvdSul50igfTW5KUCeaORLjaNRZ5g5BCFFARv2HzF0DAeDOnTvExMTg6OhItWrVdOpVrlzZtNEJUUSUtpY4dquWf0XxTBJ+vW7uEIQwiYsXLwDg6lo/n5ovHqM/Uq1Zs4bvv/+euLg4bVmFChWYMGECffr0KZTghBBCiOJAo9Gwdu0vhIWFYmFhga/vaF59taO5wypSRiUMP//8M/PmzaNHjx507twZJycn4uLi2LlzJ1OnTkWhUNC7d+/CjlUIIYQwi4MHfyM0NAQAjSabb79dQp06rri4VDRzZEXHqIQhKCiIESNGMGXKFJ3yXr164efnx7fffisJgxBCiBfWpUvROo9zcnK4fPliiUoYjFqHITY2Fk9PT4PHXn/9de7evWvSoIQQQojipE4dVwNldc0QifkYlTC0aNEiz7UVTp48SePGjU0alBBCCFGctG/fiTZt2gGgUCjo188LF5dKZo6qaOXZJbFt2zbt9y1btuTrr78mLi6O119/nXLlyvHgwQMOHTrEjh07mD59epEEK4QQQpiDUqlk/PhJpKenAfDWW8PNHFHRy3O3SldX/eaXPC+iUBAdHZ1/xeeM7FZZMkyaNI4EdaJMqyxECb9ex9G2rOyuKZ57JWFa5VPvVpm7nbUQQgghHnmRE4X85JkwyCJMQgghhMj11LtVCiGEEKLkkYRBCCGEEPmShEEIIYQQ+TIqYbh582ZhxyGEEEKIYsyohMHb25utW7cWdixCCCGEKKaMShjS09NxdHQs7FiEEEIIUUwZtfnUuHHj+Pzzz/nwww+pU6cOL730kl6dChUqmDw4IYQQQhQPRiUMX3zxBZmZmUybNi3POi/iSo9CCCGEeMSohOHTTz8t7DiEEEIIUYwZlTD069evsOMQQgghRDFm9DoM8fHx+Pv74+XlRbdu3XjrrbdYtGgRcXFxJg3o/v37fPTRRzRv3hwPDw/8/f3Jysp64jmhoaG8/vrrNG7cmEGDBnHmzBmD9RITE/H09CQiIsKkMQshhBAvOqMShlu3btGnTx8CAwMpVaoUjRo1wtraml9++YW+ffty584dkwU0duxY4uLiCAoKYt68eQQHB7Ns2bI86x89epTp06fz7rvvEhISQp06dRg5ciTx8fE69WJjYxk5ciSxsbEmi1UIIYQoKYxKGPz9/bGzs2PXrl2sWrUKf39/fv75Z3bt2oWDgwMLFy40STCnTp0iMjKSefPm4erqSvv27Zk8eTKBgYFkZGQYPGflypX07NmTN998k1q1ajF37lzKlCnDhg0btHW2b99Onz59yGMnbyGEEELkw6iE4ejRo4wbN46KFSvqlFesWJExY8Zw5MgRkwQTERFB5cqVqVq1qrasZcuWpKSkGJyFodFoiIqKomXLltoypVJJixYtdLodfvvtN8aMGcPSpUtNEqcQQoiSKTExgdTUFHOHYRZGDXoEsLe3N1ju4OBAWlqaSYK5d+8ezs7OOmW5j+/cuUOTJk10jiUlJZGamqq3BoSzszNnz57VPs5tAZElroUQQjyLjIx0lixZSETEcSwtLenbdyCDBg0xd1hFyqiEoWHDhvz3v/+lQ4cOesfWrVtH/fr1jXqymzdv0rlzZ4PHVCoVvXv3xtraWqfcysoKhUJBenq63jm5iYqhcwzVf1qOjnZYWloU+DqieLOwkD3YioKFhZLy5UuZOwxRQqxatcpkrd8ADx484P79+wBkZWWxadN/+f33Q3rvPwXRtm1b3n33XZNdz9SMXulx6NCh9OnThx49elCuXDni4uLYsWMHV69eZcWKFUY9WYUKFdixY4fBY0qlkqCgIL2xCpmZmeTk5GBnZ6d3Tu4PytA5tra2RsX0JAkJqQW+hij+srM15g6hRMjO1hAb+9DcYYgSQq3OMOnftuEPrelYWlqZ7DnU6oxi8TeSV2JvVMLQtGlTAgICWLx4MUuWLCEnJweFQkGDBg0ICAjAw8PDqCCsrKyoVatWnsddXFwIDw/XKYuJiQEMLz1dtmxZ7OzstHX+fY4sVS2EECXXoEFDGTRoqMmuFxFxnAULvtA+trJS8eWXiyhbtuTss2RUW+zOnTupV68emzdvJioqivDwcCIjI9m0aRNt27Y1WTDNmjXjxo0bOtM0jx8/jr29Pa6urnr1FQoFbm5unDx5Ulum0Wg4efIkLVq0MFlcQgghSrbmzVvx/vsfolKpsLa2Yfr0OSUqWQAjE4ZPP/1U+6Zsa2tLhQoVDHYRFJSbmxtNmzZlwoQJnD9/nvDwcPz9/RkxYgQqlQqAlJQUnbUU3nnnHbZs2cKaNWu4du0as2fP5uHDhwwcONDk8QkhhCi5unTpRsWKlXBxcaFBg0bmDqfIGZUwlC1bFrVaXdixoFAo+Oabb3jppZcYOnQo06dPx8vLi9GjR2vrrFq1Ck9PT+3jV199lblz57Jq1Sr69evH1atXWbVqFU5OToUerxBCCFFSKHKMWM3ol19+YcmSJXTv3j3P7a179epVKAGaU3EYfCIK36RJ40hQJ+LYrZq5Q3lhJfx6HUfbsvj7f23uUEQx8p//+JGQEJ9/xWIkN15Hx+frQ6mjoxPTp/sZVbdAgx6//PJLAIKDgw0eVygUL2TCIIQQovAkJMRz//59rK1M38VdWBQ8mmqfnFT4re6mkp5pmhl/RiUM+/btM8mTCSGEEP9mbWWHe70B5g7jhRYVvdkk1zEqYfD39+ett96iVatWJnlSIcSz0WRkc//YLdQ3HmBZ2hqnVpWxcTa8CqsQQpiSUYMeDx06JBs3CVEMJJ6+S+rfieRk55CZkEZc+D/kaORvU4iicPPueY79sYGIc1uIf3DL3OEUOaMSBg8PD0JCQvLcMVIIUTTSY3Q3vclWZ5GVVPBl0IUQTxYT/xfRf4bzMCWOhKTbnIoOIz2jZK0GbFSXhIODA2FhYezevZuqVatSrlw5neMKhYKVK1cWSoBCiP9RlbcnI/5/m70pbSyxLKUyY0RClAxxCf/oPNZosol/cJOK5euYKaKiZ1TCcOvWLdzc3LSPMzMzCy0gIUTeHN1c0KgzSb2RhNX/j2FQyOZZQhQ6B1v9VR0d7J6vqZUFZVTCEBgYWNhxCCGMoFRZUL5DDXOHIUSJU9mlAfFJt4mN/wulwoIald0oZV8u/xNfIEYlDLnS09M5c+YMMTExeHp6olarcXFxKazYhBBCiGLBQmlJlQr10WiysLSwpkK52uYOqcgZnTCsWbOGpUuXkpSUhEKhYNOmTSxdupSMjAyWL19eKHtLCCGEEMVBwoPbnIrern0c/+AGbd2HYWVpbcaoipZRCcOmTZv4/PPP8fb2pmPHjrzzzjsADBw4kGnTprFs2TKmTJlSmHEKIYR4waSkpJCemWayhYUKU3Jqks7jzKx0Is5vwlplY6aIjJeemYoiRVPg6xg1WmrlypWMGDGCadOm6Wwb3bVrV8aPH8+uXbsKHIgQQghRXCkVFvplypI14NioFoabN2/q7BD5b3Xq1NHZbloIIYQwhr29PTnZyudiaejMrDQiz4fyMCUOgIrl69Lwlc5mjso4UdGbsbe3LfB1jEoYXFxcOHPmDG3atNE7Fh0dLQMfhRBCvNCsLG1o1diLpOQYLC1U2NvpT7N80RmVMAwYMIDly5djY2NDx44dAUhLS2Pfvn189913DB8+vFCDFEIIIcxNoVBQplQFc4dhNkYlDL6+vty+fZv58+czf/58AIYNGwZAjx49GDVqVOFFKEQROPLbISz/0B28VKXJyzTs5g7Ar/P1B2XJceOPH/ntED179NKrI4R4fhiVMCgUCubOncu7777LsWPHSExMpFSpUjRv3py6desWdoxCCCGEMDNFjmxDmafY2IfmDkEUgUmTxpGgTsSxWzVzh/LCSvj1Oo62ZfH3/9rcoYhiZNKkcSQnqZ+LQY/Ps6jozTiUtjX67698+VIGy0vWnBAhhBBCPBNJGIQQQgiRr6faS0IIIYQoiXJycrh24wQ3757H0kJF7eqtcCn3irnDKlLSwiCEEELk407sZf66GUlmVhrq9CTOXd6LOi0p/xNfIE/VwnD37l2OHTtGTEwM/fr1IzY2ltq1a6NSqQorPiGEEMLsEpNu6zzOIYfEh3ewtSltpoiKntEJw/z58wkMDCQrKwuFQkHbtm1ZvHgx9+7d45dffuGll14qzDiFEEIIsylTqgK3YqJ1yko7lKxFnIzqkggICCAwMJDJkyezZ88ecmdijhkzhgcPHvDVV18VapBCCCGEOVVydqWqSyOUSgusLG2pV6sD9rZlzR1WkTIqYVi/fj1jx47F29ubSpUqacvd3NwYP348Bw8eNFlA9+/f56OPPqJ58+Z4eHjg7+9PVlbWE88JDQ3l9ddfp3HjxgwaNIgzZ87oHD969Chvvvkmbm5udOzYkfnz55OWlmaymIUQQrzYFAolrjXb0anV+3RoOYIqFeqbO6QiZ1TCEBMTQ6NGjQweq1y5MomJiSYLaOzYscTFxREUFMS8efMIDg5m2bJledY/evQo06dP59133yUkJIQ6deowcuRI4uPjAbh48SLvv/8+Hh4ehISEMHfuXHbu3MncuXNNFrMQQoiSQaFQmDsEszEqYahWrRqHDh0yeCwiIoKqVauaJJhTp04RGRnJvHnzcHV1pX379kyePJnAwEAyMjIMnrNy5Up69uzJm2++Sa1atZg7dy5lypRhw4YNAGzatIl69eoxfvx4atSoQbt27Rg/fjyhoaFkZmaaJG4hhBDiRWfUoMe3336bOXPmkJWVRadOnVAoFNy4cYPIyEhWrlzJJ598YpJgIiIiqFy5sk4C0rJlS1JSUoiOjqZJkyY69TUaDVFRUcyaNUtbplQqadGiBREREQAMGjSIPn366JynVCrJzMxErVZjZWVlktiFEEKIF5lRCcOgQYNISEjgu+++IygoiJycHMaPH4+VlRXvvvsuQ4cONUkw9+7dw9nZWacs9/GdO3f0EoakpCRSU1OpUKGC3jlnz54FoE6dOjrHMjMz+fnnn2natCmlS5ec6TBCCCFEQRg9rdLX15ehQ4dy6tQp7W6VTZo0wdHR0egnu3nzJp07dzZ4TKVS0bt3b6ytrXXKraysUCgUpKen652TO3DR0DmG6mdnZzN16lSuXLnC2rVrjY5bCCGEKOmMShimTZtG79698fDwoF27ds/8ZBUqVGDHjh0GjymVSoKCgvTGKmRmZpKTk4OdnZ3eObmJgqFzbG1tdcrUajUTJ07k8OHDfP3113kO4vw3R0c7LC0t8q0nnm8WFrLgaVGwsFDmuQueKJnkb6/omOLvz6iEITIykpCQECpUqEDPnj3p06ePXlO/MaysrKhVq1aex11cXAgPD9cpi4mJAdDrdgAoW7YsdnZ22jr/Puff9RMSEvD19eXq1asEBATg4eFhVLwJCalG1RPPt+xsjblDKBGyszWyZbzQIX97Redp/v7ySiyMShh2797NmTNn2LFjB9u3b2fVqlW88sor9O3blzfeeMPgm/mzaNasGQsXLuTOnTtUrFgRgOPHj2Nvb4+rq6tefYVCgZubGydPnqRv377Ao4GQJ0+eZNCgQcCjbouRI0dy584dAgMDadCggUliFUIIUXDhh/YTFXlWp6xubXfate4NwIogP71z5PjTHX+YGssbb7yhV+9pGd0e1LhxY6ZOncqBAwdYvXo1zZs3Z9WqVXTq1IkRI0YUOBB4tBBU06ZNmTBhAufPnyc8PBx/f39GjBih3a8iJSWF2NhY7TnvvPMOW7ZsYc2aNVy7do3Zs2fz8OFDBg4cCMDSpUu5ePEi8+bNw9nZmdjYWO2XRiPZrRBCCGEMRU7uOs9PISUlhf3797N371727t2Lg4MDx48fN0lAsbGx+Pn5ceTIEezt7RkwYADjx49HqXyU2yxbtoxvvvmGS5cuac/ZvHkzy5cvJzY2lvr16zNr1ixtS4Knp6dOgvFv4eHhuLi4PCEWaT4tCSZNGkeCOhHHbtXMHYpRNJnZpMekYllahVUp6/xPKAYSfr2Oo21Z/P2/NncoohiZNGkcyUlq3OsNMHcoL7So6M04lLY1+u+vQF0SAOnp6ezfv58dO3Zw8OBBNBoN7du356uvvqJDhw7GXiZf5cuX59tvv83z+NixYxk7dqxO2YABAxgwwPAv3OHDh00WmxDmlhGv5t7uP9FkZANQ1s2FMo2c8zlLCCEKzqiE4eOPP+a3334jNTWVpk2bMmXKFHr06EHZsiVr4w3xYkpJSUGTnkXCr9fNHUq+UhOTtckCQOKpu2TfSEWhLN6jzTXqLFI0KeYOQwhRAEYlDGfOnGHEiBH07t2b6tWrF3ZMQog85Gj0exBzNDkoine+IIR4ARiVMOzZs6ew4xDCbOzt7clQZj4XYxgsL90n/vgt7WPr8na81P1lM0ZknIRfr2Nva2/uMIQQBZBnwjBr1ix8fX2pUqWKzl4NhigUCtn9UYgiUKruSyislKhvJGFZ2prS9cuZOyQhRAmRZ8Jw5MgR7R4RR44ceeJFSvJ2n0IUNYeajjjUNH5JdiGEMIU8E4b9+/cb/F4IIYQQJY9RQ6WmTZvGjRs3DB77888/GTVqlEmDEkIIIUTxkmcLw+3bt7Xfh4SE0KVLFyws9DdiOnjwYL5dFkIIIYR4vuWZMMydO1e7EZRCoWDMmDEG6+Xk5NCxY8fCiU4IIYQQxUKeCcOnn37K8ePHycnJYcqUKYwZM4Zq1XSnnVlYWFCqVClat25d6IEKIYQQwnzyTBgqVKhA796PdrvKXQbaycmpyAITQgghRPFh1MJN/fr1Iz09nXPnzpGZmUnuflUajQa1Wk1ERAQTJkwo1ECFEI82nrp/7BbqG0lYlVbh1Koy1uVlQaS+fXvolb32WjdGjx4nx010fMuWHXp1RMliVMJw4sQJxo8fT0JCgsHj9vb2kjAIUQQST90l9a9EADLi04g98A+VB9RDoZS1UIQQhcuo7a2HDBlCQkICEydOJDQ0FKVSSf/+/Tl48CDr1q1jzZo1uLm5FUW8RUq2ty4Znqftre+EXSYjPk2nrFLvOliVtTFTRMaR7a2FIbK9ddEw1fbWRq3DEB0dzZgxY3jttdfo2LEjd+7coX379syaNYsBAwbw3XffGR+5EOKZqcrpdj8obSywLKUyUzRCiJLEqIRBo9FQoUIFAKpXr86VK1e0x15//XUuXLhQONEJIXQ4urtgW7U0KMCytDXlX62OwkK2qhRCFD6jxjBUq1aNK1eu0Lx5c15++WXUajV//vknNWvWJDs7m5QU2edeiKKgVFng3LEGOTk5sofLC2L37p0cOLCP0qVLM3DgW9Su/Yq5QxLCIKMShp49e+Lv749Go2Ho0KE0bNiQL774Am9vb7777jtq165d2HEKE7t48VGrkKtrfTNHIp6FJAvF340b11m7djWxsfdo3bot/ft7oVTqrpZ79OghVqz4X5fuxYvRLF++Ajs7mfkiih+jEob33nuP+Ph4oqKiGDp0KHPmzOG9997D19cXBwcHGcPwnImLi2X58qVoNBomT55JtWrVzR2S2WnUWST8et3cYRhNk5ENPGpxeB5o1FlgWzjX/s9//EhIiDf5dVNSUsjISH+mc3NyckhOTtZOQb9+/R+2bt2MtbU12dnZaDQaLC0tUavVOuelpqYwatS7WFlZPXPcKpU19vamTzgcHZ2YPt3P5NcVzw+jEgalUsm0adO0jxs1asTevXu13RIODg6FFqAwrZSUZCZPHk9y8qMZINOmTWTBgqVUrlzFzJGZj6Pj87cgWULaozdIR9uyZo7ESLaF9zonJMQTfz8OB6Vpx3JkazTkO4Usr8K5+LEAACAASURBVHOzs3l8AlpWZiYajYbMzExtmZWl/r9gpUJBjkbzjM8M2WnqZ0508pJcgHjEi8OohMEQBwcHGjdubMpYRBGIiDihTRYAMjMzCQ/fz5Ah3maMyryex09NkyY9WlBHpilSaGOobAqQgGgUClIfK1Mplaj/lSwAKHNyUFlYkJH9qMXIwdqaUgaSiOJAxqqJPH8zGzRo8FT9pOfOnTNJQKJw2dra6ZXZ2emXieIpKyuLTZv+y+3bt7CysuLu3du4uFQyd1jiMUqFgjI2NiSlpZEDWCmV2KpUqLOydOrlAC/Z25Ol0aBUKFDK2BRRjOWZMHzwwQcysOoF5O7enGrVqnP9+j8AODm9RKdOXc0clTBWSMhGgoM3AI9ah7788jOWLFleov9W7e3tscpIZ1iZ4te1lKnRoNZoKP3/rQaBWbe5m5GhPd61jCP17It/l27Qg3hUhTAuQjxf8kwYxo4dW5RxiCJiaWnJggVLmDlzChqNhrlz56FSycI/hWHDhjWcPHncpNe8c+f2Y49vMWHChwUaJPe4Fi1aMWjQUJNdrySzUiqx+lfXhlf5CkQ+fMiDrEzq2tlTy0CLnxDFldGdZRqNhh07dnDkyBFiY2OZOXMmp0+fpmHDhjKt8jmjVFrg7f0ugCQLzxkrKysy/vUJVaFQYGHxfMyUEGCjtKBtmedkoKoQjzEqYXj48CE+Pj6cOXOGSpUqcfv2bVJSUti2bRtz584lKCiI+vVlPv/zRNZfKHyDBg01+Sf1uLhY5s//jH/++Rs7O3t8fD7A07O9SZ9DiKKUnplKVPRmc4dhtKzsRwm7pcXz82ErPTMVBxPMazYqYViwYAG3b98mJCSE2rVr07BhQwCWLl3KyJEjWbJkCQEBAQUOBuD+/fvMnTuXI0eOYGVlRf/+/ZkwYQKWTxg5HBoayrfffsudO3dwdXVl5syZOjM49u3bx7Jly/jzzz8pX748gwcPxsfHp0T3+4rnU7ly5fH3/5qYmHuULVsWlcra3CEJ8cyeyynNCY/WznAoXUgLixQCB2xN8loblTDs2bOHKVOm4OrqSvb/T/+BR1Mr33vvPWbMmFHgQHKNHTsWhUJBUFAQ9+7dY+rUqVhaWua5ffbRo0eZPn06s2bNonnz5vz000+MHDmSXbt24eTkxB9//MG4ceOYMGECr7/+OtHR0UyZMgVra2u8vUvuVELxfHN2rmDuEIQoMJnS/HwxaqJxWloaTk6GsxNra2udPtWCOHXqFJGRkcybNw9XV1fat2/P5MmTCQwMzPM5Vq5cyf+1d+dxUZZ7H8c/MzCsIoIiJoIeKcBdRMgs7eSSigtGiqYt5NYxV9LMUEnRoyL5PKmIdVxLPfiQqXjc8thiZUaAWpiCouYRRVlEQHaYef7wYR4JlDFnGJbf+/Xy9Wqu+75nvjeJ/Liu676uYcOGMWbMGFxdXQkNDcXW1pbo6HszyW/evMlrr73GpEmTcHZ25sUXX6R3796cPHlSL5mFEEKIxkCngqFz585ERUVVe+zQoUN6m78QHx+Pk5MTzs7O2jYfHx/y8/M5f/58lfPVajWnTp3Cx8dH26ZUKvH29iY+Ph64t5vm/PnzteefPHmSuLg4nnvuOb1kFkIIIRoDnYYkZs2axZtvvom/vz/PP/88CoWCw4cPs2HDBr755hs2bdqklzC3bt2iZcuWldoqXqelpdGtW7dKx3JzcykoKNBuvX3/NYmJiZXabt++TZ8+fSgrK6NPnz4EBAToJbMQQgjRGOhUMHh7e7N161ZWr17NJ598gkajYfPmzXTo0IENGzbwzDPP6PRhqamp9O/fv9pjZmZmjBgxAnPzypO4VCoVCoWC4uKqa6MXFRUBVHvNH8+3sLAgOjqa1NRUli5dyvz581m9evVD89rZWWFqKo+sCVGXmZjodw8JUT0TEyUODjbGjmF0FX/fGuPXQqeC4dSpU3Tv3p1du3ZRVFRETk4OTZo0eeQd0RwdHTl06FC1x5RKJTt27KgyV6G0tBSNRlPt8sUVhUJ111haVp7BamVlRadOnejUqRPl5eUEBQUxb968Kr0T98vO/uNq8EKIuqa8XDZGqg3l5WoyMvJqPrGBq/j71pC/Fg8qhnQekpgzZw4jR47EwsICCwuLPxVCpVLh6ur6wOOtWrXi+PHjldrS09MBqv3B3qxZM6ysrLTn3H9NxfmJiYmUlJTg5eWlPe7m5gbcGwJ5WMEghBBCiHt06sszNTWtlS2svby8uHbtGmlpadq22NhYrK2t8fDwqHK+QqHA09OTuLg4bZtarSYuLg5vb28AvvjiCxYvXlxpq9lff/0VlUpFu3btDHczQgghRAOiUw/D1KlTCQkJITk5GTc3N5o3b17lnB49ejx2GE9PT7p3705QUBCLFi0iMzOT8PBw3nzzTe0Sxvn5+RQUFODg4ABAYGAgU6dOpWPHjvTq1YutW7eSl5fHqFGjABg3bhxffPEFK1eu5JVXXiE5OZlVq1YRGBhI06ZNHzuzEEII0RjoVDCEhIQAsG7dOoBKKyRqNBoUCkW1jz0+KoVCQUREBIsXL2b8+PFYW1szevRopk2bpj1ny5YtREREkJycDEDfvn0JDQ0lMjKSsLAwOnbsyJYtW7TrRri5ubF582ZWr17Nrl27sLe3Z8KECUyePPmx8wohhBCNhUJzf1/9A/z88881vtH9ayE0FA15UosQDcW7786kJPt2ndzeuqHYkXMbMzv7Rrm64R81hpUeH2vSY0MsBoQQQgihO523txYNx82baRw5coCSkhL6938RV9enjB1JCCFEHScFQyOTl5fLggXvkpeXC8C3337FypX/hYtLO+MGE0IIUadJwVDHRUfvJC4uVm/vd/dunrZYACgrKyM0dKFet5n19n6agIDxens/IYQQxidrqjYySmXVpa6raxNCCCHu98g9DDdu3CA9PR03NzcUCkWVJZiFfgUEjNfrb+tqdTl///tiEhN/AaBNG2dCQ8NqZWEuIYQQ9ZfOPQxff/01gwcPpn///owbN44rV64wd+5cgoODKS8vN2RGoUdKpQkLF4bi6NiKli0dCQ9fK8WCEEKIGulUMHz99ddMmzaNJ598kmXLlqFW39t8o3fv3uzfv59PPvnEoCGFfikUCiwsLLC0tMTERIYjhBBC1EyngmHdunW89NJLREREMHLkSG37+PHjmTZtGjExMQYLKIQQQgjj06lguHTpEr6+vtUe8/LyqrRZlBBCiAcrUpeTVVqCDovsClGn6DTp0c7Ojt9//53nnnuuyrHff/8dOzs7vQcTQoiG5nReLt/mZFOm0dBCpWJUC0dsTOXpdlE/6NTD4Ovry5o1azh27BilpaXAvXHwpKQkIiMjGTx4sEFDCiFEfVdQXs43d25T9n89C5mlpfyYe8fIqYTQnU6l7ezZs7l48SLTp0/H9P+q4cDAQPLy8vD09GTWrFkGDSmEEPVdXnkZf3yeLLus1ChZhPgzdCoYzM3N2bRpEydOnODkyZPk5OTQpEkTfHx8+Otf/1ppu2shhBBVOajMsDUxJae8TNv2lKW1ERMJ8WgeafDs2Wef5dlnnzVUFiGE+FPuqtXsyLlt7Bg1MrO0QFVUhFqtxkql4lx5KefrQe67ajWyebjQqWB4//33H3hMqVRiZWVFu3bt8PX1lQmQQohapc99UAxNXVREeUEBarUajUqFqpkdSmXdX6Hfnvr1dRaGoVPBcPPmTU6dOkVxcTFOTk44ODiQlZVFamoqSqWSFi1akJWVRWRkJFFRUbi4uBg6d52yfPlisrPr/m8J96vI++67M42c5NHY2dkTHLzY2DFEHVJf/j6UlBTzt79N0C58d/fuXV580ZexY181cjIhdKNTwfD8889z6dIlIiIi6Nq1q7Y9KSmJGTNmMHHiRIYOHcrUqVNZvXo1a9asMVjguig7+zZZWVkoVPVnXw3N/z0gczu3wMhJdKcpLTR2BNGI6Hun2OLiYu7ezavUdvBgDAkJP+vtM0B2ixWGo1PBsG3bNubMmVOpWADw8PBg9uzZfPjhh4wdO5bAwEAWLlxokKB1nUJlSZMnRxg7RoN2N2W/sSMI8aepVCoUCkWlBZvMzc2NmEiIR6NTwZCTk4ONjU21x8zNzcnOzgbA1taW4uJi/aUTQggj0fdOsQCnTyewadMGMjMzefrpXkydOhNLSyu9foYQhqJTweDp6UlERASenp6VJjXm5OTw8ccfa3seTp8+TZs2bQyTVAgh6jlPTy/Wr9+EWl2OUikbv4n6ReenJF599VX69etHz549sbe3Jysri1OnTmFubs62bdv48ccf+eijj1iwYIGhMwshRL0mxYKoj3R6nuepp57i8OHDBAYGkp+fzy+//EJpaSkTJkzgyJEjuLu706RJEz788EPGj5fJNkIIIURDo/PCTfb29g9dArpr165VJkUKIYQQomHQuWBITEwkLi6O0tJS7SxftVpNYWEh8fHxREVFGSykEEIIIYxLp4IhKiqK0NDQavdvVyqV1W57LYQwnKSkcwB4eHQ0chIhRGOh0xyG7du307dvX2JjY5kwYQIBAQGcOXOGNWvWYG5uzogR+lt/ICsri1mzZtGzZ0+eeeYZwsPDKSsre+g1+/fvZ9CgQXTt2pWAgAB+/fXXB567ZMkS+vXrp7e89YFGo6m22BP1V3T0TvbsiTZ2DCEaldLSUgoLCygqKjJ2FKPQqWC4du0a48aNw9bWls6dO5OQkICFhQWDBg1iypQpfPbZZ3oLNGPGDDIzM9mxYwcrV65kz549rFu37oHn//jjjwQHBzNhwgT27t2Lm5sbEydO5Pbtqks1f//99/zzn//UW9a6TqPRkPOfk9xM2MatMzvITz+nPVZWnEdpYbYR04k/IyMjnRkzpnD2bCJnzpzif/5np7EjCdEo5OTc4Z13ppGens6tWzcJC1vW6H4R02lIQqVSYWFhAUDbtm25evUqpaWlqFQqvLy82Lp1q17CnD59moSEBI4dO4azszMeHh7MmzePpUuXMm3aNMzMzKpcs3nzZoYNG8aYMWMACA0N5aeffiI6Opq//e1v2vPu3LlDcHAwPj4+XL9+XS9567rCrIsU3PoNAE1ZOblXf0TVpBUFtxIpzLwIgJlNa+yfehGFySNtXCqMJCpqO7du3dS+3rv3c/z8/LGwqD/LkgtRG/S9tPedO9nk5ORoXyck/MzMmVP0+r1X15f11umnhIeHB99++y1PP/00f/nLX1Cr1fzyyy/07NmTW7du6S1MfHw8Tk5OODs7a9t8fHzIz8/n/PnzdOvWrdL5arWaU6dOsWjRIm2bUqnE29ub+Pj4Sud+8MEH9O/fn5YtW7J79269ZQbIz89HU1pU55YuLirIq9KWm3KY0uL/35OhJO8G2ed3Y2ZR91eb05QWkp/fuCr6P7pxo3Kxq1aruX07i9atZcG0+kLmn9RPFZuG3a+8vGpbQ6ZTwfDGG28wa9Ys8vLyWLZsGf3792fevHkMGTKEmJgYvLy89BLm1q1btGzZslJbxeu0tLQqBUNubi4FBQU4OjpWuSYxMVH7OiYmhnPnzhETE8O2bdv0krU+MDFVVSoOAFBUHYVSq8trKZF4XF5e3ly+nKJ93aKFA0884WTEROJRRUVtp7i4mJCQZVhZ1f1Cvb7S99LeKSkXWbRoHuXl9/69tLOz57/+K6JR9e7pVDAMGjSI9evXc/nyZeBet/+cOXPYuXMnXbp0ISQkRKcPS01NpX///tUeMzMzY8SIEVU2Y6nYsKW6PSoqJp5Ud03F+WlpaSxfvpz169c/8jennZ0VpqY1r8jWtKkNxeWKOrn5VN71UxSkn0NhYoaNUw/MbJ4gI/FzNOqKiaQKmrYfgFkTx4e+T11wN2U/TZs2wcGh+n1NGoO33pqIpaWK3bt3Y2FhwUcf/TctWzY1diyho9DQUM6fvzdM+PbbE1i7di2urq5GTiV04eDQgzVr1nD48GGsrKzw9/enVauWNV/YgOhUMBw+fJinn35a+8Pezs6OLVu2PPKHOTo6cujQoWqPKZVKduzYQUlJSaX2inUfqvthX1EoVHeNpaUlGo2G+fPn4+/vT8+ePR85b3a2bls/1+VuKRunHtg49ajUZu/uS/7NX9Goy7Bq2bFeFAsVysvVZGRUHWppTIYNG8WTT97rzjYza9rovx71RVZWJt988432dUFBAVu2fMrs2e8aMZV4FC1buvDGG29pXzfU770H/VKmU8GwZMkSlixZwqBBgx4rhEqlemg13apVK44fP16pLT09HaDKsANAs2bNsLKy0p5z/zWOjo7cuHGDn376iTNnzrBr1y7gXjFRVlaGp6cnGzdu/FOFRH1n1qQlZk8OMHYM8Rhk/Lv+yc/Pr9KWl5drhCRC/Dk6PVbZrFkzCgsLaz7xMXl5eXHt2jXS0tK0bbGxsVhbW+Ph4VHlfIVCgaenJ3Fxcdo2tVpNXFwc3t7eODo6cvToUfbv38++ffvYt28f48ePp2XLluzbt4/OnTsb/J6EEALAxaUtTk7Oldr69RtopDRCPDqdehheeeUVlixZws8//4ybmxvNmzevcs7w4cMfO4ynpyfdu3cnKCiIRYsWkZmZSXh4OG+++ab2kcr8/HwKCgpwcHAAIDAwkKlTp9KxY0d69erF1q1bycvLY9SoUZiamtK2bdtKn2Fra1tt++P68ftvMImvvGBU6/aedPAZCsBXu5ZVuUaOP9rx8sIshg0dWuU8IeqLpUvDWLBgLiUlJUyePJUePbyNHUkInelUMKxYsQKAPXv2VHtcoVDopWBQKBRERESwePFixo8fj7W1NaNHj2batGnac7Zs2UJERATJyckA9O3bl9DQUCIjIwkLC6Njx45s2bIFe3v7x84jhBD61KRJE6ZOnQnIsJKofxQaHZaq0mWhIyenhvdol64TWt59dya3cwvq5FMSDcndlP3YN7UiPHytsaMIIUSD9ViTHhtiMSCEEEII3em8HnBaWhobNmzgxIkTZGRkEBUVxYEDB3B3d2fkyJGGzCiEEEIII9PpKYlLly4xcuRIvv32W3x8fCgtLQXg7t27vP/++xw+fNigIYUQQghhXDpPemzfvj2ffvopSqWSvXv3ArB06VKKi4vZtGkTQ4YMMWhQIYQQQhiPTj0MCQkJTJo0CTMzMxQKRaVjL730knbJaCGEEEI0TDoVDCqVqsryyxVyc3Or3XZaCCGEEA2HTgVD7969WbduXaUlmBUKBUVFRWzdupVevXoZLKAQQgghjE+nOQzz5s1j7NixDBo0iE6dOqFQKAgPD+fKlSuUlJSwatUqQ+cUQgghhBHp1MPQunVrYmJieP311yktLcXFxYXc3FyGDBnC3r17cXFxMXROIYQQQhiRTj0MqamptGnThqCgIEPnEUIIIeqspKRzQONc2lunHoYBAwYwbtw4oqOjyctrmPt/CyGEEA+TnHyeNWs+JCLiv8nISK/5ggZGp4IhLCwMGxsbQkNDefbZZ5kxYwbHjh3TLuAkhBBCNGQpKRf44IP3ycrKJD39FvPnB1FUVGjsWLVKpyEJPz8//Pz8yMnJ4csvv+TgwYPMnDkTGxsbBg8ezIgRI/Dy8jJ0ViGEEMIovvvuG9RqtfZ1Xl4ep08n8MwzzxkxVe3SeS8JAFtbWwICAggICCAzM5OPP/6YqKgooqOjOX/+vKEyCj0rK8oh/9ZvaNRlWDl4YNakpbEjCSFEndakSdUdHG1smhohifE8UsEAkJyczKFDhzhy5AhXr17lqaeews/PzxDZhAGoS4vIPLcfTXkxAIVZKbTo6IfKqrmRkwkhRN01aJAvX3/9b27fzgLA3b0DnTp1MXKq2qVTwfD7779z8OBBDh8+zKVLl2jevDnDhw/Hz88PDw8PQ2cUelR05z/aYgEAjZrCrEtSMAghxEPY2jYjIuIfhIS8j4mJCUuXhhk7Uq3TqWAYPHgwlpaWDBgwgPnz59O7d2+USp3mS4o6Rqmy0KlNCCFEZaamKl5/fYKxYxiNTgXDypUrefHFF7GysjJ0HmFg5rZtMLdtQ3FOKgCmlnZYtXA3ciohhKgfGuP6CxV0KhhGjhwJQHZ2NqWlpWg0GgA0Gg0FBQUkJCQwevRow6UUeqNQKLF3G0zJ3XQ06jLMbFqhUEhvkRBCiIfTqWBITk5m7ty5pKSkVHtcoVBIwVDPyJMRQgghHoVOBcOqVau4c+cO7733Ht988w1mZma88MILfPfdd3z33Xd89tlnhs4phBBCCCPSqS/6zJkzzJo1i8DAQHx9fSksLGTcuHF8/PHHDBgwgO3btxs6pxBCCCGMSKcehpKSEtq1awdAu3btSEpK0h7z9/fngw8+MEi4+kRTWsjdlP3GjqEzTXkJAAoTMyMn0Z2mtBCQibdCCGEMOhUMrVu3JjU1lZ49e9KuXTvu3r3L9evXcXJywtzcnJycHEPnrNPs7OyNHeGRZWcXAWDXtD79ALaql19rIYRoCHQqGAYMGMCHH36ItbU1AwcOpH379qxZs4a33nqLbdu24ezsbOicdVpw8GJjR3hk7747E4Dw8LVGTiKEEKI+0GkOw/Tp0+nevTvR0dEAvP/++3z55ZcMGzaMEydOMGPGDIOGFEIIIYRx6dTDYGlpSUREBCUl98a9+/Tpw4EDBzh79iydOnXCxcVFb4GysrIIDQ3lxIkTqFQq/P39CQoKwtT0wVH379/P+vXrSUtLw8PDg4ULF9K1a1ft8VWrVrF58+ZK17i4uPDvf/9bb7mFEEKIhuyRNp8yM/v/CXLOzs4GGYqYMWMGCoWCHTt2cOvWLebPn4+pqSlBQUHVnv/jjz8SHBzMokWL6NmzJ1u3bmXixIl8+eWX2NvfG+++cOEC48ePZ+rUqdrrTExM9J5dCCGEaKjq1BJ/p0+fJiEhgZUrV+Lh4cHzzz/PvHnz2L59u7Z34482b97MsGHDGDNmDK6uroSGhmJra6sdPgG4ePEinTp1wsHBQfunophojNTqcgoLCykoyH/g11UIIYS4X50qGOLj43FycqrUc+Hj40N+fj7nz5+vcr5arebUqVP4+Pho25RKJd7e3sTHxwOQl5fHzZs3cXV1NfwN1ANlZWUsWbKA9PRbZGRkMHfuDHJzc40dSwghRB33SEMShnbr1i1atqy8ZHHF67S0NLp161bpWG5uLgUFBTg6Ola5JjExEbg3HAGwZ88e5syZA0Dfvn155513sLGxMch96FN09E7i4mL19n4FBQVkZKRrX9+8mcbcuTOwtbXV22d4ez9NQMB4vb2fEEII46vVgiE1NZX+/ftXe8zMzIwRI0Zgbm5eqV2lUqFQKCguLq5yTVHRvbUEqrum4vyK/S+aNWtGZGQkqamphIWFkZKSwmeffYZCoXhgXjs7K0xNjTvXwdLSDBMTfXYEaapt0+dnWFqa4eBQ94sxIYQQuqvVgsHR0ZFDhw5Ve0ypVLJjx44qY+oVu2NWt7V2RaFQ3TWWlpYABAQEMHDgQO2cBXd3d1q0aEFAQAC//fYbnTt3fmDe7OwC3W/OQIYPH83w4frb2Cs//y5BQdO4cycbuFeohYT8nTZt9DuBNSMjT6/vJ4QQonY86Be+Wi0YVCrVQ+cStGrViuPHj1dqS0+/133+x2EHuNdrYGVlpT3n/msqzlcoFFUmOLq5uQFw8+bNhxYMDZG1dROWL/+Qo0cPU1JSTL9+A/VeLAghhGh46tSkRy8vL65du0ZaWpq2LTY2Fmtrazw8PKqcr1Ao8PT0JC4uTtumVquJi4vD29sbgLCwMPz9/Stdd/bsWYBGOxGyRQsHxo17ncDAybi4tDN2HCGEEPVAnSoYPD096d69O0FBQfz2228cP36c8PBw3nzzTe0aEPn5+WRkZGivCQwMZN++fezcuZNLly4REhJCXl4eo0aNAmDgwIEkJSWxatUqrl69yg8//EBwcDDDhw/nL3/5i1HuUwghhKhvFBqNprpZcEaTkZHB4sWLOXHiBNbW1rz88svMnj0bpfJebbNu3ToiIiJITk7WXvPFF18QGRlJRkYGHTt2ZNGiRXTq1El7/Pjx46xbt46UlBSsra0ZNmwY77zzTpXJklWzyDi8EEKIxuVBcxjqXMFQl0jBIIQQorF5UMFQp4YkRO1JSjpHUtI5Y8cQQghRT9SphZtE7bhy5RIffbQKtVrN9OlBdO3qaexIQggh6jgZkniIhjgkcedONtOnT9auXaFUKlm2bBVPPulm5GRCCCHqAhmSEAAkJMRVWuhKrVZz8uQJIyYSQghRH0jB0Mg0b95CpzYhhBDiflIwNDLdunnSpcv/b+Ll4tKWfv0GGDGREEKI+kDmMDxEQ5zDUGHJkmDUajVLlqw0dhQhhBB1SJ3YS0LUHWPGvGrsCEIIIeoR6WF4iIbcwyCEEEJUR56SEEIIIcSfJgWDEEIIIWokBYMQQgghaiQFgxBCCCFqJAWDEEIIIWokBYMQQgghaiQFgxBCCCFqJAWDEEIIIWokBYMQQgghaiQFgxBCCCFqJAWDEEIIIWoke0kIIYQQokbSwyCEEEKIGknBIIQQQogaScEghBBCiBpJwSCEEEKIGknBIIQQQogaScEghBBCiBpJwdCIxcbG4u7uzs2bNwHo168fkZGRRk4lHsbd3Z2YmBgACgsL2blzp5ETCVH/3f99VZ358+cTGBj4SO8ZExODu7v7YyarW0yNHUDUHbt378bCwsLYMYSOtm3bxueff8748eONHUWIBm3BggWo1WpjxzA6KRiElr29vbEjiEcga64JUTtsbGyMHaFOkCGJBmDPnj0MGTKEzp0788ILL7B27VptNfzVV1/h7+9Pt27d+Otf/8q6desoKyur9n3uH5LIzMxk+vTp+Pj40L17dwIDAzl//nyt3ZN4uD179rBmzRquX7+Ou7s7sbGxxo7UaLi7u3PgwAFeeeUVunTpgq+vL2fOnOGf//wnzz//PD16JiqMEAAABuhJREFU9OCdd96hpKREe01UVBTDhg2jS5cueHp6MmHCBK5evao9npmZyZw5c/Dx8cHb25uZM2eSnp6uPb5v3z6GDx9O165dGTRoEHv37q3Ve24MUlJSGD16NJ07d8bPz48zZ85oj90/JBEbG0uXLl2IjIzEx8eH1157DYCTJ0/i7+9P165dGTNmDKmpqca4DYOSgqGeS0pKIiQkhKCgII4ePUpwcDCbN29m//79HD16lBkzZjBkyBBiYmKYN28e27dvZ8WKFTW+75IlSygrKyMqKoo9e/ZgbW3NjBkzauGOhC58fX2ZPHkyrVq14ocffsDT09PYkRqVFStWMHnyZGJiYmjSpAlTpkzhq6++YuPGjaxYsYKjR4+ye/duAI4cOcKKFSt4++23OXLkCJ988gnXr18nLCwMgLKyMiZMmEBqair/+Mc/2LFjB5mZmcycOROAQ4cOsWDBAkaNGsW//vUvJk2axMKFC/nhhx+Mdv8N0WeffcbYsWOJiYnBy8uL119/XTu/649KSkqIjY3l888/Z+HChVy9epUpU6bQo0cP9u3bx9ixY9m4cWMt34HhyZBEPXft2jUUCgWtW7fW/tm6dSutWrVi5syZDBkyhMmTJwPQrl077ty5w9///ndmz5790Pe9evUq7u7utGnTBnNzc0JDQ0lJSUGtVqNUSp1pbBYWFlhZWWFiYoKDg4Ox4zQ6o0aNol+/fgD4+fkRGhrK4sWLcXZ2xs3NjU2bNnHx4kXg3lDf8uXL8fX1BcDJyYmhQ4eyf/9+4N5vpsnJyRw7dgxnZ2cAli1bxp49eyguLubTTz9l+PDhvPHGGwC0bduW/Px8GVPXs9dee42XX34ZgIULF/L9998TFRVFUFBQtedPmjSJtm3bAhAeHs4TTzxBcHAwSqWS9u3bc/HiRTZv3lxr+WuDFAz1XJ8+fejWrRsvv/wybdu25bnnnsPX15fWrVtz8eJFRo4cWel8b29vysrKuHz58kPf9+233+a9997j6NGjeHt707dvX0aOHCnFghCAi4uL9r8tLS1RKpW0adNG22ZhYaEdkvDx8eHChQtERERw+fJlrly5woULF3B0dATgwoUL2Nvba4sFgPbt2zN37lzt8REjRlT6/EedsS9qdn8vnVKppGPHjtqirzr3//+6ePEiHTp0qPTvY/fu3Q0T1IjkX/96zsLCgh07drB79278/Pw4d+4cr776Khs3bqz2iYfy8nIATE0fXisOHjyY77//nmXLluHg4EBkZCQjR44kMzPTIPchRH3yx+8fhUKBQqGo9tyYmBj8/f25ceMGPXv2ZNGiRdpev+req6bPEoZhYmJS6bVGo8HMzOyB59//76tCoagyCVmlUuk3YB0gBUM9d+LECdavX0+XLl2YNm0au3btYuzYsezduxdXV1dOnTpV6fyEhARUKlWl35D+qKysjLCwMK5fv87w4cNZsWIFBw8e5Pr16/z888+GviWhowf9gBJ1S8XY+PLlyxk3bhw9evTgP//5j/YHjKurK7dv3+b69evaay5dukSvXr1ITU3F1dWVs2fPVnrPefPmsWzZslq9j4bu3Llz2v8uLS0lMTGRJ598UqdrPTw8OHv2bKUJ5X/8f9YQSOlaz6lUKtavX4+NjQ0vvPACmZmZxMbG0r17d3x9fZkyZQodOnRg4MCBnD9/nrVr1zJ69OiHPiZkamrKb7/9Rnx8PAsXLsTe3p5//etfqFQqOnXqVIt3Jx7G2tqanJwcLl++jJOTE+bm5saOJKphb29PQkICSUlJWFhYcODAAQ4dOkTz5s0B6N27Nx07duS9995j/vz5mJqaEhoaiqurK23atGHSpEnMnj2brl278uyzz/LTTz9x8ODBBjmpzpg2bdqEi4sLHTp0YOPGjdy9e5dx48bpdO3YsWPZvn07ISEhTJw4keTkZLZv327gxLVPehjqOR8fH5YvX050dDRDhw5l2rRpeHt7s2DBAvr06UNYWBj79u1j2LBhhIeH8/rrr7NgwYIa33f16tW0adOGt956C19fX44dO8b69eu1k3yE8Q0aNAgnJydGjBjBt99+a+w44gEWLVqEjY0NY8eO5ZVXXiExMZHQ0FCysrK4ceMGSqWSDRs2YGdnx2uvvcYbb7zBE088wdq1awEYMGAAISEhbNu2jaFDh/Lpp5+yatUqevfubeQ7a1jefvttNm7ciJ+fH1euXGHz5s06r03zxBNPsG3bNi5fvsxLL73Exx9/XGnYqaFQaGT1FyGEEELUQHoYhBBCCFEjKRiEEEIIUSMpGIQQQghRIykYhBBCCFEjKRiEEEIIUSMpGIQQQghRIykYhBBCCFEjKRiEEEIIUSMpGIQQQghRo/8FfrfXYgbPFXUAAAAASUVORK5CYII= name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_dni for model in MODELS.itervalues()
], axis=1)
dni.rename(columns=MODELS, inplace=True)
dni.insert(0, 'cs', ['dni']*len(dni))
# DHI
avg_dhi = pd.concat([
pd.Series(data=station_atm_params_3min_clear['dhi'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
dhi = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_dhi' % model] - station_atm_params_3min_clear['dhi'])[is_bright[station_id]],
name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_dhi for model in MODELS.itervalues()
], axis=1)
dhi.rename(columns=MODELS, inplace=True)
dhi.insert(0, 'cs', ['dhi']*len(dhi))
# GHI
avg_ghi = pd.concat([
pd.Series(data=station_atm_params_3min_clear['ghi'][is_bright[station_id]], name=station_id)
for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0)
ghi = pd.concat([
pd.concat([
pd.Series(
data=(station_atm_params_3min_clear['%s_ghi' % model] - station_atm_params_3min_clear['ghi'])[is_bright[station_id]],
name=station_id
) for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems()
], axis=1).mean(axis=0) / avg_ghi for model in MODELS.itervalues()
], axis=1)
ghi.rename(columns=MODELS, inplace=True)
ghi.insert(0, 'cs', ['ghi']*len(ghi))
# melt
data = pd.concat([dni, dhi, ghi])
data = pd.melt(
data, id_vars=['cs'], value_vars=MODELS.values(),
var_name='model', value_name='mbe'
)
meanlineprops = dict(linestyle='--', linewidth=1.5, color='black')
sns.boxplot(x='model', y='mbe', hue='cs', data=data, whis=[5, 95], showmeans=True, meanline=True, meanprops=meanlineprops)
plt.title('All SURFRAD stations by model')
plt.savefig('all_stations_all_years_by_model_and_cs.png')
Use of Measured Aerosol Optical Depth and Precipitable Water to Model Clear Sky Irradiance by Mark A. Mikofski, Clifford W. Hansen, William F. Holmgren and Gregory M. Kimball is licensed under a Creative Commons Attribution 4.0 International License.