{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"
![](\"Logos/NOAA_logo.png\")
\n",
"
\n",
"
![](\"Logos/NWS_logo.png\")
\n",
"
\n",
" \n",
" Time Series Data Analysis with Python\n",
"
\n",
" Deanna Spindler
\n",
" IMSG at NCEP/EMC
\n",
" Verification, Post-Processing and Product Generation Branch\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"What is pandas?
\n",
"\n",
"The name comes from panel data, a statistics term for multidimensional datasets.
\n",
"A high-perfomance open source library for tabular data manipulation and analysis\n",
" developed by Wes McKinney in 2008.
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"What does it do?
\n",
"\n",
"- Process a variety of data sets in different formats: time series,
\n",
" heterogeneous tables, and matrix data.\n",
" - Provides a suite of data structures: \n",
"
\n",
" - Series (1D => think columns), \n",
"
- DataFrame (2D => think tables or spreadsheets)\n",
"
- Panel (3D => a matrix of tables, I prefer to just use xarray)\n",
"
\n",
" - Missing data can be ignored, converted to 0, etc\n",
"
- Facilitates many operations: subsetting, slicing, filtering, merging,
\n",
" grouping, re-ordering and re-shaping\n",
" - Integrates well with other Python libraries, such as statsmodels and SciPy\n",
"
- Easily load/import data from CSV and DB/SQL\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Example
\n",
" \n",
"Using pandas DataFrames to validate WAVEWATCH III model output with NDBC buoy data.
\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\n",
"- Creating DataFrames from ASCII files and remote data servers (OPeNDAP)
\n",
" - DataFrame manipulations (selecting, merging, grouping)
\n",
" - Basic descriptive statistics (RMS, Bias, Cross-Correlation, Scatter Index)\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"![](\"ww3_prod_val_41008.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "notes"
}
},
"source": [
"\n",
"Python libraries provide a complete took-kit for data science and analysis:
\n",
"\n",
"- NumPy: general purpose array functionality with emphasis on numeric computation\n",
"
- SciPy: numerical computing\n",
"
- Matplotlib: graphics\n",
"
- pandas: 1D and 2D array-like types (series and data frames)\n",
"
- Statstool: statistical analysis\n",
"
- Scikit-Learn: Machine learning\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"First things first
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy\n",
"import pandas\n",
"import tarfile\n",
"import xarray \n",
"import netCDF4\n",
"from datetime import datetime,timedelta\n",
"from dateutil.relativedelta import relativedelta"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Start by choosing a buoy and period of interest
"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"buoyID='41008'\n",
"year=2018\n",
"month=6"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Get the quality controlled NDBC data
"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"url='https://dods.ndbc.noaa.gov/thredds/dodsC/data/stdmet/'+buoyID+'/'+ \\\n",
" buoyID+'h'+str(year)+'.nc'\n",
"ncdata=xarray.open_dataset(url,decode_times=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Select the specific month of the year
"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"startat=datetime(year,month,1)\n",
"stopat=startat+relativedelta(months=1)-relativedelta(days=1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Subset the data to this time period, and make it a DataFrame
"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"data=ncdata.sel(time=slice(startat,stopat)).to_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Take a look at what is there
"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Index(['wind_dir', 'wind_spd', 'gust', 'wave_height', 'dominant_wpd',\n",
" 'average_wpd', 'mean_wave_dir', 'air_pressure', 'air_temperature',\n",
" 'sea_surface_temperature', 'dewpt_temperature', 'visibility',\n",
" 'water_level'],\n",
" dtype='object')"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data.keys()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Another way is to look at the first few rows of the DataFrame
"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" | \n",
" wind_dir | \n",
" wind_spd | \n",
" gust | \n",
" wave_height | \n",
" dominant_wpd | \n",
" average_wpd | \n",
" mean_wave_dir | \n",
" air_pressure | \n",
" air_temperature | \n",
" sea_surface_temperature | \n",
" dewpt_temperature | \n",
" visibility | \n",
" water_level | \n",
"
\n",
" \n",
" | latitude | \n",
" longitude | \n",
" time | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 00:50:00 | \n",
" 172.0 | \n",
" 8.2 | \n",
" 9.4 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
" 00:00:03.950000 | \n",
" 100.0 | \n",
" 1014.500000 | \n",
" 25.700001 | \n",
" 26.100000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 01:50:00 | \n",
" 178.0 | \n",
" 7.5 | \n",
" 8.5 | \n",
" 0.89 | \n",
" 00:00:07.140000 | \n",
" 00:00:03.940000 | \n",
" 109.0 | \n",
" 1015.200012 | \n",
" 25.700001 | \n",
" 25.799999 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 02:50:00 | \n",
" 186.0 | \n",
" 7.0 | \n",
" 8.2 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
" 00:00:03.900000 | \n",
" 113.0 | \n",
" 1015.200012 | \n",
" 25.700001 | \n",
" 25.700001 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 03:50:00 | \n",
" 200.0 | \n",
" 6.9 | \n",
" 8.8 | \n",
" 0.86 | \n",
" 00:00:07.690000 | \n",
" 00:00:04.160000 | \n",
" 103.0 | \n",
" 1015.299988 | \n",
" 25.799999 | \n",
" 25.600000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 04:50:00 | \n",
" 205.0 | \n",
" 6.3 | \n",
" 7.3 | \n",
" 0.79 | \n",
" 00:00:03.700000 | \n",
" 00:00:04.080000 | \n",
" 171.0 | \n",
" 1014.799988 | \n",
" 25.600000 | \n",
" 25.600000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
"
\n",
"
I like to index on a datetime stamp, so let's reset the index
"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"data=data.reset_index()\n",
"data.index=data['time'].dt.round('1H')\n",
"data.index.name='datetime'"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" latitude | \n",
" longitude | \n",
" time | \n",
" wind_dir | \n",
" wind_spd | \n",
" gust | \n",
" wave_height | \n",
" dominant_wpd | \n",
" average_wpd | \n",
" mean_wave_dir | \n",
" air_pressure | \n",
" air_temperature | \n",
" sea_surface_temperature | \n",
" dewpt_temperature | \n",
" visibility | \n",
" water_level | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 01:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 00:50:00 | \n",
" 172.0 | \n",
" 8.2 | \n",
" 9.4 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
" 00:00:03.950000 | \n",
" 100.0 | \n",
" 1014.500000 | \n",
" 25.700001 | \n",
" 26.100000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 02:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 01:50:00 | \n",
" 178.0 | \n",
" 7.5 | \n",
" 8.5 | \n",
" 0.89 | \n",
" 00:00:07.140000 | \n",
" 00:00:03.940000 | \n",
" 109.0 | \n",
" 1015.200012 | \n",
" 25.700001 | \n",
" 25.799999 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 03:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 02:50:00 | \n",
" 186.0 | \n",
" 7.0 | \n",
" 8.2 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
" 00:00:03.900000 | \n",
" 113.0 | \n",
" 1015.200012 | \n",
" 25.700001 | \n",
" 25.700001 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 04:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 03:50:00 | \n",
" 200.0 | \n",
" 6.9 | \n",
" 8.8 | \n",
" 0.86 | \n",
" 00:00:07.690000 | \n",
" 00:00:04.160000 | \n",
" 103.0 | \n",
" 1015.299988 | \n",
" 25.799999 | \n",
" 25.600000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | 2018-06-01 05:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 2018-06-01 04:50:00 | \n",
" 205.0 | \n",
" 6.3 | \n",
" 7.3 | \n",
" 0.79 | \n",
" 00:00:03.700000 | \n",
" 00:00:04.080000 | \n",
" 171.0 | \n",
" 1014.799988 | \n",
" 25.600000 | \n",
" 25.600000 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
"
\n",
"
The parameter names are long, and there are columns that are not used.
\n",
"\n",
"Let's fix that
"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"params={'wind_dir':'udir',\n",
" 'wind_spd':'u10',\n",
" 'wave_height':'Hs',\n",
" 'dominant_wpd':'Tp',\n",
" 'longitude':'lon',\n",
" 'latitude':'lat'}\n",
"dropkeys=[key for key in data if key not in params]\n",
"data.drop(dropkeys,axis=1,inplace=True)\n",
"data.rename(columns=params,inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" lat | \n",
" lon | \n",
" udir | \n",
" u10 | \n",
" Hs | \n",
" Tp | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 01:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 172.0 | \n",
" 8.2 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
"
\n",
" \n",
" | 2018-06-01 02:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 178.0 | \n",
" 7.5 | \n",
" 0.89 | \n",
" 00:00:07.140000 | \n",
"
\n",
" \n",
" | 2018-06-01 03:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 186.0 | \n",
" 7.0 | \n",
" 0.80 | \n",
" 00:00:07.690000 | \n",
"
\n",
" \n",
" | 2018-06-01 04:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 200.0 | \n",
" 6.9 | \n",
" 0.86 | \n",
" 00:00:07.690000 | \n",
"
\n",
" \n",
" | 2018-06-01 05:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 205.0 | \n",
" 6.3 | \n",
" 0.79 | \n",
" 00:00:03.700000 | \n",
"
\n",
" \n",
"
\n",
"
The Tp column does not look right...
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "notes"
}
},
"source": [
"Pandas represents Timedeltas in nanosecond resolution using 64 bit integers"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"data.Tp=data.Tp.astype('timedelta64[s]').astype(float)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" lat | \n",
" lon | \n",
" udir | \n",
" u10 | \n",
" Hs | \n",
" Tp | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 01:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 172.0 | \n",
" 8.2 | \n",
" 0.80 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 02:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 178.0 | \n",
" 7.5 | \n",
" 0.89 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 03:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 186.0 | \n",
" 7.0 | \n",
" 0.80 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 04:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 200.0 | \n",
" 6.9 | \n",
" 0.86 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 05:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 205.0 | \n",
" 6.3 | \n",
" 0.79 | \n",
" 3.0 | \n",
"
\n",
" \n",
"
\n",
"
Suppose we want to look at the data in a specific column,
\n",
"or a Series:
"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"Hs=data['Hs']"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/export/emc-lw-tspindle/wd20ts/anaconda3/lib/python3.6/site-packages/numpy/lib/function_base.py:4291: RuntimeWarning: Invalid value encountered in percentile\n",
" interpolation=interpolation)\n"
]
}
],
"source": [
"quartiles=numpy.percentile(Hs,[25,50,75])\n",
"hs_min,hs_max=Hs.min(),Hs.max()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Min: 0.290\n",
"Q1: nan\n",
"Median: nan\n",
"Q3: nan\n",
"Max: 1.130\n"
]
}
],
"source": [
"print('Min: %.3f' % hs_min)\n",
"print('Q1: %.3f' % quartiles[0])\n",
"print('Median: %.3f' % quartiles[1])\n",
"print('Q3: %.3f' % quartiles[2])\n",
"print('Max: %.3f' % hs_max)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Get rid of the NaN's
"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"Hs=data['Hs'].dropna()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Min: 0.290\n",
"Q1: 0.460\n",
"Median: 0.535\n",
"Q3: 0.630\n",
"Max: 1.130\n"
]
}
],
"source": [
"quartiles=numpy.percentile(Hs,[25,50,75])\n",
"print('Min: %.3f' % Hs.min())\n",
"print('Q1: %.3f' % quartiles[0])\n",
"print('Median: %.3f' % quartiles[1])\n",
"print('Q3: %.3f' % quartiles[2])\n",
"print('Max: %.3f' % Hs.max())"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"It's nice to take a quick look at the data
"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEQCAYAAABLMTQcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXmYHFd5Lv6eWnqbnkXSSCNbki3Lki3LG8byvjVgsLEJDvwCAbJBAOcmQBb4ketcAiGsDrmBQEJIHJOwhj2Ag41tbGu8CW+yrM3aF2udGW2zT29Vdf+o+qpOVZ2qrpmununpqfd5/MjTU1N1uurUd97zfhszDAMJEiRIkKC1IM30ABIkSJAgQfxIjHuCBAkStCAS454gQYIELYjEuCdIkCBBCyIx7gkSJEjQgkiMe4IECRK0IBLjniBBggQtiMS4J0iQIEELIjHuCRIkSNCCUGbqwt3d3cby5ctn6vI+jI2Noa2tbaaH0dRI7lE0JPcpGpL7VBuie7Rhw4YThmEsrPW3M2bcly9fjhdeeGGmLu9Db28vCoXCTA+jqZHco2hI7lM0JPepNkT3iDH2SpS/TWSZBAkSJGhBJMY9QYIECVoQiXFPkCBBghZEYtwTJEiQoAWRGPcECRIkaEEkxj1BggQJWhCJcU/Q0ihVNfxs4xEkHccSzDXMWJx7ggTTgX94eBfueWIfOnMqXnP+opkeToIE04aEuSdoabxycgwAUCxrMzySBAmmF4lxT9DSKFd1AEBKSaZ6grmFZMYnaGmUEuOeYI4imfEJWhpk3BN/aoK5hsS4J2hpkCxT0fQZHkmCBNOLxLgnaGmUqqYjNTHuCeYaEuOeoKVRspl7osskmFtIjHuClkYiyySYq0iMe4KWBrP+rSbMPcEcQ2LcE7Q0JMk07+WEuSeYY0iMe4KWhmwZ92pi3BPMMSTGPUFLQ2amcU8cqgnmGhLjnqClQcy9oifMPcHcQmLcE7Q0bONeTZh7grmFxLgnmBNIQiETzDUkxj1BS0O3isokskyCuYbEuCdoaeiWGpPIMgnmGhLj7kGpquHJ3cdnehgJYoJuWfdqwtwTzDEkxt2DT973Mn7v689hZ9/ITA8lQQywZZlEc08wx5AYdw82HRoE4NQkSTC7YcsySZx7gjmGxLh7MFauAgDSanJrWgEJc08wV5FYMA/GSqZxTwpNtQZIc0+Me+MxXKzASFpeNQ0S4+7BqGXcNT2ZpK2ARJaZHgwMF3HJJx/G1x7fO9NDSWChpnFnjP0HY2yAMbY14PeMMfYVxtgexthmxtir4x/m9KFYMRleEl3RGkhkmenBwEgJAHD/5mMzPJIEhCjM/RsAbg35/RsBrLL+uxPA1+of1syjmjD3lgA9xkRmayxU2TQlySLaPKhp3A3DeALAqZBD7gDwLcPEMwC6GGNnxDXAmUJiDFoDxNyTeu6NhSIn1TebDXFo7ksAHOJ+Pmx9NquRyDLTg9+991n83YM7GnZ+Mu5JPffGQpUaz9y/++wr+M2vPt2w87calBjOwQSfCZdvxtidMKUb9PT0oLe3N4bLx4PR0VHXeDa+tAnakThuT+vAe4/iwFN7xvDUnhO4KtMX63kJ5XIFAHBqcHja5lsj7lOz4/i4adRHx4uRv/tk79PHHhwDANz38Dp0pERmp/VQz1yKw3odBrCM+3kpgKOiAw3DuAfAPQCwdu1ao1AoxHD5eNDb24tCoQA8eD8AYM2FF6OwpmdmB9VksO9RnLDud6PmgrzuIaBaRSbXhkLhxoZcw4uG3Kcmx8GT48AT6yArauTvPtn7tPDpR3B8pIT2sy9E4fxFUxvoLEM9cykOWeY+AL9vRc1cDWDIMIxZ7zJPZJnGYzpiorUkWmZaYKDxvo15ORUAcHK03LBrtBJqMnfG2PcAFAB0M8YOA/gbACoAGIbxrwAeAHAbgD0AxgG8p1GDnU4k0TKNx32bhBu8WOGEQibPs5Ew7HyCxhl3CnLQEuIVCTWNu2EY76zxewPAB2IbUZMgiZZpPP7s+y81/BpOKGRiEBoJelsa+d7QriBZqKMhyVANQMLcWwOGHQqZPM9Ggu5zI98b2hUkC3U0JMadA68BJxOoNaAl9dynBdPBhYixJ8QrGhLjzoGvJ5NMoNaA04kpMe6NRX3vy8tHh7G7P7yHAj3DRJaJhiSQmwNv0JPCYdOHBW2phpyX34lVkufZUNQb+HTbV54EABy4+/bAY6gPbrKrjoaEuXPgPf1J6Fzj0ZaSAQDndLc15Py8Pa9oelKOtoGYjjubyDKTQ2LcOfDbvYS5NxaGYWCiogFwYtHjBj3DlCLBMJJn2kg0et2saHriP5kkEuPOgd/uJeygsShVdZtZ6w261xTjnlbMaZ4808bBaDB3d/nDEs09EhLjzqGSTKBpQ9Fi7UDjmDudNq2Y8k8pcao2DPU8wihyGb8wJw7VaEiMOweeuSdZcI3FeJkz7g261cTc29KWcecWlATxoh7jzs+FILgj2ZJ3MwrmfLTMNZ9/FAzA566R3Q7VZAvfUExwhrZRsgztCNpS5jSPYkQSTA31yDLU2jIMSZjy5DHnjfuxoaL1f22JQ3Uawd/fhsky1lqdT5vTfCJh7g1DPY8winHn2XoSChkNiSzDgdfZk1DIxoIkE0ViDXeokiyTMPfmxNhkmXuTau47+0aw/K77sWdgdKaHAiAx7i6UXZp7c06gVgERMUVmDWPuZNxzxNwT494w1PMIj1vNtc3ziE/kIl5N+m727hwAAHzr1wdmdByExLhzSEIh40NV0/HOe57B+r0nhL8nw6tKUsMWUlo08qlElmk06tHc+4cd4x4UCaPPgrpPZy/IAQC2HR2e4ZGYSIw7B565N+sEmi04OljEr/edxEd/tFn4e3pXFblxsgxdoy1NDtXa2/8EU0M9zL1vuGj/f1Czj+oscKjSPRiaqMzsQCwkxp1DsZIw97jw6ftfBgDMa1OFvydWrchSw2WZvKW5FxPm3jDU8wQHOOMeFK7q1tybk3jRPG4UWZksEuPOgX/5m9VpM1vwq5f7AQDzcuKiYGR4U7LUwDh381+HuSfGvVGop24Pz3SDEs3497FZiRctQHqjazFExKw07t/+9QG8cnIs9vOScU8pjdOB5xrmB1R8NGzmzhr2MhCDyiWhkA1HPU+wzBn0coBx11wZqs3J3GkeN4vpmHXGfaKs4eM/34Z33PNM7OcuWhMrn1aSLLiY0J4Rp1LQC6BIrGELqR0to1qyTMLcG4Z61mdeZw/S3DWXQ7VJrKcHNPRmIYazzriT0T09Hn8HdNL78mmlaSfQbIEqMwDBLyKxalWWGsfcrdPKEoMqs6YNoWsF1CPL8FJMECunciCM1ZZl+oeLeHR7/5THM1XQnG6W0tKzzrg3smgQTbK2tNK0ut5sgGEY9nMK3GbzskyDmTtj5iKSdGNqHOKSZYJYL5GEtCLV3FW/855n8N5vvjDtDFpLZJn6QJ5yBhb7uYsVDYwBWbX2BEoQDH5yB22zidyojYyWsQYiMQZVlpIFu4EQPUJdN7B+z4maTLZc1e2dXhB5I0OdUeWau+rDgxMAgLFpDn2lMTZqPk8Ws864BxmLOFCsaMgoMhRZSmSZOsAvjEHMnU9iatQ6SracjHsj585ch8iAf+vXB/Cue5+1I6eCUNZ05KxEs6AwR1qYM4pc06Gaszp8RSlrECeMJBSyPtDKPlHRYr+JxYqOjCpBlRvn5JsLiFKjx3aoTkP5AVkCUjJLZJkGQvQE958wI9qOWkw6COWqbrdcDNpd0fuYVmvvwKgK6HQb9yQUsk7wxuLfntgX67lLVQ1pRYYsSYnzrQ7wL18QW9b5JKaGa+4MqiI1bQhdK0Bkz6I+1XJVt8NVgx2qDnOvtasm5j5amt7oKBpWs5iOWWfc+W3+8wdOxXpum7lLLGnWUQf4rXWgLEPRMhJz/Rwn6BGSLDMXOvhUNB3j5eq07zzrqS1T1nQ70SzQoepi7tFkmdHi9DJ3msPNIsvMunru/Moet0vVYe4s0dzrgOZi7kGFoMx/VdnkF5phQIr5iRJzl5gZTz8XNPfLP/0rDBeruPbcBfiv9189fRcOeV1qvUnlqm7nItR0qEZg7lmbuU+zLGMkskxd4Lf8LGbrXtUMqEoSWVEvXLJMDYeqYkVJNIJp2sZdYkjNEVlm2GKr6/eenNbrip5e1NezXNXtmvtBrJw+T6u1n+PMa+7TetlAzDrj7naKxWvdq7oBmTHIDcyanAvgmVW5KtY9KbKAmHsj2I43WmYuGPeZwlQfn2EYnmiZGg5VRa5JvIi5T3copJ44VOsDv7WWYmbumm5AlhgUiSWGoA5UuWzC4Prc5r+KNA3MnZkZs3NBc58piDR3+iSoGBjgvM/E3Gs5VKNEy9Ccmu7nncgydYJ/YHHLMppuQJEkMzwvYe5TBr18OVVGKYC50/1ViLk3YC31JjElC3bjwNsz2pXRZ2Gllkm2I+Ye9N65Nffw58gYEYbpfd4035pF0p2Fxp13qMZr3Ym5y9LciKxoFGhr3d2exunxijDBxSn5a72IDZZlUolxbyj4p0ePknZwfJ8EL8i4U5x7UAiyncSkStCN8IgUIn3T/bg1blFrBnI4u4173A5VXYciM0gMODFawrajQ/FeYI6AJvaSrizKVR2D4/7ONE4nJsn1N3HCcMkyEirVmX/hZjMqmo6/+MFL2Hvc3wCaX8Bp4SajHrR7AxxZJlsjQ5UvPwAAlRBWTqRvuuURfujNQCRmnXHnoy+kmK27phuQGLMz6j70vY2xnn+ugF68JV1ZAED/SNF3DL14xNga0bxa46JlkiSm+rH58BB+uvEIPvLDTb7f8WZU98gxkZg7RcsE7JjtOHdFCj0O4Jn79Bp3fjFphrkWybgzxm5ljO1kjO1hjN0l+P1ZjLF1jLGNjLHNjLHb4h+qCZeeFbfmbhhQJGaHk2UUOd4LzBHYzH2eZdy5BsgEeozzrGYejeg76cgyZrJUq8e5T6bU7BO7juO+TUcndX4ymsKrcB86zN007kGt8wC/5h7EyHXbuIeXKQAcszDtVSFdrQBnfpdYM4mJMSYD+CqA1wM4DOB5xth9hmG8zB321wB+aBjG1xhjawA8AGB5A8brWhHj7qVY1UzN/dSYWSs+H9BoIkE46BmdScx9SMDcrReB2vANTsRfn9+JlpkbDtXJ+Il+/z+eAwC8+dIzI/+NzaUEiwgfLUO/pvGEtTcseTR3LeA70MJM2adR3v3pl2U45t4EGe5RmPuVAPYYhrHPMIwygO8DuMNzjAGgw/r/TgCTowSTAC/LhIVYTQWabkCRGU6Mmkwzn06M+1TAa+6A2TzBC3rxGsncDd64K60fChmmbccJ0V3kbRk9W5oHYfHm9A5najhUS1Y57ly6NnOnX82kLNMEtj1S+YElAA5xPx8GcJXnmE8CeJgx9iEAbQBuFp2IMXYngDsBoKenB729vZMcLrD3gGME+o+fnNI5RBgdHcXIqIQT0gTevlLG17cCQ6fjO38rYHR0NNL92HzcfJl3bN2EdhV4ccc+9MpHXMfsOGQ+x/3bNwMAnntpG/KndsU63k0D5jg2vrgBA8eqGC9Wp+V5Rr1PcWO47DZmUcYwmXHuPm0uHsPDI76/2zLgGPAnnnwSWYXh5GnTd3V04JTwOqOjo9j+wosAgO1bt0BmwL79B9Db6+eGu/aVoTJg725zjjz59Hp0Z8Xc9FifSc4OvHIQvb19kb9fvTh02JEfn1ofPL7JoJ65FMW4i5Rt75L4TgDfMAzjHxhj1wD4NmPsIsMwXOuXYRj3ALgHANauXWsUCoVJD7gA4PMAfvfeZzFR0VAoXDvpc7jw4P0AgHw+j1TGwJmLu/Dxd1yGjf/yNHJpBYWCdx2bu+jt7UWUZ1Z5uR/Y8AKuXLsWSw5sgpzPoVBY6zrm8DOvANu24uYbr8Enf70OPcvOQaGwMtbxlrb1AS9uwNq1a9G/5RgeO7Qv0vjDsH7PCVy4pBOdWTXwmKj3KW4cHZwAHnvM/jl0DNa8n8w4U3tOAM8+i3x7HoXCDa7fla17DQDXXnc9OrMq/mn7euD0aciZHAqFm3zn6+3txZrzLgSefw5XXn4ZUi89hzOXLkOhcIHv2MeGtiI3cBQXrbkA2LoJa6+4Csu724Tj/FnfRuDoUSxZuhSFwprI369ePHRqC3DoIADgiiuvwtkLxOObDOqZS1GWlsMAlnE/L4VfdnkvgB8CgGEYvwaQAdA9pRFFRFqRYt+GmnHu5i1pSysYmeaqcq0CSh5RZIaF7WmcGvM7VEkyyagyUoqE4WLjZBlZYmbyi27U5ac5dGoc77r3WXzyvm1xDTFWBNXxie381r2rWd7XI4uEVWekMacUM3kwyC9CjXRkK/s0miwTeEhDwMfeN4MEGMW4Pw9gFWPsHMZYCsA7ANznOeYggNcBAGPsApjG/XicA/Uik5Jd4XMDw0UcH/EbkcnAzFA1J08+rUx74aFWAU1sRQoutUsvPiUYNSK6gE9iyqjmVK/HT7PliJn30Aj/QByg73b2gpzdtq4WJhNhE7Z48Kch7Zn+PTpUDHyXXMY9pKaT00jHCoUMEbXpDNOdocon4s2KJCbDMKoAPgjgIQDbYUbFbGOMfYox9mbrsI8AeD9jbBOA7wF4t9HgFuA51W3cr/zco7jis4/Udc6qbkC2Xoq2tBLq5U8QDJrYYXV6+DBFs8Ry/C+is4A4yS9hqfC1sHfATN45d2H92+1GgAxlPq2gohmRDPdkFjtapMWn9Scx8Qbuz3/wkvCcZc18HinZNNxBC0ixYpbjJvIVRgbs609ztAzP3JuhB3OkcBDDMB6AGd7If/YJ7v9fBnBdvEMLRzYl48RY2dyuqfXHoxuGAU3XIVvBvClFij0aZ65A5+SQoOQhvhyvKrOG1OPgOzFR8kuxjmc6YrFPFndqdEwgmbLdCuEta7odFx6E4yMlLO7M2Iw4DGSIRU/KzdzNf3nj/tx+cWMdMuZpVUY2JQcuvqWqm7mHhbUa9uISeEhDMOuYe7Mim5JRruq47ctPxnI+A05tGcBkEkHlahOEg5dcVElsuA2buTMoUmNkGbqGLLFYmDv5YBqtbU8VNK72jOnsjUJObvjCOnzguy9O6vyiHYG7toxblgGA83ryoedMyRKyqhy4Wy5WNKRV2a7/H6q56/TvDCYxJcZ96sipJjvZZzXhrRdU7Ie2fWZzh5l/QLMRtuGWGJQAPd3VJUlmDUn64K9Bmns9xp104+mKJ58syJh359MAgOGIvoGHX+6PdFxQVy1AzNyruoE3rOkBACzuzAr/rsRp7tmUjImA51Os6siosr0TKYWUNJgxWSZh7vEgm3KGHkf0ggG35p6SpZZPV28U3HXUxffRrvtiZY82grnzO4i0xdynKrX9cssxO12/WeU6WnQWtZvGParjN6rK5DB3/+8Mgeau6wYyqow1Z3RgIiCRieZGWpGQ8wRJ8ChVNGQUKdIiTXa1Ecx9/d4TWH7X/djZN+L7XbOVH5i1xp0vGvaN9QfqPp9uMXfS3FVZgqYbTbECzza4DbfYWUoGgpFDtQHMnb8G1QmaKnP/2UtOElbzGndzXIs6LOMuqMYpApWAqAXSuYWNOQTRMpphypxhjHyysowtr4XunhrH3B/eZu5yntztDwbkp3kz2I1Za9yDXrCpBukQc+dlGaA5qrvNNvAhiEF6OrEqmVFETeMcqrLEkKZQyIqOUlXDxZ98CPdvPhb5XF1ZxwCGSQKNQu/OAaz++C8xEpIPQO9ET3sGADAYkbl3hSRk8SBDLLJbbs3d/FfXzXufSwUb7XJVhyIxSBJDNqUELgKlqm4yd3uRDpNlzH/jNLBVTccbvvS4TSRFJIGXZZohWqbljPtUpRSaB5TERMa9WVlaM8PpgGS2txM9E29/00aEQrri3DnmfmqsjJFiFX/7P9HlvM6cYwBnQq770q92oVjRsWfAX0udUPYy94jGfXFnJtJxdH6RYRPVc6edcEYNlluI3QNAVpUCj6toBlRFshfpcFnG79CtF2MlDbv6nXsvWoQ03ZixcsMizFrjLgcIhT/ecHhK5yMbrtiau/lvs0ZGNDNcoZCyuOelE6Zo3vNGRBdo3DVsrbaq2XNnMgs3H24bVsK2YbDGHHaXbFnGYu6nx4MrbZLTFYje3JoWtVqZ23Q+zTAgWcw9iJHrVg8FwCz7Ox6gzVd1k+HTIh327GyHboy7Qe+CPlEWh/em7CSrxLhPGe+74Rzh5x/76dYpnY9PvAESWaYe0LxmjNn9aL1ymW6YLIcxBrVhoZCO9t9l6cqnxyp25cHJRL3wMeQzwdypGXyYk5CISFdOhSKxUCOsG4ZdPjfq9xm0FovRUtX3PIUZqroBWUKoo1Q3nHcuTJvXrHLcUZi7KBSzXnjtgGic42XNDkNNmHsdaEsruPKc+bGdb8sJ82HxDlUgYe5TgVuWoUXSb9zpXjfKocrr+vNyKlKyhP7hoi0BTYa5lyo62jMKrjpnwYxo7sRuw+YjLUApWUJ7RgnV56uajrevXYbXrl4UeZHrs5quaLqBv/zxZtffuaNlzH81o7YsQ4s8AGRVGRXNEBIq8oelFQmMhe+eRBmy9cJLPkTRP4PjZXTnU/Z4Zxqz1rgDYufpq8/qmtK5vr7VZCVe5p6EQ04eblnGvJ/eF1Y3HINlFoxqhEPV/FdiDIwxLOpIo3+4WCONXoxS1cyETqszEyJLzD2sJEa5qkOVTedke0bF8EQYc7cczUpwyr8XA1xd/h9tOOxySPP30s4Q1UxZZnFHBiOlqt0Eh4fBzYOw+j9U1I+yjcMyjR1ZKNLXigTvMz8qaEAzNFHBAsu4T3ddGxFmtXGnl/fsBTn7MzLOU4XCxbkDCXOfCvhQSEUS97zkGZupyzcuiYlZs3xxRwb9w6UpMToqXJWW469GGgVU8iBItgBMo0hJPjWZu6VhT6bMRt9wEV2cY3lnvxPrLSw/YDH3i5d2AgA2Hx70ndPsW2z+f9g7V9V1myikleAyBeb1HVkoLnjJybYjQ655ZBgGhiYqmN9m+jKSOPc6QQ9xQZsTphY2+aOAFgc1Ye5TBl9awGbuHuPNMzZFYg2qCuksMgDQlUvh9Hh5Sn4UKlyVmgTTjRNEWcIaiZerul1DxzTuwcydSm1EZe6abmBwvILlXI3y41xvXFcXIj5aRmJYucgsPXDw1LjvvDoXLaMG+Ll03XBp8xlVqpGh6ow5LnjHNFbWXDuR8bKGimbYtiiRZeoEzaf5nHGvt5Ij6cDpSTB3XTdw36ajTeFEaQbwmrsSUOiJZ2xBETV1j4PCW20/iuncncq1zAJ1UmAJYy/6h4tYv+fEpK8TBFqgHtoW3FmoVNVsObEjo4Ya96pleKMyd2qVR60TAbdB5+8IH4ooSQxtVvNr0cKkG86uJIi5005QsY27HJrEZMtCDXSoAm6nLoWdJpp7THjPdcsBAH9007kAgBXdbaHMJgpqsQgRfrThEP70exvxnWdeqevarQJvaQFALMtIkqO5NybO3QmFBMxnq02xYYeZRCNHbrR9+1eewrvufXbS1wmCpW7h0R0DgXJDycXc1UBZRtcNGLbmLkciMNRwg4+JdyUSuTR381+q1ZS1wkhFxEvnZZmA3bITyWb+PhvioAUaU36gXPWfi5fnRq26QxSVpTXBjn9WG/c7XrUEB+6+HVcsn48Dd9+O61d1C2WZFw+ejvygyRgRiwhLGiEcHTSdKycFDqO5CNuRGeJQ5WUZuVEZqlxtGcB8thVdn9K1qLS02Wi79otLTdYbgSBGPl7WkLVYcpgswzPhVMSOZlQ07QzOuPPvmre2jGEYttNcsuSfoKxOuQZzJxZMzL0trYQ23Z4+5u58RmOmhSxh7jEjK0hzfnL3cbz1X9bjPyPWnyGHardVfClKUhQ9eGJNcx08G7MdqrqAuZMsIzXKoWr+K/HMXZtavSByqKYiyjKEuHrW8Brz4ISYREyUNWRVkmUUjJarQlLDM+G0Vf20FvkZtY27I8vw7NnrUHUyvilBSVyCwCXLKAGyjObOQWlLK6Gt+xpSfkAwP/lFkWxA1sodaAaJtqWsUVY1t5j8jT10yuzAvktQxU0EMkZLurK4dFm0sEp6sFFbm81mHB8p4dCp8VBj4HKSBWT6msadl2Ua51Dlx1LVjSmVFy5WTYfqZAvKxbUj4TXmoLICExUNOZu5qzAMYFTAcKu2cXdLIceGJgKvP1Yyr0+aMl2P4NXcyRg6xl1cN8aU58z/D2rEUeV68gJAPi3bi40IBufQjQsiWYZn7nRPKZwzYe4xg7ZE/CQivTVqWVPa+gHA0nnZSNE39AJH6WYz23HFZx/BDV9Yh689vjfwGJ6NBUUdaTp3TEQde7LgM2UBSpYyprSQlCo60hE7AfGI63sNT1RtSWQwoNqjKcs4oZCAWMJxmLBkh05+55lXcM3nH8PGg6eF5x4tmdfMZ5zmbXypAHecu9Mwg49hFztUnUU+kLl7ssfbUoq92IgwHdEygNuhSr+n8ggJc48ZdsmAiKFqoi2zwrFvb5/WIJDh4heGVgQfHfHU7uBIEF5y6bDSsQc9dU6qmm7X71EkhuEa9UqmAl4eMq9jFiibmkPV0twDfAhBiGNHMjRRwcFT47jKysgOZO7lqk1wKA1e5FT1au4A8IT1PHf3i31Mo5YxbU87ce4uWYbj7oZh2NcgvhNUN0Y3BFnhPubu1tzzmfDm9Y0oHEbP++cfuA7fv/NqAG7mTgQvnTD3xmCyWaWiZ8+zb5GGLwItJmGdaloB/JpJ208RdK4ufqdVTvYPv/GCyyhVNN1m9fRC/JyrmR4HeFYIWPH0umHXlpkMihUnWgaILrfE0WFqx7FhAMB1K7sBBJfyNWWZ2sydl0zIT1SrGBo9O565u2QZj+bON0oBrAiXAFmGHlG6JnM3f59Pm/6EIH9GI5n7wva0vYPimXvVlmYlK28jiZaJFSJvOz1/cfcYP3j2HVbIiActJq2ezeo27sGNlzXOqPIZjcMu427YhvIPrzOLwD2z72ScwzWjNbjnKVua+1RSw/k4d2B6mTvpy0vmmc7MoK5G42WnWXxHNoS5c0w4HZEQbT82jAVtKczLqXj0Izfh9kvOcBHsdTbVAAAgAElEQVQfr+aue6SUbEDxML4qZBA5E0XLGEZwTotTOCz0K00KZU56zQi6etFir0hSw2TGyaK1jPukmbtIluGYu8BBKwKt4C1v3Lnb8MutfVh+1/3CkD+DM6qdXCMI/naXNd02lAvyadywqhubDg3FOl7DcMsyqmQ6Q3nWHSVEtqrpqOqGS5YJe9YbBxzjG8dLTucgiSvIqE2U/cxdVF+GFhyJN+70fQKUxa1HhnDJ0k4wxnDuwjxWLsyjxL8bngxVR5ZxmLs4WoZ3vgclvLmds21p87sFSTM0lDgjsGh3npIlYVcvupYqm+G/zdB/ubWMu4C5hzlURe81H/FCL0ot9r7bioUva83ZODkuiBa5LYf9BtmbfUrgJYoKp7kDwCVLO7Gzf6SuBtZeeGUZSmLijUeUWGgqUpVWpEiloNcddIxOHNorMUSzQbQknI8VawEizZ3+Fd1PPnSXHKr2OxMw3MOnJ7C82yk9QOzV7qvKHWs6VN3GPajgWpRQSIcVW87ZGo10SGuPMwLLjohTuLLDglBIVZaQUuSmaPLTUsZdVKY3jHWLekFSKCQALrMu2HkzWqpi3/Ex33VbEaKvJ7q/XqNK4A1ihWPuAHDJ0i5ouoELPvFgqLNsMtB0+DR3wO0Ii6LLkoE0mXttzX1Zu/O94tBebWddSBNpYsUULRO2iy1xLNTbcUy0GI2WqhgtVdHT4SQweQ0xv0ZqOsfcuQQl0fuhu4iAeFcU1GuhlnGP08A60pBk73ZEDlVFZlbuwMzbgpYy7iJWFWZwRaSNj5ahbL+ioOsK4TSXlTonjbvgJnq17jtedSYAoMLpOpWq4TLuly41cwoMA9jRNxzLeHWPLEOSm3s7Xdu4O8w5mubOtySNY3tuN5FWgptI0zFpi5A4jtJg455WJZ8sI2L6/Vap354Op3uTfX6LvfISp8aFm0qcQRYad06WScvUPMR9z2zD6qnYGvQMaChxdswqVTQwZi5AjDG0exKpqpwmr8qsKWxBSxp3/sYSc4kaFeVyqBJzrwQzyZLgWq0KQR6H8AXzhiC+9dVLAbjvT1nTXQspX7MkrnwBg6tfAzjPln/ptQjGl2fuNOawZ83/Kg4GV7YMqCpLgU5+ug5JXWHM3V4sZNl3nIjt9g+Rcfcz95JAlilVdSeBjJNcROeOlKFKSUye/sZBBpSuHef7WNJ0pGTJHmtHVnVlCvOyjCwx3LfpKPafGIvt+lNBSxl3UZxsWMx7rVBIivQ4PRZcF5t/eZtBZ2skRBKGSELh64UAEMaGV3XdZmCE16/p8R1XD/iGIIDD/PhGD1GcbmTc04rssMaQZ82vF3E49Yj9pxQpMOqENy6Aw25F7JXYtsnc3dq8iLlTGGRX1slOTXuNO/edK5oulFKCmLtXlvFlqHo091p+DxpKlLIKUVGq6K7yIl05FUPj7tBewJxjey2Z9sM/fCmWa08VLWXcRXGyZOhFCQ2iz3g2SdvQ/mF/1xX7/Py1Wty4i5j7qCBTUDMMm+EA4igIrywDOCGRcS2SmleWsTV3jrlH0twdGcOuchnyd27mHoMsw5W3yKlKJOPOmFUULJS5O5o7ST1ffnS3L4qMHId8boP3XeP/oswxd8mWXEyHqqiXrlOGQoLEamvutVpg8u91XOy9rOlIKU74b2dWdeVt0HxQOZ/dxoODof66RqOljLuIudMEEK3y4jh355bQNjTMuEfV91sBovdEVMDJDIV0fhYady6JiVBruz1ZGN4kJmscvHGMpLmTLKNwoZBhsgzP3OPU3GUJmZSMcQG7ptonfLRXWhE3tRBp7jy87fDoHHxuAzF+keZOkTuAs6CSL8B733Qdnl2eP6rGp7nTriHgGfCbpbj63YqYO59MVuWYO49/emxPLNefClrKuAsdqpRgJDLuAubOvxz5tIJcSkZfVObe6pq74OuJSq8GyTJ88SUzzt39IqQiRKJMBronWoaYn0uWiXAtkUP18Z3H8c+P7cbRQX+xLd7JHEeGapnTe3OqjE2H/IzQy9wB0wCLdkH8YiEy7t77z/scnHMHR8uUNd12nntlIlEBOT6wSiTfeDNUw6SxZ/edxBHumZRiCk8uVTU7BBIAOrMpF3Mvc9IR3+rz9AyWAW9J4y6SSkRsUPRa8w+GMbO578BwcG1uYg+KxGJjCc2KqvUGn9eTtz/z1owBKM7duY8pW8pwM3ev5h43c/dFywhkmclo7hlVtqWJb6w/gP/78C6895sv+I5vBHOnhZD6Bf9yi7sjk9i4i3VuR3OX0ZZWfH2HvfXd+Th/gt+hyi3cVd0lJfHHh1UHtcdcI0M1zFn82/c8A8Ap8xzXO1muuuert0dt1SIrjDEXaay37Wc9aCnj7myZ3VtE72cEQ/DcvTrwoo50uCxjTdZ8RpkzzP2Lb3+V/dnASHiGKhAky/g19ygJQl70DRWx01POec/ACI4MTvi0f1EoJK+57xkYxeHT/j6fjuYs2yV1CaK5wQ8/njh33b43f3nrauRSsq/ZtKgyaTqgEUeZM9YZVcaqRXnh7wli5u5JYvIwd/reZBCDDDLfG5XG72Xk3gzVWqGQYdebKkpV3ZaWaAz8wl3VDVvS5ZW+ejvD1YOWMu52nCw3OUo2c/ffZHESk5vFLO7IhMsy1uRpSylzRnNPKRLyVgp4v2BXo3lCIe2Wha44d91n3KOk9ntx9ecfxS3/+ITrs5u/+ASuu/sxs80bJ/04oZB8tIzB/d3juP7v1vmuYTtUFclVOAsAFubTvuN5HjGVImVelLl7JUsM5/W021nR9nXImCrcjimQuev27wHgwjM7AZhFsfjf88ersltuCAuFLFedblf07Ccly3gMsjdDVY2ww6NnFpvmXtXsvsqA0/eXonH43RWPpmfujLFbGWM7GWN7GGN3BRzzdsbYy4yxbYyx/4p3mNGgKn7jQKF6oqgIYhu8IfJuUXssWSaoAh29VO2Z1jfuZJtVWcLmv3kD3nXVWULm6t1qi5yQZU23nxehlqNssjA7Ezlsi4zDRI1omR19wzh0ymHwPHNt8zB3ETPmNfc4mHu56jB3wJxrQxMV1xiDZBmR5l7yyCxndpmBAwusRvM+Waai2fVU+HPzx3pDIX2hmUGyjO6dK/4FyanbEr5QiBClhWAtHB8pYWCk5NLciTSQT6Wq+8kKIA4tnS7UNO6MMRnAVwG8EcAaAO9kjK3xHLMKwF8BuM4wjAsB/HkDxloTtFXkX14qeSpy0tEntJ0izYzHwvY0ypouLMAEOBMsn259WYYMoSqbfTF72jM4NVb2bY+98eWibXRVN1xhY4Cz84paj78WJipO8wqAi3OvkaF66z8+idd/6XH7Z96hyi/+i9rT6Bsu+hZ+TXfCBmNJYtLcxj2XkrHt6DBu+MI6u2KkyLjXZO7WsVSQjJ6Zd0EoVtySBJ0b4EMhrdBHFq65e8/tlWVSgoqKtnNW8SwU3qbrgmcZB+F67zefx77jYy7N3ev8r2qGL1IGaH7mfiWAPYZh7DMMowzg+wDu8BzzfgBfNQzjNAAYhjEQ7zCjQZYYFrSlcJzTgSmaQ+Q487ZhUyT/7SCNNegh0eRpS88B5k6yjDWx57WZRsHbPIJvnQY4Wje9tNSmzifLKH6GHwaedYt2Vl7mTtEWfMGnoPK/fN2QIhcKyeP8xe0oVnRfoxHNcLKb4yo/wBsW/judGDF3lWWh5i4LWfh4qYqU4mRb0kJEI/XJMhXNF1VDP1NTeLr9tKBUImvubglPlMnqWygCmHtRwNLrzZkYL1ex9ciQPTaCnXDFlW0QlcFuds19CYBD3M+Hrc94nAfgPMbY04yxZxhjt8Y1wMliUUcGA5xUQO24RFELNCFpxRWtvDTxg7ZX9FLlM0rLZ6hW7ftl3hMq5+tt+xacoWpY/1ovq1eWiZD9yeNv/2eb/f+ie883r+DHMcHVCooSzVKsmMZV8kh2qxe3A/A7VTXdMcBxO1QBp+YRABT+by9+/tJR+57xum9G9fcaXf3xB3HvU/tdxvrchaZDdc0ZHQDERtPbnIVIz5cf2e36PK3IKGuGbyeRDjDImu52eguZe8BC4T1OZEjrlUW2Hhm2HaR88IDqGQM1cwGAi5Z0cNefOZug1D5EWOHZ+0YoAFYBKABYCuBJxthFhmG4XPqMsTsB3AkAPT096O3tnex4a0KtFrHn6Kh97lMjpi45NDLqu95gybrxmvkCpKD5jtnTZ/7uyV8/iwPt/rVwxz6TuYycOo7xYrUh36lZMFEsAWB4Zv3TyKkMrxw3703v08/i8DyOTZ4oYrxq2PeCdki79uxDLzuMCWuVOLh/P3px2HUNBmDXvv3olWt3ZfrWr53aHY+sewL5FHNt1U8OjqBNH7PHsW/QfNEHx52XdMPGlzBxUBYy/8fWrYPEGPYeKEFmuu/Zto8fBQA8/ORzONrtfP9ytQrdeql37NqD3sorNb9LGPpPTKCswb7+yX63E/u+9Vtxdoc5N1947lnsz1oGtVjG3oEKHn50HVKyO0TPO9f/+uoMsspJ/ATAxs1bkT6+w/7dwWNFVMuG7/uf3SHh5IR5X/Za7wHTqzh4+Chy42ao5gvPP4sDWQm7T5v3/oUXN6J40DQ7o6OjGB2VcMoYt889NjKBEvddAWDHfpM8PLP+aWQsQiAxYPde9zw5Pu43pM+/tAUp7rtMFg/ud4jLtsOnnbl02Pz8iafXozsr4Wh/EeWKeY/+ZLWBP7GGNTg2UZdNGB31262oiGLcDwNYxv28FMBRwTHPGIZRAbCfMbYTprF/nj/IMIx7ANwDAGvXrjUKhcKUBh2GB09uxiPbB0Dnrjz2IAANqUwW3usNDBeBdY8ik0ljpFLConntKBRucB2jbe8HXnoBF7/q1XjVsi7f9TZVdwO7dmHV8mV4rv8V3zVaCQ/sfxhABYWbbkAupaDr0CC+uOFpfObZIr71h1fixvMWAgDu3fMsUuUqCoXr7L9VfvUAzlx2FgqF1WYG5CO/wgXnr0Lh2uWua6Qe+SXOXLIMhcIFtQf04P32/776yqvxs5eO4Au/2ml/xtQ0li/tRqFwKQBg6cAI8MwTrmSsiy6+BDeet9BklA/90nX6K6+9Afm0godObUbbSWdO0XXveO01+OKGXiw+53wULl9q/93fP/9LzO/Mo298GGefswKFwrm1v0sIvrpjPTolCYWC2bvzhdJOPPyKk/l41UUrzR3Ktm248frr7KiXiQXHcP/+F9Fz3mW4dFmXed8f+hUAoMcz1wuA6aB9ah1Wnrfa/j4DI0VsffBRXLK0E4XC9a5xvWZ4G3668QgKhQK2aLuB3buQz2Uwf+E8rFgxH9i6FTdedy0WdWQw//Ag8OzTWL3mYhSsGkK9vb3I5hh6FuVRKFwOAPjm/udwfLTkGts2Yw+wcydeU7jR9qulH30Qi5csRaHguP9u5aKmqK3f4rNXonC9WdbigS3H8CfffREAsO9zt/l2YiL8vP8lLGw/geMjJfR05ew5MLjxCLD1JVx+xVU4p7sN/7rr18joQKFwjfmHj5pzZKIK3HjjTZGuJUJvb++UbUoUWeZ5AKsYY+cwxlIA3gHgPs8xPwPwGgBgjHXDlGn2TWlEdaItrdhtyDTdsGtmCGUZ61/VuvHtGf9alwlpegCY2zKJOV2bWhkkT5OPoourbfszrv8pX8aV0J5xSqSKnH8EUSiceCzu5zle1tC787jvM16fpg4+PEi3F+m1tM0vVXSXLPHLP7sBv/jQ9XYhLa/Poao7ddVjiZbRDFepBt5JDFB0ilVcjLuni6zaSDQ+Xj7i2x8SvBEwANA/ZO4SbrIWbh58r1B6GmnFjFN3ZKLamrtLllEkV8gs4Eg5vAN+flsKJ0fdCXQ7rHyHd1+7HI9+5CYwBgxxSXbf+vUB+/9PC5LvRBgpVrAwn8b33n81/uv9V9ufe3M3TKezfz7rBjA6Q/Vlahp3wzCqAD4I4CEA2wH80DCMbYyxTzHG3mwd9hCAk4yxlwGsA/BRwzDibYgZEeSQ+cB3X3RpsmEOVdKQ2wUvP73UQXr6P6/bY19XN+J5mcPwwxcO4ZYvPVH7wAaA1kdyPPMt9NKcs9Gro9KxVIujLNCHnfOIw/e8oPoqVy6fb/5crmLbEXdXqImy5tKn89zz9UaziBZvMu6m5ux8vwvO6MBFSzrRljY/89bX0Qzu/DHFuQc5VAFzERP5MTJ2sxnze3zhQUeeEJUd8HVlgtNdbK11n3moimQvKo5DVTbLD0SOc/dEVll/z6Oq61Ak5mK/PR1p/PfGI/jDb7jEAQBmV68zu7LoyLiLe/EBE2G5KzzMZDuGa85dgCVdWee7eypYlqp6YF/hmSpBEEWWgWEYDwB4wPPZJ7j/NwB82PpvRpFWzOSC+7ccc30e6lC1Jo03QcU8XzBzp4l6Xk+7q9GwImCkceEvf7wZgBUfPMWt3lRBt5Au2+Ey7s53Ngx/vkBnzqnF4STc+O9TVy4lLGngBeUvkARx5PQExjwOtbKmu5k7Z+hzKQXFStk2QqJkF4qQGilWXQsDQZHNzkjeJtSaQY2SWWwOVf7+5jzMnXcc87shJ9KrCsMwsIVb/PoEyWeicEVRQTKCKjE7zptCIXMpGcWKNqnyA/ypRY0uRNnMVNTvsR1OYN6i9jQGRkp40yVmcxhvcS9+Tg4Ml3Dhmb6v5EMl4H32duQqBUTLAMCu/lGcvaBN+LtGoqUyVAE3gyS0p5XQqpDyFGUZYnZvW7ss9rootTATkTm6YS6ExMr5l4Wf2JonFBIwJRzaItvlUQUvzcJ8Gg9s6cOLB0+HjmXUY9z3WY0RvEaYN4SSxOyf6d8XXjmFt/zL0/jIjzb5rkHFuQbHK0IZg67njUjRdNNIKZIUSw9Vb/ajV5aZKGuuwlX2cdYzmSjrODZUxInRMv7Y0v9FpWht417hmbs74YmHIkswDHOnRkQpl5LdOwlrIthVJIWhkJ5Klt5QSEH2J984hB/r7119tv09urJe5u6cQ8TcP/OLl7FupzuKu2oxdy9o7n75kV0AKNHLf49kiflKRUwXWs64i9hgW1oRvmQUPVA4fyHW9sj4rcuX+Y4Ja1dGHZpyKXnajftMZL5pAkZO4DVpUQ9Vvv413SNvqQfAMdp3/WRz6FhIClllFTF7Zp+pAi6dl3Ud1+kxyqS70yLw+K7j2HhwEM/tP+W7BjH3oYmKy7/AI59R7EQ5AvkclJjarXkzVM/scn/H8bJmxVlLLjmMFoHxchUHTpqL3w0ru/Gh167EV9/1at91ZMmsAc/ndPBdm7zgdWd6u7KqydyrmuGSUoIzVOEac0dGxfBExVdC2Ptez29LuX42DAMjxSo6ss7i3pFVXWG6fKizt6wxANz71H685z/dMk8lIPOUDP46y89TrLo192+85wrc/daLsWpRHpsETeSnAy1n3EUMI58JMu7mv+f1tOODl2WE0TA2c/c43IYmKvjYT7cCsIy77N/SNhIzwtwFiUcEfsfk1VEB9xbZ0Yf956IFYOm8XOhYSJZZ0Z3HovY0ntx9wvo7t+Fb7GF45Fch5u41zADwB9ecDcDZmQ2Ol9GVS/mOA4D2jIoRD3OvGiarNeuPxCPL8Pf9wjM7XL+fqGhC6Yi+Y7GiOX1QOzP4yBvOx0VLOoXXMnci/g5D3pwEwKM7G44sQ8zdWwoBCJBluGnQ05FBWdNx2tPlyDvveOe4YRiYqGjQdAPtGWcR7sqlfM5ugs8JHiCfVTRdSEL4uavphq9EQ+H8RXjHlWfhkqWd2Hx4MLB8SSPRcsZdxNzzaUX48Oh+e51/PGyHqoe5f+lXu2y9ry2l2JMqaDLFjWZj7vzHVcELQcxd1w1hZAfhi283wxaXeYy0F8Tw82kFyy09M61IPlbr3b4Tkyfj4NXLAWB+mxVKaBmM4WLV5Tzm0Z5W/Jq7burRZjRJTA5Vl+au4C2XOXmEE2UNY6WqLxpIlSUoEsN4WbMLvInkDB5tadlO/KNrA+JnRc+4qpmKu8Rg93j19sh1yI973np3eaIGOSLNPZ92DGmpqtuLNC+tdmYV1/vIfy+vX2dM0FGMvpswqov77MRoKTBDdeWiPAbHK74s5ulAyxl3EXNvzyjQDX/4HDmBwtySQZo7P0lzKdluyTcwEs0LXy9EoXuNBmnuIvA7o1LVHxbWmVVhGMBIqRoaCrl2+Xx059PCEs08qKxEPqPY17rwzA5fJImXudPPDqv1L/rz8yZLHy9rtuEONO5ciCdB45h7XG32vKTls2+5yP7/iYpl3FN+n1HWYtL9w0Xk04rQMcyjLeX2IZRDnhVfVoJCGrOqguMjJbx8dNhlACVrsRNFy/DkanGn+R7xmrh3oQDczH2Ce04u5p41nfP03o+WqrhhVTdWL273kbCgcEXRrgFw34+vP7UfugFfFi8Ae8c3PE2kj8ecMO406b1dcRzmHnw+6mbuNab8LiuXVjjGEdzYI07MRFqziLn/6etWAXBHI5Wrus+xTZN8aLzii6TwIqOK65DzoPC+XEq2n/mqRe1421onmWhJV9alwQIOM/TWZSekFAnXnrsAABkNPyPkIXSoGpZDVWZ1yzJm3Rh/Y5OMImOR5Uw+NjiB0VJVGO1FGvj+E2M+yUqE9ozianrO1373wi4XYTlUGYBsyvzs2f2nfP4OUSEz3VMemnZNfPhgper//vwiVaxqNjPmn1NXTnXFmQ9NVNCZVX1aPCBu9A5YddoF85T/7J4nzJSeeW1+6Y5IwXTt6Hm0oHH3b41olfdukeknrz7sRUaRfMaUN+5tKRkL29NgzGweMVVouoGHtvVF0uemU5Y5MjiBTYcGhcz9w68/D4y58whKVc33MvKT3Jvg4kVGlWvW4Xa3ijOfeTYlY+Widhy4+3YcuPt2PH3Xa32SGxn3sie8kLDrM2/EGZ3mMRMVTdiogkd7RvU7VHXToap6GjpMBjv7RrBnYBRVy3B676ckMTz3sZvxsdsuwNGhIl45OS5k5aSBbz48hIsDdHYebWkF6/eetKNp7PssjJYhWcZ0qDLmblTivZ4oQc2b8ObUcuLq/wh8PXkfczfH25FxO1QB4N4n9uHUWBl9Q0X0dGTQlVXxwiunsbvfafIi8r0A4r4DgDhwQnR/uwLqL00HWs64iyYhreZepyolMdWw7chY7CcIubQCVZbQlVVxcmzqzP0/n96PP/r2Bvxi87Gax06nQ/WmL6zDHV992mTuwnhnt/xQFsgyFEo4OOHEloueFWCyxFqLF++UpfMEGWAeSyz2enRwwjbuHRkFqsxw9QozUYccY2YUih567nzGZO685Fe14twViU255O8t//gEbv7i46HOZwC4fPk8AKaMIcrAzWcUbD82jFNjZVwiCBjwYscx0+B96VdmiF+YhMbLMiZzZxjm/A+rF7sdvylBrXav8z0jyCupCPrt8jsv06Hsl2WoINpXHtuDd97zjFmOoCODrpwKTTfwei4ZMIi5V3RxKKRX7gPMKqFe2DvWhLnXD2G0jDXpdxwbdlWOi+rAFsXe8l2cctaLb26Bp250KVzt5GjtBWI6mTstirphCMsiK55knZJgG93FMXdv8wUvMqpc06fAM3d6+UTP3ouLrEiTl48OI2UZkoXtaWz55C347vvM9HJJYsio5gJTslvsic9N0Te8ZkuyDHXrqQdhDk0AuGxZF1ZabfJ4JyNhUXvG7tp06dLazJ207pGiWeo2jLlT2Y6KZpjvA3Nryws8MgUvy2w9MgTdMDsZ8eSKFlF3IpWfPfOySLHi9FvgZZnLz55nk4qdFktf1JEW+k+2Hxv23wxQcID/u89rS+HA3bfbEXZ/eev5QtXArpw6YcpMo6Uq9h4f9R3XCLSgcfffYNIif/ueZ/Deb/JxrMTca8gyNZm7HOm4WrC7vIdkuNJQZyJaZqIqjpZRJOYyYmLmbr7op8fKoeUHAFiGtYYsY9cvFydSBYEia9540WJ7MUgrMjKq7PpuuZSC8XK1JnMnY8I7VTXdZLWKPHXmTijXYO6MMbzbKr4mioShz2SJCZmlF2++1EzbfHL3CfzGPz+Fg6fGwZjYkU4Gt6oZsGw7rl/Zbf9epLmXNB17Bkbwpn96Cj/ZVfGVh6Zn4mXu3sWFb29YDGDu/Pn4vxOFtX732YMATImVR1ATDgLdFm+tfwL5fGjx+dPvbcTr/uHxWDpE1ULLGXeKx+VXZ367un6vU/JG1GZPhLSAkfOsn1iV6DgeJ0dLwrZuBNJnw8ZDvwpr2t0ojJSNwJec2LimG6jqhm+Rnd+Wgiwx9A+XhM2ceWSU2osk1VthjNn1toPYNQ9JYtj8yTfg7992qf3ii1hpVpUxUdYDG3UQiDiQZmsYhsncJWbJVe75cFzQUDwMtkMzZMF/15Vn4ZEP34QPvmal73cUxTUvpwqJjxf/YIWiHhmcgGEAW44MubKSeZDRK3Oa+29fscxecL2JXyTL0Dvy9NGqKctwc0qSGFKy5Nq5iUIh57Wl8O33XgnA1Nz3WLsTr3H2QlUkV9kM5xq6fS3+HRU5s3k4DU/E1/XW69loZV57m7o3Ai1n3CkU7u1rl+IKS4/MBdx4eoYsNBiydvQG39Em6LhiRcPln3kEn/j51sDzEPsN695ChuifH9sTeEyjMFL2V3sESJYxxx60jZclhoX5NPqHi6E6LmA5VGv4FPjYb8M27rWNF2BmQaqyo9WL5BwzXrtqGxlRxT/AYYqU+ENzSpYk130BgPs2HcUVn30ELx2Kno7ubWYtgiQxrFyUF9ZAIedwGKngocoSujlWvKNvJDCc02HuOgzDAIO5CJAM6mXutLOlsQyWDGieBtmAea95h3pQIhHp3mPlKn604bCwTSZVs6RnLDEmfN529rSm458ecxqQBEXLEGhYQZKgLJmNxakAG/kBpiNrteWM+9kL2vDgn9+Av3rjBfjWH16Fdf9/wfdwiI3Zce61HKoCJima7mlFCozyIG85bf9EIKMXlMkVi18AACAASURBVFABOAYs7GVvFIYDmLvCOVRpcRNN9p7ODPpcxl1849NqNIcq/T09x6AEqyCEMXeKMrFlmSDmnnYzd/puisygyJKrKuR9L5ltEPbV0Fz5aCk6b6349CBcYHVXOj2JaA2Rdi+Cbdx1w5SiuKQmAHZJZEJ3Po3jIyXXbqYiCvP0yJtmM3X/M6J3gXax77jiLN8xn/5NMx/AMe7u31c0HQMjRTu0FgA2vGKya91qBynS3Am1mDtgznN6P8gH8PoLegKPjwstZ9wB00svSQzZlIxzutt8k4cmgx3nXuN8JiM3J+Tyu+7HX/33FmEz3jBH4NWff7TmuGmCjYXUf6ZJH7YARMWLB09j+V331yzSRYtfSQvQ3Ll47nf++7MAxL6PnvY0BoZL9r0MK2VQU3PnszZteW2Sxt16IUVjzagyJsq1HaodHlnGKYrGoEpuR/MJy1HujYv3gt+10FwNirOvBYpYOb+ntt5OEEXdiMDLMpqu25FUF1pt5vw5BubOjQ+HFIV58u8bIF4AzOPM59Zn1Zy//Ox5vmPSiozze9rt+SRLzM4PAIAv/moXrvys+92k9omUFxO6ayLNPUQS5KOERktVXLF8HhZ3hmcKx4GWNO5eeKvoUaJRlPIDgMMk6EX/3nMHhbJBkFbsjVsP0l3p/EExt6Ml08GXs1K8o261g/C4VfSod0d4P3O+SYIwWoZzqFLUgeiF6OkwmfuJ0RLyaSWQ7WRUCaUIzN2WZazPJlsAmbbIKUHdFLrHZBTSIaGQgGOwyZibJX/dce40D2rlQvCyHB3rdRRGRUqR8JM/vhbfft+Vkf/Gu0v444BOUjQvqprpZyHH6L/97lr84M6rfYliizsyOD1ewbiHmHjnivc9CqrMSO81LYBBWcR8fL3EGG46byHeeNFiAMBLBx2J7GO3md2/aNz07IKysul8ou/gvr4jM46VtCnvwiaLOWHcvTeTwr2ixrmnrSQm3iiLtPWgKI8Jj6HackSsudLWLYjZ0SQmo1SLAdYCsXB+jTAMA7d9+Um87V/X25/x90fE3E0j5v7eolu6uDODoYkKDp4ct7sEieDdAX32/pfxln952nVMiUsLp+cYsnsWgkIDRYlGWYu5O0lM4pPnPTVqKlwUjyIzV1Y0ZVHWymL+3APb7f8/OjQBYOrMHTAZ7aL26EyxPaPajlgAuOocf6MOwAleqGg6NN3xx3TmVFy1YoHveIrcOTI44T5PDVkmqAQAldg9Zt0jr8ZP4CVCiZm6/G9YUUGHTo/bv+vMqpAlZxdqG/cI0WthYdUpLmpqVFADqFGYnqvMMLw3k4oGkVGQIzJ3PkJFyNwDQiG9CQyHT0/4jgEc1hcU595vsbgVC9uw5cgQxkrBBa2igF5GjZuZ42UNL1vsmxqC8AY9KBXbayAHBUkbtB3efGQIq6zYbBEyimxHLcgSw78/ud93DN+dyJHXJsfdb7vkDBwdnMAtFovjQTVZShUNjAXHmbelFDDmhELa4awC5k7zQFSsjMePNjhNw/cfN3MfOqbI3KeCP33dSgxPVPG7XzclNm+9HgKfcGRq0+H3n2Sak55yu17WS2SKIIpzB6h+D7MXy6CyzPz5aS5TSOzh0xNYe/Y8vOmSM/Cbly3Bx3++1X5mtiwT6lAlghRs3fn4/tGSuPFLIzAnjLv3ZpJzk6SEMG84wBt3x+hSiVkeQZmV3tTjoDBGMgxBv+8fcTP3oKy6qKA17Wu9e/Hua5fjQ/+10aX3j1iLB7/4BTpUPRKRaGdDzO34SAnXnetndgQnBV3zlXYlCY2XZej4Ws/Ri46Mig+/4Xzh76gmy1ce2xMYCgiYkSr5lGKzcpdDldPcdd2wSYV3JxeGPZbzVVQ3plG4ZKk7kzWoDg/ld4yVTeMuyl7mQS0PvUW0RA5Vfh5WNCNQ9sgossPcA4y7m7mb/57Xk7eNbmdWxbuvOwcAXMXeHtzaByCcuVPdqjBnfkoxQ2JfOjSI4yOlhLnHCd6459OKY9xtTU1CmJnMpxWMlTVhU90f3Ok0zTXlBD+jJ8Z2w6pu7BkYtR1AXtCkCtq209+tWGiyjnplGd4B+Zn7t+O5A+6GFcNWoSWelcgC7YPayfE+gD+4ZrnvOBo3AFe4nRdOo2YdbdxhxYpu66w8c/8/t12ABfk0br3Qz8CnilxKtp/bWfPDa8u3Z5ziYbxDVZElO9lqtFy15a+wUFfALHhG0sW+42PIp5VJRwLFCSoG5gUZtvFS1QwZrKGL0Q7AZ9w9hjubkm3nMyAuP0DIpGSMlKpoTyuYF1Bzn3eYU0y9Ikt4/Zoe3L/5GI5xPhBZYtAsxr7JClm9YVU3gvCZt1yEcxe14dpzg49RLYfq03tMQnirYKfYCMwJzZ1fKfmOQKSt1WJ8VIui3+MIu+XCHpe2mFZllKu6L5KGFpO73rgaizoygcycxjNaqgoNN5VtJQbMb9+nAp6RixyYNG4+ukHE3GWrbjnJDZ940xohO1nSlbW14zA5KajMMqVwD01UsH7vSVsa68ql8L9vXR1r79qsKtuG+vet5h1BMLsxEWHgHaqOfvsAVy9ovIZxL1U1vPPKs2xjWI/eHgeyAcydxjduMfda6w+VWfbKlF7j3s5V2tStpLighYPGcPHSzsCewi5Zhpvzd926GgBcC4kqM3sX2jdcxKVLO0Mbx3Tn0/joLatrMvdHdwzgX9btQUdGwRWCZuONwJww7vzDNY27aSSIadbWCk1DdNhiU2daYUzeaA+avN5tN1XYa0spmJdTXcWVePD6rGgBGBgpYlFHGudZYW276sxy418GkcEZnCi7mmsA4sJhuZSC8UrVfkkW5MUMijFmL0xBPUmBYONORuGnL5qLGt/wOW7wBq3WNjqbUjBRcTIcAa6HqvXz01Zm9CVLO0Nj+A3DMNv65VTMs+7RjBv3AM1dkhiyqozxcrVmPDjgPFfv/Bd1WSLjXisckaawV0ri4XWoEpbOy+L3rj4bX/0dp+Wg+czMaw4Ml2o2N4kC2mGOlbVpCYEkzAnjzqMrp/pkmVpbXnLUHB2cQC4l23WbvYk6ZAS8WjhdJ6VIoan1FV23s+5Exr1vqIjFHRl0ZlW88aLFdVea40M0TwicuEMTFXzveXfSlWghpEbEUbr9kKHoDNhCA+Kyr4BZ7Ovbvz4grJsdN/jG2rUcYCmucJq9G7SYe0XT8f3nDuJ/Nh3FlefMx+rF7aHM3WxRZ6Azq9oS1FTDIONCLiSlvy0tY6ysmaGQNd4jOg/VWSF4DXdb2qkpzy+WItCcCyuKxp+fX38YY/j0b17kYtJ8cEDfcDEWY8xfP47FIirmjHG/5cIefOy2CzyyTHiNEwJJCEcGJ5BPK/bP3uQXu0Kgx7jz3WzCimJpmoEzu4KN+9BExdYV+Z6kUwUvtxw6Ne77/elxp08sQRSRQM0PKCY7bAKT4e4IYaMUU+5NCPvwDzfh4z/fhgHrhf6DGnJJPeDZai3mbjrhyLhT8TenKuRd/70FgCkxUEGyIPDtA8mROVPM/Z/fdRkuXtIZWm0zl1IwXqpC0/3dkvzHBjF399/l02a0VKmq2Ytm0Dv67uuWY8XCNmHoJYF/T2stQKqVVazrhut9qwe8wzjM1xQ35oxx/7ffW4v337gCXTkVuwdGsbt/xJ44NZm79YAPnTKNO0kKwczdbZToxafGEsHM3bD1PZFTVePqXHRmUxgar+DAiTHc++S+QKknDJWqw9zHODb5yd9YAwD42/u2+f5G1Fi5K2c2rPiPp82QxZ6QGPYoxatElQF5kPb+wdeuqnmuqYJPfIti3MlxSqxPtWrL8E7m/pGiFXkVnH1L3zmrypzmPjPM/U2XnIn/+dD1oUl+VKahqhu1m96oYs096D36/AM78PC2fgDBxv1/37oaj32kgPkhu7lUgCwjAjlUyZ/jTYCcCvjrh+V3xI05Y9wJnVa9i9d/6QlXZEMYiE0D5sSzmbvqnZTmRBgpuScvHx6XCambUtXMsKxcShZmsfJb386sirKm46M/3oTP3L8dj20PzzIVIagcLS1mfBnft12+FBIDrhYwJLof246a8fFBoXMA8N4bzJCzNWd0BB5j1/QOMIInR03jHseLFwSesdVizqos2d2lqtyz9hqkD71mFXIpGWVNFzZsB8BlxEr29wvb5cw0yLhHiXNPKxIk5vfvpGT3cyTj/o31B/CXP9lsHTN1UxWkuYtgNlgx7Hc0So+AWuDnwf/36qUhR8aLOWjcHRZE+qgovI9He0a1U5PHy04SgpeF5gOZuyP/BIVLAk7t6C5Bj0cArheIdg/PHzDrwohip7ceGQoNlwwy7l5H8Q//6Bp84bcuwb1vyGGZICwwzDnqxWvOX4QDd9+ORWHSjUING8SL4Aky7hGrQE4FF57pLD61mbuprZ8cLWG75eRWrTh3wkdvOR9vv2KZPeagWHe+xPBMM/coaEsrGLMcqrV2wGYDbf8zUz3lH0T5G95jJgNRElMQKNu6aNcUio+5v+/6c+xgiOnAnDPu/AsXpXYE4Ror6WbNmZ12JIU3Ky0f4FDlm1NQuKSoT2pFNzPxOnMpobPUZO7mI/OGEnpDGYsVDW/6p6fwx9/ZEPidvP0sCV4H2rL5WTDGAlkPX/2PGkfUA2KspnPRP8YToyWkFamhsd/z2lJY0pWFLDE7aiUIpLnf/pWn8OlfvAzAarPHMTZigPTcaPfhBd+zlW63twBXMyGXkjFe0lySYRhEC6WXlV8g2NXV8ouFQZTEFASzCJ7BNWmp30TSgtbInaYIc86482zQDoWMMCkvWtKJ5/7P6/DFt19qP3Av+8oHOFQpCYMxZv+tqHwBsZ/OrGKHa3p/bzN3j3H3GmqSdShxQoRyVUdHRvFtFfnF7uVP3YIzOrOB5wDgan7wiTetCT02Cvj7KLpPJ8dK0/KiPPwXN2LDX98cKjMBTlZjH+cEz6VkYZeoNdaOYMuRIRiG4XNk064uo0q2H6VnEnVhphsUBssTjzCINGev4b56xQJs/9SteO5jr7M/qxVmGQaXca/F3K0GK3Y10Ag+olqgBX26E9HmnHHnm/ZSskLUm76oIwNVluyV2KsJ0/bZ69zku7eLGgADZlhiRTOgSgxdWTFzr2i6qzgTD+9Y+iJ0aqpoOnIpBcvmu413R1a1HVS1DBvglmVqvTxR0M6V0R0XbNFPjpYDG7DEiba0ImzJ5kVK8bfTy6iyyyCRcT9/cTtUmeHlY8P42uN7ccMX1rnqu5c45k4hg9MZGz1Z8Mw9SvUHUWNpUd2cbEp27QhFlTujwpWhGsmhyjP3+ueZqBXjdKB593sNws1renDNigXYe3wUGoVZTZIV0AP3ppFnUzLyacXnDOULHzkJOm5j4OwiJHTlVLxy0h+ayDN3ryzjZe5R2vBRzQ5a8K49dwE+95aLsby7Db0fLdhOwlqop3iZCBlVRko2meuAwLE8XtbsDkPNAD4UkuBn7pJ9bE9HBn1DRTywxfSX8I00eOZOWa9h0Uczjba0gqGJCha2pyMyd/dzO2t+LrCao+v+1cGg+cCHWkUCFZlhomLY0mocDlXa2U4loq0ezDnmDgCrz2jHhBW+BYizLsNwtuVUPGuB37m4yGpIwIMvWcoXxeLBFzHrzKkoVXX8z6ajvmNorF5G6ZUvBmqUlQWsDjcyw0VWc4Wz5uewvNus/9KRUbEgYkxu3MYdMHXmkWI1cJGKsqOYLpAswzvucinFJTfwxqnHKkFBCzg/F5woDdnWnhfmm2ch8yKXMss07OgbieS78kpMYan4fAhmUJvDKOA1/VrrD/UD/p17raYzMTD3FdY7tSykjEEjMCeNey4lY7yi2VEsUSYlj2tXduP7d16N/3WTv4nB4o6ML0bd7CTjbsflNcYVbhfxvutXAAB297vLC/DMnW8E3JlVbactwa7NYUDYNQpwdhRL5+Xwkz++Fh+fol5ej7MrCO0ZFcNc1us33nOF6/fdASUOZgJmnLvuigRJK5LLl8Nv7xd3ZFw1zXkfDckyaVXC37/tUvzsA9cFMttmAO98j0KSvM7hoFIVXkTJjwj8W5V3qNYOheTLgMThUL1qxQJ8/86rA5ueNApz1Lgrlq5mvkhTcdZcvWKBUKunLTdg1o1fftf9+MXmY3bkRJAzVuP0/4XtaSxqT7sWCcOg+ubm37tYjeJvzM2f/7avPCn8DnzZ3MvPnjdtpUijoD3jMHfGgOtWuqvuTWcady1QKKQqux13bs3dncjCy258dBWv9ebTCl61LLhmSjOA5w21JA/AHy1DrLYW6tG++Tj6WmM8MVrCDq5mUxyaO2DaizgL20VBpKsxxm5ljO1kjO1hjN0VctxvMcYMxtja+IYYP+iBUTu7yTL3MPR0ZDAwUoRhGNh6xEzqGS9rtn4436ph623I4a2hQS3pCKIiZ999n9kAPCVozM37A3YEFBgL6k05FfzgzqvxyIdviuVcABn3CoYmKshbEsd333eV/fvmMu4SDAN2qVjnczFzn++R1NzGPb4ojekA33gkyntEkVCL2tN470UpvH3tskjXqYdBu5h7jTG+eNDdJS0u4z4TqHnHGGMygK8CeCOANQDeyRjz7d8ZY+0A/hTAs3EPMm6QpPGjDYcgsXgiPAg9HWlUNAOnxsp2mjzgSBfkHHvvN1+wNfXHdx3H9X/3GACnMQC/AwDEjUWuW9mNc7rbTOauBRt3APj6U/5uRpWqEZukctWKBVgZ0l1psmhPmyUNxkpVu1EFz96nM427FugeDnuiIURx7oATXUSsdYQz7hMVDRKrnTXdLOBlDlHYqhdk3GWJ4YalauR3rz7mbt77qbzmcRK/6UaUN/tKAHsMw9hnGEYZwPcB3CE47tMAvgCgdpjGDKNw/iIAJnOvJ35WBKeqY8kVzkjyB1846EPf2wjAjEWnF4Mm07kL27DvxKgtt4SVJ04rso+5j3tkH0qu4VHWdKgxRAM0AiTLjJXFPScXtDWTcbfaFXp8Gyr3rHLcd6DoiZSViMUz933Hx7B0Xq5m0/Zmwf+66VwsnWeG0kaJBqFnOdnm7vVErRBzjxLy/LMPXGf//5KurN0acjYiyh1bAuAQ9/Nh6zMbjLHLACwzDOMXMY6tYVjcmbF7eMadWEChXk/sPu6qqEjszsuU//2Jfa6IEPr9JUu7UNEM7LQkFTuyR7AYjZaqeGR7P97zn8/hoz/ahJ19IzW7/QDUzag5jUh7RsVIsYLRkiY07pMpedBoBNUa55l7novuoUgnxhjaUrJdrmJn3wge3NaHS0LK1zYb2tKKnbg2EiGOm5h7WM9REeJh7rXn+quWdeH2i88AAHzotStnzSIrQhQPmujb2U+GMSYB+BKAd9c8EWN3ArgTAHp6etDb2xtpkI0Aq1jRCoaG3t5ejI6OxjKevjGTQd/9yx2uz8dGhoTn/+wD27F6vjn5Ll0oo3JsJ3qHdmNw1Org8+QLOHWGgpGyecv3792D3uorrnMctLIc1+08bv778hEsyjGc0ylhYZbhuT7TeKxbt841WQdHxpHTxyJ/77juURSc7CtjrKzhcP9JZBT4rrt72yYUDzbHrmPfIYexruyScPVCc069MuwssC8885RNJPYOmp+Pjo5CgYE9rxxGb+9xPHTAPM858qkZfTcmi/2nzO8zcEo8x3mcmDDndbFUxuhoKfL3fPrJJ6Y8vpPWNQ1Dj3S9E8dNsrVvzy70ju+b8nXjQD3vXBTjfhgA7/VYCoAPwG4HcBGAXstwLAZwH2PszYZhvMCfyDCMewDcAwBr1641CoXClAYdB77zygvYdbofmZSKQqGA3t5exDGegZEi8OSjvs8XLpiPQsF0CHY/9YirOcaOUzpuv+QMfPVdTkeYoYkKPvbUwxjN9mAHy+Gt1ywBHnsUq88/D4WrPXXMH7zf9WPZkJFpy2NBWsG333sV/v2JffjsA9tx+TXXu7IBU8+vw5IzulAoXBbpu8V1j6Jgj7wPP9+7HRNIYdXiThQKlo/e+q4333hNzbII04UTGw4D2zYBAD54yyXoGtqNQqGAQ6fG8Tfr1wEAXvfa19jHLzg8BDzzFNra2pA1DLTPy6NQuBwvPrwTbOce/MXbXherH6jROKNvBJ9/7gnocrrm/Dg1VgYe/xUkWUE+X/t4et71zLsToyXg8UeQUpRI5/lZ30ag7yhWr16NwjRWcRShnncuCvV5HsAqxtg5jLEUgHcAuI9+aRjGkGEY3YZhLDcMYzmAZwD4DHuzgRJv4g5Pak+L5QJ+W3nvH6z1yUHetOyOjIKsKuM7zxzE3b/cgSEri1GkuX/qjgtdP4+VNYyXNTvumjoXDY55SxHH51CNG7QIHRsqimWZbPPEufN6MB+THiQdrerJ46pz5uPzb73Y1VJu0GpIPpsMO+B8zygN2zuzKq5f2e1qbddokGwW9baSfDNZv0CzoeabbRhGFcAHATwEYDuAHxqGsY0x9inG2JsbPcBGgSZk3N7wjCrZk+h915+DcxeaERF8oa9XLevCT//kWtffeVPMzX6jzmcU/yzyEfz+Ncux49O32j9ruoGJimYX16LJXdbM7fPxkRJ+suEwSlW9aY07X0Nd1OYujuSSuMCPj3/OQe35MqqMH/zRNbjsrHnIpxVsOjSIR17ux+B4pSHZvo0GjTmKcZclhu+87yrcsGpho4dlgxbfqIsmHTdJt0DTIVLWimEYDwB4wPPZJwKOLdQ/rMaDjEcUJ8tkwBhDSjFb6XXlVJuxe1ncioXusEFR3PaijgwOWMku1JotqIKl1+E0VqraKfo0uSki58++vxHrrYbNzepQXdXTjo6MgmJVd3V/+oubz8NPXjzcVI4ufmfRlUuB2nZHGWNbSsFwsYoP/NeLWLt8nrCNYbMjo8poTyv48BvOi/3cv3/N2djwyum6zkEO1ShJVgDwO1edhR9vOIxrVwa37psNaJ6UxGkGpU03wkbQdq4zq7o6J/HIpxV86o4L8Ymfm63sROnVvFRDGadRijMBpmbvbQdIxn0nl9TUrMx95aI8Nn/yFt/nf3bzKvzZzY1rrzcV8Ax9ssybFoZSVcfz+0/jqhXBtVaaGVv+1v+s4sCn7rio7nMQ4YrK3C87ax4O3H173dedacxZ455tYOEpClvszKXsjMNOQenYtWc7LzLf+YewkIuxpdDGqDJSRTNsFmjLMpZx53ebQWF8CaKDN+7elnjd+VRoNi1fzbOs6bE0ZE7gR1qWppTENJsxd417A9OKSavrzKq2Vi5idGvO7MC+z90GQKwH8r0m6f8nE5dv93pVnGJlVU3H6XF/5myCqYN65wJ+B/3zH7s5VJ7Zf2LU9XMzl/edzUirUmRZplUwZ99sbyu5OHHlOSYjX9KVsRtBL+kSszdJYoHbxVef5RSNGq/UZu6UfEHwyjLlqo4To2WXoyhh7vUjrOBaLd39unPNkgrLrfLRzVQzp5WQkqWm8tNMB+Yuc2+gcf/6H6xF/3AJKxfl8cXfvhQfPn2enRE7GfzW5WaM7Ud/vBkTlkM1jLn/4ztehcvO6sJn7t8OwElzdzR3zdehabbUMGlm1JMa/9Fbzsf7b1yBu36yBQdOjkfq/JRg8kir8qwPbZws5ixta6Qs055R7SJauZSC83rap8QaGGNYZjUGGS/XLk9MtdkJpN/ymvuAZdypaFUiy9QPerZnzZ98MwZFltCdT+Om80wGv3AW1zJpZqTkxjZUb0bMWeZOskyzx7KS8XWiZcInKB8fvsKKsec1dzNtAXjVWV3Yd2IMx4aavs7brMBjH7mprmJmv3v12bhwSSdefda8GEeVgJBWJVT0aG0jWwVzlrY1krnHCdry29EyNWQUcqK2ZxTbqBNzL1U0jFnyDlXGjDvOf65ixcJ8XR2TGGOJYW8g0srcc6jOWeZOBq/Zn7ddK9wqH1yricOaMzrw+bdejIu5xB/boarpqFbMrcob1vTgnt+7HNev6haeJ0GCVgKVV55LmLPGnYxms2uctAhRbfhajmBJYnjnlWe5PrMdqhUdxaoGWWJIKxLecOHiBow4QYLmQ0aRE+M+V3BmVxaff+vFeN3qRTM9lFB4jftUQjgVK4GjrOkYLVbRlpLnXFhYgrmND752ZaR6862EOWvcAfgYbjOCQhVt5j5FX0FakfGjFw6jb7iIMzuTWOoEcwuXzUF/xpx1qM4WpK3O7UMTJuuYanx+W1q2Y9zzmTm9pidIMCeQGPcmh6oQcy+DsaknzLRzTTrCMioTJEjQGkiMe5Mjo8jIpxVUNAM5depaeZFrmK3G3BQ8QYIEzYfkLW9ySBLDRUvM+jTjldpNr4PAG/dMA0svJEiQoDmQGPdZgP/X3r0HyVXWaRz/Pkm4BBKIQETuLIElMchlg8GgWSIM6Ja7EF1AClDuyqIBwi6yXAKKLGt5W5USBDYCgogQoAxaZQB1uQQEEiTLAgEDZDdaqCDLkgABQp79432b6QyZ9JnpmT49p3+fqq7p6fN29ztPTf/O6fec856uCVsCzZ1NW5udEmDDmCwshMqLT/kQcOLUnZp+jZWr6rbch8jZuSGE/oviPkRM3nEz9thuTOOGvTijq/sSaIfvvd1AdCmE0MbisIkh4saTpzT1/BkH7MKMA9rr8nQhhMETW+4hhFBBUdxDCKGCoriHEEIFRXEPIYQKiuIeQggVFMU9hBAqKIp7CCFUkFzSFaIlPQ/8dylvvnZbAC+U3Yk2FxkVEzkVEzk1traMdrA9ttETSyvu7UbSAtt7l92PdhYZFRM5FRM5NdZMRjEsE0IIFRTFPYQQKiiKe7cryu7AEBAZFRM5FRM5NdbvjGLMPYQQKii23EMIoYI6qrirvxcg7SCRUd9EXqFddVRxB+ISRI1tULsThat3kv5F0gTHuOZaKWn+EmKh3zqiuEs6StJ9wFclnVR2f9qRpCMkLQa+JekMgChc7yTpSEl3A6cAR5fdn3YkaTgwD/i+pIYn23QqSZ+R9GVJjOl1KQAADAFJREFUIwfj9St7Jaa81bkxcBGwG3AOsDlwnKQ7bC8tsXttRdIOwKnA8cD/AnMkvWD7B+X2rH1I2gT4GrAjcDYwAdg0L1OsCNcwAliftPH4IUm32V5Vcp/aQq5LI4ATgbOAlcDtwD0D/V6V3HKXtL6TFcANtrts/wepcP0BWFZqB9tE3bDLRsCTwGO2nwBOB/5R0malda7N2H4ZuNL2R2zPBwwcnpd1dGGXtF7d/WG2XwduA24FTgDeXVbf2kldXXoTeJi0gXA5aYNz84F+v8oVd0kXANdLOk7SFrZ/ncf/PgFcCewMfFPSp3P7ymXQiKRzJO1TV5RGAGNJRR7bdwBPAV/I7TsuI+jOKd8fbntB3eKbgVWSdi+nd+1B0tnAbEnHSBpre7WkrYEu4NvAc8DhkqZLGl1qZ0tUV5eOlbSZ7QdsvwZcBmwLdA3056xSH1pJM4EPkQLbH5glaatcxH4PvN/2NOAW4CxJY2yvLq3DLSZpK0k3k4r2dbXHbT8KvAJ8tq75PwNHdFpGsPacbL/Vo9m7gGep2GeoKEnj836sicBNwGGkIr4h6dvxg/n/ZhnwFeDzQM8MO0KPunQAcIGkrQBsrwSuAo4kDfkNmMr8Y+adOHsBX7L9C+DLwKvATIC8pnwxN38SWEjeUu0g/wfcZHsM8FJtx2l2AfBxSXsD2H4auBMY1fpulm6tOUl6ex+V7WeB7YE987LKfJYKWg7caPto27eRNpim5GK1NXC0pLuAjwJzgQdJ48sdZR116fRaG9s/Al4G9pP0fklHDcR7V+IfMu/Qegv4I2mMD2AJ6R9uV0mT6tsCs4D1gD+1uq9lsv0q8LP860zgXEnr52WLgWuAUySdJekyYBzw51I6W6LecrK9StKw/IGFtMV6YH5OR327sf170jBnzQPAppI2sv074EZgru19gWNIK8HtWt/T8jSoS+Pr6xLwA+DSvGzDgXj/IVncJU3ORy8Aa+zQuhzYVtKk/GFbCjwE7JGf9zngEeAN4Pgq78HvmVGN7eX5n+5e4C7ge3WLv51vWwMrgL/N44KV1decbK+uG6J5A7i16ucDrCOjV+p+3R9YlleM2D7b9jfy/deAg2230/UbBpykgyWNq/1eoC7VvvXtTNqivw7Y1fbsgejPkCrukvaT9DjwGWCTusdrf8f/kIYSvgBg+3nSnvra8gXAobZnVrVorSMj1RWh2pbnPwDTJW0haSIw3vYi4AzbZ/b48FZKEzlNkPS+/Pj3bc+p6tEyRTKqG6oaB9ydH9tL0va1tgAV35DqknQ/MBvYqu7xRnWp5kVguu2TaivHgTBkinveUXMacKHtE/NXv9pRDLWvxJsC1wKbSzovr0V3Bd6Et8fdf1tC91uiQUa2baWTSkYA2P4j6Wvgn0g7dciPV3rHV5M5XQOszo+/Wcof0AJ9yKh2RvPGwFhJV5HOLRkO1T1MNK/fRkm6DTgv334N7JCXjyhQl1YB2H7R9h8Guo9DprgD2wB/tn2DpJGSPpH/uYYBSLoU+HfSB+90YCTwY2C+7WvK6nSLNcrou6Rx0p3y2PGnSIesnWV7su3Hyut6S0VOjRXNaGtJW5LO1v0ssMj2x/IO58rK67cVwA9tT8s7S+cBh+TlqySNkHQJJdWltj1DVdKppLHfBbbnkLa+PyxpKmkt+Rrw98BiSVeQjuo41nbtRKVzJV3odEJFJfUjo9HAcTkjJD0B7Gn7pVL+gBaJnBobgIzOAa6pOyKtkupyWmj7Jts35MeHkU6SXCZpg1x3dqPMumS7rW6ASEcozAcOBZ4ATsjLvgEsBrry7+8F/gvYpe75w8v+G4ZARiPK/hsip/a4xeetqZyOBcbWtdkXWNzL81ueU9ttudu2pA8D59n+laQVwN9I+iTwLWAG3WOhj0u6lzTeVzv1udLjxTAgGVV251a9yKmx+LwV00tOHyGdmHVtbnOfpN9JOsT2T6Q051BZOZU65l53VELt9/qjWqYC2P45aethEulA//OAMyRNlDSL9NVnWW5buWONI6NiIqfGIqNi+pDTU8BESeNzu01I2b2R2zj/LCWnsneorjHVZV0IS4DR6j7k7C7SHuedbH+VdDzo50jzxBxmu8on2kRGxUROjUVGxfQ1p1G53cukeWK2bFE/16mU4i7pA0pzd3xX0kHKZ/yp+5jZB0lfdw7MhxQ9Tjp+dAqA01S0p9k+xvZzJfwJgy4yKiZyaiwyKqafOW0D7F33MkfYvrqV/e5Ny4u7pGl0n2b7JOkQqnfVj3HaXkI6g2tn0gRWAK8Dz9Rex9U+xngakVFDkVNjkVExTea0tPY6TnPrtIUyttx3Bx6y/UPS1731gBW1rz6SLpI0mzSx13eAyZIWks7iur2E/pYhMiomcmosMiqmcjkpj/kP3htIHwBetP1U/n1P4JekOUxOIR1S9AzpBIBlpKlBz89rSSSNIh2SVuVjjCOjAiKnxiKjYjohp0Hbcpc0RtLPgDtI8zzXdjo8QpoGdAfgFKf51eeTzgB8xfaRtpfU9lDbXtHOATYjMiomcmosMiqmk3IazGGZjUlrvRn5/tTaAtsPkq78U5sl7pfAGNIZXrXjZyt5mFUPkVExkVNjkVExHZPTgBZ3SZ9WmkluE6f5nq8gzeu8EthH6fJbSNoAuI/09QfS1Uk2y+0qe/wsREZFRU6NRUbFdGpOTY+5SxLwHuB60uQ4T5PWiKfZfiG3+SDpYsILbF+bH5tIuvrPe0jzWHze6eLMlRMZFRM5NRYZFRM50dzcMuT5EoC/BK7L90cAlwC39Gg7kzQV6BhgZH5sJOlEiZbOudDKW2QUOUVGkVMZt34NyyhNZXkxcLGk/UhzE78Fb0/KfyowJS+ruZJ0JtcdwFJJ29h+zfYzVFBkVEzk1FhkVEzktKY+F/cczELS1d+XkC4PVZsedDK8PafChcAX6576MdJY1iPA+5zGviopMiomcmosMiomcnqn/swKuRr4urvHqPYC/gI4H7gMmJQPF7qVFOyOtpeSdkp02b57QHre3iKjYiKnxiKjYiKnHvozLLMQuFHdV4CfD2zvNJ/CcEkznPYqbwu8lQPE9k+qGGAvIqNiIqfGIqNiIqce+lzcbb9q+3V3z098IPB8vn8cMEHST4EfAQ/DO6fQrLrIqJjIqbHIqJjI6Z36fbGOvIY0aXrLufnh5cA5pDmfn62NX+Wxro4TGRUTOTUWGRUTOXVr5iSm1aTJdV4Ads9rxVnAatv3VmnHRBMio2Iip8Yio2Iip6ypk5iUJt+5L9+usj17oDpWFZFRMZFTY5FRMZFT0mxx3xb4FPBNt+Jq3kNQZFRM5NRYZFRM5JQM+pS/IYQQWq/sa6iGEEIYBFHcQwihgqK4hxBCBUVxDyGECoriHkIIFRTFPVSGpC9K+qd1LJ8u6b0FXmeNdpIulNQ1UP0MoRWiuIdOMh1oWNx7trN9vu07B61XIQyCKO5hSJN0rqQnJd1JujgDkk6S9JCkRZJulrSRpH2Bg4GvSXpE0rh8+7mkhZLukTS+l3ZXSzo0v/ZSSRdLul/SAkl/JWmepKclnVzXrzNzH/5T0pdKiCZ0uH5PHBZC2SRNAo4A9iL9Lz9Mmvr1FttX5jYXASfYvkTSXOCntufkZb8ATrb9W0n7AJfa3n8t7Xq+9TLbUyT9G3A18EFgQ+Ax4HuSDgJ2ASYDAuZK+uuqTi0b2lMU9zCUTQVutf0qQC7KALvloj6GdAm1eT2fKGkUsC9wU13x3qDg+9be51FglO3lwHJJKyWNAQ7Kt9/kdqNIxT6Ke2iZKO5hqFvb/BlXA9NtL5J0LDBtLW2GAS/Z3rMf71mbr2R13f3a7yNIW+v/avvyfrx2CAMixtzDUHY38HFJIyWNBv4uPz4aeE7SesBRde2X52XYfhl4VtJhkC7cIGmPnu36aR5wfP52gKRtJL27idcLoc+iuIchy/bDwI9JFze+GbgnL5oFPEC6ov3iuqfcAJwp6TeSxpEK/wmSFpHGyw/ppV1f+3U7cD1wv6RHgTk0t7IIoc9iVsgQQqig2HIPIYQKiuIeQggVFMU9hBAqKIp7CCFUUBT3EEKooCjuIYRQQVHcQwihgqK4hxBCBf0/HOR5AfqTrw0AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Hs.plot(grid=True);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In the above example, we used the OPeNDAP server to get the quality controlled NDBC data.
\n",
"\n",
"If we wanted near-realtime data, we could get it the same way...
"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"url='https://dods.ndbc.noaa.gov/thredds/dodsC/data/stdmet/'+buoyID+'/'+ \\\n",
" buoyID+'h9999.nc'\n",
"ncdata=xarray.open_dataset(url,decode_times=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"or another way
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Unidata has been developing Siphon, a suite of easy-to-use utilities for accessing remote data sources.
"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" wind_direction | \n",
" wind_speed | \n",
" wind_gust | \n",
" wave_height | \n",
" dominant_wave_period | \n",
" average_wave_period | \n",
" dominant_wave_direction | \n",
" pressure | \n",
" air_temperature | \n",
" water_temperature | \n",
" dewpoint | \n",
" visibility | \n",
" 3hr_pressure_tendency | \n",
" water_level_above_mean | \n",
" time | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 80.0 | \n",
" 4.0 | \n",
" 5.0 | \n",
" 0.8 | \n",
" 12.0 | \n",
" 9.1 | \n",
" 100.0 | \n",
" 1017.9 | \n",
" 27.5 | \n",
" 29.0 | \n",
" 22.0 | \n",
" NaN | \n",
" 0.9 | \n",
" NaN | \n",
" 2018-09-11 15:50:00 | \n",
"
\n",
" \n",
" | 1 | \n",
" 90.0 | \n",
" 3.0 | \n",
" 4.0 | \n",
" 0.9 | \n",
" 11.0 | \n",
" 8.6 | \n",
" 97.0 | \n",
" 1017.9 | \n",
" 27.6 | \n",
" 28.9 | \n",
" 21.9 | \n",
" NaN | \n",
" 1.4 | \n",
" NaN | \n",
" 2018-09-11 14:50:00 | \n",
"
\n",
" \n",
" | 2 | \n",
" 10.0 | \n",
" 3.0 | \n",
" 3.0 | \n",
" 0.9 | \n",
" 11.0 | \n",
" 8.1 | \n",
" 101.0 | \n",
" 1017.4 | \n",
" 27.8 | \n",
" 28.8 | \n",
" 24.3 | \n",
" NaN | \n",
" 1.8 | \n",
" NaN | \n",
" 2018-09-11 13:50:00 | \n",
"
\n",
" \n",
" | 3 | \n",
" 320.0 | \n",
" 3.0 | \n",
" 3.0 | \n",
" 0.9 | \n",
" 11.0 | \n",
" 8.2 | \n",
" 89.0 | \n",
" 1017.0 | \n",
" 27.8 | \n",
" 28.8 | \n",
" 22.9 | \n",
" NaN | \n",
" 2.0 | \n",
" NaN | \n",
" 2018-09-11 12:50:00 | \n",
"
\n",
" \n",
" | 4 | \n",
" 340.0 | \n",
" 3.0 | \n",
" 3.0 | \n",
" 1.0 | \n",
" 11.0 | \n",
" 7.8 | \n",
" 88.0 | \n",
" 1016.5 | \n",
" 27.8 | \n",
" 28.8 | \n",
" 22.7 | \n",
" NaN | \n",
" 1.7 | \n",
" NaN | \n",
" 2018-09-11 11:50:00 | \n",
"
\n",
" \n",
"
\n",
"
Next, we need some model data for the same time period.
\n",
"In this example, we are going to use the archived monthly buoy files:
\n",
"NCEP_1806.tar.gz
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"MemberName:
\n",
"\n",
"- NCEP_1806_000\n",
"
- NCEP_1806_024\n",
"
- NCEP_1806_048\n",
"
- NCEP_1806_072\n",
"
- NCEP_1806_096\n",
"
- NCEP_1806_120\n",
"
- NCEP_1806_144\n",
"
- NCEP_1806_168\n",
"
\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"ncep_cols=['id','year','month','day','hour','u10','udir','Hs','Tp']\n"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"w3data=pandas.DataFrame()\n",
"tar=tarfile.open('NCEP_1806.tar.gz')"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"for fcst in range(0,169,24):\n",
" memberName='NCEP_'+str(year)[2:]+\"{:02n}\".format(month)+ \\\n",
" '_'+\"{:03n}\".format(fcst)\n",
" member=tar.getmember(memberName)\n",
" f=tar.extractfile(member)\n",
" frame=pandas.read_csv(f,names=ncep_cols,\n",
" sep=' ',\n",
" usecols=[1,3,4,5,6,7,8,9,10],\n",
" skipinitialspace=True,\n",
" index_col=False)\n",
" frame['datetime']=pandas.to_datetime(frame[['year','month','day','hour']])\n",
" frame=frame.drop(['year','month','day','hour'],axis=1)\n",
" frame=frame.set_index('datetime')\n",
" frame['fcst']=fcst\n",
" w3data=w3data.append(frame,ignore_index=False)\n",
"tar.close()\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" id | \n",
" u10 | \n",
" udir | \n",
" Hs | \n",
" Tp | \n",
" fcst | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 00:00:00 | \n",
" 13130 | \n",
" 8.4 | \n",
" 21 | \n",
" 1.4 | \n",
" 13.6 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 13130 | \n",
" 8.2 | \n",
" 26 | \n",
" 1.4 | \n",
" 12.9 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 12:00:00 | \n",
" 13130 | \n",
" 8.3 | \n",
" 26 | \n",
" 1.4 | \n",
" 12.5 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 18:00:00 | \n",
" 13130 | \n",
" 8.5 | \n",
" 25 | \n",
" 1.4 | \n",
" 12.0 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-02 00:00:00 | \n",
" 13130 | \n",
" 9.7 | \n",
" 25 | \n",
" 1.5 | \n",
" 5.9 | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
So now we have the NDBC data for buoy 41008 for June 2018, and the WW3 buoy data for all buoys for June 2018.
\n",
"Our next step is to subset the w3data to just the data for buoy 41008.
"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"buoyID='41008'\n",
"model=w3data[w3data.id==buoyID].copy()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" id | \n",
" u10 | \n",
" udir | \n",
" Hs | \n",
" Tp | \n",
" fcst | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 00:00:00 | \n",
" 41008 | \n",
" 6.9 | \n",
" 178 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 5.7 | \n",
" 227 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 12:00:00 | \n",
" 41008 | \n",
" 4.8 | \n",
" 277 | \n",
" 0.7 | \n",
" 7.3 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 18:00:00 | \n",
" 41008 | \n",
" 3.1 | \n",
" 182 | \n",
" 0.6 | \n",
" 7.3 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-02 00:00:00 | \n",
" 41008 | \n",
" 7.5 | \n",
" 190 | \n",
" 0.7 | \n",
" 7.2 | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
If interested in a single forecast:
"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"m000=model[model.fcst==0]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" id | \n",
" u10 | \n",
" udir | \n",
" Hs | \n",
" Tp | \n",
" fcst | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 00:00:00 | \n",
" 41008 | \n",
" 6.9 | \n",
" 178 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 5.7 | \n",
" 227 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 12:00:00 | \n",
" 41008 | \n",
" 4.8 | \n",
" 277 | \n",
" 0.7 | \n",
" 7.3 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-01 18:00:00 | \n",
" 41008 | \n",
" 3.1 | \n",
" 182 | \n",
" 0.6 | \n",
" 7.3 | \n",
" 0 | \n",
"
\n",
" \n",
" | 2018-06-02 00:00:00 | \n",
" 41008 | \n",
" 7.5 | \n",
" 190 | \n",
" 0.7 | \n",
" 7.2 | \n",
" 0 | \n",
"
\n",
" \n",
"
\n",
"
Recall the NDBC data
"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" lat | \n",
" lon | \n",
" udir | \n",
" u10 | \n",
" Hs | \n",
" Tp | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 01:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 172.0 | \n",
" 8.2 | \n",
" 0.80 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 02:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 178.0 | \n",
" 7.5 | \n",
" 0.89 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 03:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 186.0 | \n",
" 7.0 | \n",
" 0.80 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 04:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 200.0 | \n",
" 6.9 | \n",
" 0.86 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 05:00:00 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 205.0 | \n",
" 6.3 | \n",
" 0.79 | \n",
" 3.0 | \n",
"
\n",
" \n",
"
\n",
"
Now we can merge both the model and the NDBC data for the same buoy:
"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"both=pandas.merge(model,buoy,left_index=True,right_index=True, \\\n",
" suffixes=('_m','_b'),how='inner')"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" id | \n",
" u10_m | \n",
" udir_m | \n",
" Hs_m | \n",
" Tp_m | \n",
" fcst | \n",
" lat | \n",
" lon | \n",
" udir_b | \n",
" u10_b | \n",
" Hs_b | \n",
" Tp_b | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 5.7 | \n",
" 227 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 6.1 | \n",
" 217 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 5.3 | \n",
" 228 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 48 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 6.0 | \n",
" 214 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 72 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 6.6 | \n",
" 208 | \n",
" 0.8 | \n",
" 7.5 | \n",
" 96 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" id | \n",
" u10_m | \n",
" udir_m | \n",
" Hs_m | \n",
" Tp_m | \n",
" fcst | \n",
" lat | \n",
" lon | \n",
" udir_b | \n",
" u10_b | \n",
" Hs_b | \n",
" Tp_b | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 06:00:00 | \n",
" 41008 | \n",
" 6.1 | \n",
" 217 | \n",
" 0.9 | \n",
" 7.4 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 216.0 | \n",
" 5.6 | \n",
" 0.74 | \n",
" 4.0 | \n",
"
\n",
" \n",
" | 2018-06-01 12:00:00 | \n",
" 41008 | \n",
" 5.0 | \n",
" 275 | \n",
" 0.8 | \n",
" 7.3 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 267.0 | \n",
" 5.7 | \n",
" 0.73 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-01 18:00:00 | \n",
" 41008 | \n",
" 3.4 | \n",
" 164 | \n",
" 0.6 | \n",
" 7.3 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 161.0 | \n",
" 1.5 | \n",
" 0.61 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-02 00:00:00 | \n",
" 41008 | \n",
" 7.7 | \n",
" 179 | \n",
" 0.7 | \n",
" 7.2 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 174.0 | \n",
" 4.4 | \n",
" 0.53 | \n",
" 7.0 | \n",
"
\n",
" \n",
" | 2018-06-02 06:00:00 | \n",
" 41008 | \n",
" 6.9 | \n",
" 233 | \n",
" 0.8 | \n",
" 7.2 | \n",
" 24 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 217.0 | \n",
" 4.3 | \n",
" 0.62 | \n",
" 3.0 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(dpi=100)\n",
"ax=plt.axes()\n",
"both[both.fcst==0].plot(ax=ax,y=['Hs_b','Hs_m'])\n",
"ax.grid(which='both')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"So we have both the model and the NDBC data for the buoy
\n",
"\n",
"in the pandas DataFrame \"both\"
"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Index(['id', 'u10_m', 'udir_m', 'Hs_m', 'Tp_m', 'fcst', 'lat', 'lon', 'udir_b',\n",
" 'u10_b', 'Hs_b', 'Tp_b'],\n",
" dtype='object')"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"both.keys()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Ready to calculate some basic stats...
\n",
"First, see what methods the object has built-in
"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"['Hs_b',\n",
" 'Hs_m',\n",
" 'T',\n",
" 'Tp_b',\n",
" 'Tp_m',\n",
" '_AXIS_ALIASES',\n",
" '_AXIS_IALIASES',\n",
" '_AXIS_LEN',\n",
" '_AXIS_NAMES',\n",
" '_AXIS_NUMBERS',\n",
" '_AXIS_ORDERS',\n",
" '_AXIS_REVERSED',\n",
" '_AXIS_SLICEMAP',\n",
" '__abs__',\n",
" '__add__',\n",
" '__and__',\n",
" '__array__',\n",
" '__array_wrap__',\n",
" '__bool__',\n",
" '__bytes__',\n",
" '__class__',\n",
" '__contains__',\n",
" '__copy__',\n",
" '__deepcopy__',\n",
" '__delattr__',\n",
" '__delitem__',\n",
" '__dict__',\n",
" '__dir__',\n",
" '__div__',\n",
" '__doc__',\n",
" '__eq__',\n",
" '__finalize__',\n",
" '__floordiv__',\n",
" '__format__',\n",
" '__ge__',\n",
" '__getattr__',\n",
" '__getattribute__',\n",
" '__getitem__',\n",
" '__getstate__',\n",
" '__gt__',\n",
" '__hash__',\n",
" '__iadd__',\n",
" '__iand__',\n",
" '__ifloordiv__',\n",
" '__imod__',\n",
" '__imul__',\n",
" '__init__',\n",
" '__init_subclass__',\n",
" '__invert__',\n",
" '__ior__',\n",
" '__ipow__',\n",
" '__isub__',\n",
" '__iter__',\n",
" '__itruediv__',\n",
" '__ixor__',\n",
" '__le__',\n",
" '__len__',\n",
" '__lt__',\n",
" '__matmul__',\n",
" '__mod__',\n",
" '__module__',\n",
" '__mul__',\n",
" '__ne__',\n",
" '__neg__',\n",
" '__new__',\n",
" '__nonzero__',\n",
" '__or__',\n",
" '__pos__',\n",
" '__pow__',\n",
" '__radd__',\n",
" '__rand__',\n",
" '__rdiv__',\n",
" '__reduce__',\n",
" '__reduce_ex__',\n",
" '__repr__',\n",
" '__rfloordiv__',\n",
" '__rmatmul__',\n",
" '__rmod__',\n",
" '__rmul__',\n",
" '__ror__',\n",
" '__round__',\n",
" '__rpow__',\n",
" '__rsub__',\n",
" '__rtruediv__',\n",
" '__rxor__',\n",
" '__setattr__',\n",
" '__setitem__',\n",
" '__setstate__',\n",
" '__sizeof__',\n",
" '__str__',\n",
" '__sub__',\n",
" '__subclasshook__',\n",
" '__truediv__',\n",
" '__unicode__',\n",
" '__weakref__',\n",
" '__xor__',\n",
" '_accessors',\n",
" '_add_numeric_operations',\n",
" '_add_series_only_operations',\n",
" '_add_series_or_dataframe_operations',\n",
" '_agg_by_level',\n",
" '_agg_doc',\n",
" '_aggregate',\n",
" '_aggregate_multiple_funcs',\n",
" '_align_frame',\n",
" '_align_series',\n",
" '_box_col_values',\n",
" '_box_item_values',\n",
" '_builtin_table',\n",
" '_check_inplace_setting',\n",
" '_check_is_chained_assignment_possible',\n",
" '_check_label_or_level_ambiguity',\n",
" '_check_percentile',\n",
" '_check_setitem_copy',\n",
" '_clear_item_cache',\n",
" '_clip_with_one_bound',\n",
" '_clip_with_scalar',\n",
" '_combine_const',\n",
" '_combine_frame',\n",
" '_combine_match_columns',\n",
" '_combine_match_index',\n",
" '_compare_frame',\n",
" '_consolidate',\n",
" '_consolidate_inplace',\n",
" '_construct_axes_dict',\n",
" '_construct_axes_dict_for_slice',\n",
" '_construct_axes_dict_from',\n",
" '_construct_axes_from_arguments',\n",
" '_constructor',\n",
" '_constructor_expanddim',\n",
" '_constructor_sliced',\n",
" '_convert',\n",
" '_count_level',\n",
" '_create_indexer',\n",
" '_cython_table',\n",
" '_deprecations',\n",
" '_dir_additions',\n",
" '_dir_deletions',\n",
" '_drop_axis',\n",
" '_drop_labels_or_levels',\n",
" '_ensure_valid_index',\n",
" '_expand_axes',\n",
" '_find_valid_index',\n",
" '_from_arrays',\n",
" '_from_axes',\n",
" '_get_agg_axis',\n",
" '_get_axis',\n",
" '_get_axis_name',\n",
" '_get_axis_number',\n",
" '_get_axis_resolvers',\n",
" '_get_block_manager_axis',\n",
" '_get_bool_data',\n",
" '_get_cacher',\n",
" '_get_index_resolvers',\n",
" '_get_item_cache',\n",
" '_get_label_or_level_values',\n",
" '_get_numeric_data',\n",
" '_get_value',\n",
" '_get_values',\n",
" '_getitem_array',\n",
" '_getitem_column',\n",
" '_getitem_frame',\n",
" '_getitem_multilevel',\n",
" '_getitem_slice',\n",
" '_gotitem',\n",
" '_iget_item_cache',\n",
" '_indexed_same',\n",
" '_info_axis',\n",
" '_info_axis_name',\n",
" '_info_axis_number',\n",
" '_info_repr',\n",
" '_init_dict',\n",
" '_init_mgr',\n",
" '_init_ndarray',\n",
" '_internal_names',\n",
" '_internal_names_set',\n",
" '_is_builtin_func',\n",
" '_is_cached',\n",
" '_is_copy',\n",
" '_is_cython_func',\n",
" '_is_datelike_mixed_type',\n",
" '_is_label_or_level_reference',\n",
" '_is_label_reference',\n",
" '_is_level_reference',\n",
" '_is_mixed_type',\n",
" '_is_numeric_mixed_type',\n",
" '_is_view',\n",
" '_ix',\n",
" '_ixs',\n",
" '_join_compat',\n",
" '_maybe_cache_changed',\n",
" '_maybe_update_cacher',\n",
" '_metadata',\n",
" '_needs_reindex_multi',\n",
" '_obj_with_exclusions',\n",
" '_protect_consolidate',\n",
" '_reduce',\n",
" '_reindex_axes',\n",
" '_reindex_axis',\n",
" '_reindex_columns',\n",
" '_reindex_index',\n",
" '_reindex_multi',\n",
" '_reindex_with_indexers',\n",
" '_repr_data_resource_',\n",
" '_repr_fits_horizontal_',\n",
" '_repr_fits_vertical_',\n",
" '_repr_html_',\n",
" '_repr_latex_',\n",
" '_reset_cache',\n",
" '_reset_cacher',\n",
" '_sanitize_column',\n",
" '_selected_obj',\n",
" '_selection',\n",
" '_selection_list',\n",
" '_selection_name',\n",
" '_series',\n",
" '_set_as_cached',\n",
" '_set_axis',\n",
" '_set_axis_name',\n",
" '_set_is_copy',\n",
" '_set_item',\n",
" '_set_value',\n",
" '_setitem_array',\n",
" '_setitem_frame',\n",
" '_setitem_slice',\n",
" '_setup_axes',\n",
" '_shallow_copy',\n",
" '_slice',\n",
" '_stat_axis',\n",
" '_stat_axis_name',\n",
" '_stat_axis_number',\n",
" '_take',\n",
" '_to_dict_of_blocks',\n",
" '_try_aggregate_string_function',\n",
" '_typ',\n",
" '_unpickle_frame_compat',\n",
" '_unpickle_matrix_compat',\n",
" '_update_inplace',\n",
" '_validate_dtype',\n",
" '_values',\n",
" '_where',\n",
" '_xs',\n",
" 'abs',\n",
" 'add',\n",
" 'add_prefix',\n",
" 'add_suffix',\n",
" 'agg',\n",
" 'aggregate',\n",
" 'align',\n",
" 'all',\n",
" 'any',\n",
" 'append',\n",
" 'apply',\n",
" 'applymap',\n",
" 'as_matrix',\n",
" 'asfreq',\n",
" 'asof',\n",
" 'assign',\n",
" 'astype',\n",
" 'at',\n",
" 'at_time',\n",
" 'axes',\n",
" 'between_time',\n",
" 'bfill',\n",
" 'bool',\n",
" 'boxplot',\n",
" 'clip',\n",
" 'clip_lower',\n",
" 'clip_upper',\n",
" 'columns',\n",
" 'combine',\n",
" 'combine_first',\n",
" 'compound',\n",
" 'copy',\n",
" 'corr',\n",
" 'corrwith',\n",
" 'count',\n",
" 'cov',\n",
" 'cummax',\n",
" 'cummin',\n",
" 'cumprod',\n",
" 'cumsum',\n",
" 'describe',\n",
" 'diff',\n",
" 'div',\n",
" 'divide',\n",
" 'dot',\n",
" 'drop',\n",
" 'drop_duplicates',\n",
" 'dropna',\n",
" 'dtypes',\n",
" 'duplicated',\n",
" 'empty',\n",
" 'eq',\n",
" 'equals',\n",
" 'eval',\n",
" 'ewm',\n",
" 'expanding',\n",
" 'fcst',\n",
" 'ffill',\n",
" 'fillna',\n",
" 'filter',\n",
" 'first',\n",
" 'first_valid_index',\n",
" 'floordiv',\n",
" 'from_dict',\n",
" 'from_records',\n",
" 'ftypes',\n",
" 'ge',\n",
" 'get',\n",
" 'get_dtype_counts',\n",
" 'get_ftype_counts',\n",
" 'get_values',\n",
" 'groupby',\n",
" 'gt',\n",
" 'head',\n",
" 'hist',\n",
" 'iat',\n",
" 'id',\n",
" 'idxmax',\n",
" 'idxmin',\n",
" 'iloc',\n",
" 'index',\n",
" 'infer_objects',\n",
" 'info',\n",
" 'insert',\n",
" 'interpolate',\n",
" 'isin',\n",
" 'isna',\n",
" 'isnull',\n",
" 'items',\n",
" 'iteritems',\n",
" 'iterrows',\n",
" 'itertuples',\n",
" 'ix',\n",
" 'join',\n",
" 'keys',\n",
" 'kurt',\n",
" 'kurtosis',\n",
" 'last',\n",
" 'last_valid_index',\n",
" 'lat',\n",
" 'le',\n",
" 'loc',\n",
" 'lon',\n",
" 'lookup',\n",
" 'lt',\n",
" 'mad',\n",
" 'mask',\n",
" 'max',\n",
" 'mean',\n",
" 'median',\n",
" 'melt',\n",
" 'memory_usage',\n",
" 'merge',\n",
" 'min',\n",
" 'mod',\n",
" 'mode',\n",
" 'mul',\n",
" 'multiply',\n",
" 'ndim',\n",
" 'ne',\n",
" 'nlargest',\n",
" 'notna',\n",
" 'notnull',\n",
" 'nsmallest',\n",
" 'nunique',\n",
" 'pct_change',\n",
" 'pipe',\n",
" 'pivot',\n",
" 'pivot_table',\n",
" 'plot',\n",
" 'pop',\n",
" 'pow',\n",
" 'prod',\n",
" 'product',\n",
" 'quantile',\n",
" 'query',\n",
" 'radd',\n",
" 'rank',\n",
" 'rdiv',\n",
" 'reindex',\n",
" 'reindex_axis',\n",
" 'reindex_like',\n",
" 'rename',\n",
" 'rename_axis',\n",
" 'reorder_levels',\n",
" 'replace',\n",
" 'resample',\n",
" 'reset_index',\n",
" 'rfloordiv',\n",
" 'rmod',\n",
" 'rmul',\n",
" 'rolling',\n",
" 'round',\n",
" 'rpow',\n",
" 'rsub',\n",
" 'rtruediv',\n",
" 'sample',\n",
" 'select',\n",
" 'select_dtypes',\n",
" 'sem',\n",
" 'set_axis',\n",
" 'set_index',\n",
" 'shape',\n",
" 'shift',\n",
" 'size',\n",
" 'skew',\n",
" 'slice_shift',\n",
" 'sort_index',\n",
" 'sort_values',\n",
" 'squeeze',\n",
" 'stack',\n",
" 'std',\n",
" 'style',\n",
" 'sub',\n",
" 'subtract',\n",
" 'sum',\n",
" 'swapaxes',\n",
" 'swaplevel',\n",
" 'tail',\n",
" 'take',\n",
" 'to_clipboard',\n",
" 'to_csv',\n",
" 'to_dense',\n",
" 'to_dict',\n",
" 'to_excel',\n",
" 'to_feather',\n",
" 'to_gbq',\n",
" 'to_hdf',\n",
" 'to_html',\n",
" 'to_json',\n",
" 'to_latex',\n",
" 'to_msgpack',\n",
" 'to_panel',\n",
" 'to_parquet',\n",
" 'to_period',\n",
" 'to_pickle',\n",
" 'to_records',\n",
" 'to_sparse',\n",
" 'to_sql',\n",
" 'to_stata',\n",
" 'to_string',\n",
" 'to_timestamp',\n",
" 'to_xarray',\n",
" 'transform',\n",
" 'transpose',\n",
" 'truediv',\n",
" 'truncate',\n",
" 'tshift',\n",
" 'tz_convert',\n",
" 'tz_localize',\n",
" 'u10_b',\n",
" 'u10_m',\n",
" 'udir_b',\n",
" 'udir_m',\n",
" 'unstack',\n",
" 'update',\n",
" 'values',\n",
" 'var',\n",
" 'where',\n",
" 'xs']"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dir(both)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Hs_b, Hs_m, T, Tp_b, Tp_m, abs, add, add_prefix, add_suffix, agg, aggregate, align, all, any, append, apply, applymap, as_matrix, asfreq, asof, assign, astype, at, at_time, axes, between_time, bfill, bool, boxplot, clip, clip_lower, clip_upper, columns, combine, combine_first, compound, copy, corr, corrwith, count, cov, cummax, cummin, cumprod, cumsum, describe, diff, div, divide, dot, drop, drop_duplicates, dropna, dtypes, duplicated, empty, eq, equals, eval, ewm, expanding, fcst, ffill, fillna, filter, first, first_valid_index, floordiv, from_dict, from_records, ftypes, ge, get, get_dtype_counts, get_ftype_counts, get_values, groupby, gt, head, hist, iat, id, idxmax, idxmin, iloc, index, infer_objects, info, insert, interpolate, isin, isna, isnull, items, iteritems, iterrows, itertuples, ix, join, keys, kurt, kurtosis, last, last_valid_index, lat, le, loc, lon, lookup, lt, mad, mask, max, mean, median, melt, memory_usage, merge, min, mod, mode, mul, multiply, ndim, ne, nlargest, notna, notnull, nsmallest, nunique, pct_change, pipe, pivot, pivot_table, plot, pop, pow, prod, product, quantile, query, radd, rank, rdiv, reindex, reindex_axis, reindex_like, rename, rename_axis, reorder_levels, replace, resample, reset_index, rfloordiv, rmod, rmul, rolling, round, rpow, rsub, rtruediv, sample, select, select_dtypes, sem, set_axis, set_index, shape, shift, size, skew, slice_shift, sort_index, sort_values, squeeze, stack, std, style, sub, subtract, sum, swapaxes, swaplevel, tail, take, to_clipboard, to_csv, to_dense, to_dict, to_excel, to_feather, to_gbq, to_hdf, to_html, to_json, to_latex, to_msgpack, to_panel, to_parquet, to_period, to_pickle, to_records, to_sparse, to_sql, to_stata, to_string, to_timestamp, to_xarray, transform, transpose, truediv, truncate, tshift, tz_convert, tz_localize, u10_b, u10_m, udir_b, udir_m, unstack, update, values, var, where, xs\n"
]
}
],
"source": [
"methods=[x for x in dir(both) if not x.startswith('_')]\n",
"print(', '.join(methods))"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" u10_m | \n",
" udir_m | \n",
" Hs_m | \n",
" Tp_m | \n",
" fcst | \n",
" lat | \n",
" lon | \n",
" udir_b | \n",
" u10_b | \n",
" Hs_b | \n",
" Tp_b | \n",
"
\n",
" \n",
" \n",
" \n",
" | count | \n",
" 116.000000 | \n",
" 116.000000 | \n",
" 116.000000 | \n",
" 116.000000 | \n",
" 116.0 | \n",
" 116.0 | \n",
" 116.000000 | \n",
" 116.000000 | \n",
" 116.000000 | \n",
" 114.000000 | \n",
" 114.000000 | \n",
"
\n",
" \n",
" | mean | \n",
" 4.934483 | \n",
" 210.379310 | \n",
" 0.546552 | \n",
" 7.294828 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 209.034483 | \n",
" 4.875862 | \n",
" 0.562018 | \n",
" 5.947368 | \n",
"
\n",
" \n",
" | std | \n",
" 1.573924 | \n",
" 59.233166 | \n",
" 0.137955 | \n",
" 1.670555 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 0.000000 | \n",
" 56.027777 | \n",
" 2.001505 | \n",
" 0.144051 | \n",
" 2.663560 | \n",
"
\n",
" \n",
" | min | \n",
" 1.100000 | \n",
" 15.000000 | \n",
" 0.300000 | \n",
" 3.100000 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 15.000000 | \n",
" 0.100000 | \n",
" 0.310000 | \n",
" 2.000000 | \n",
"
\n",
" \n",
" | 25% | \n",
" 3.700000 | \n",
" 170.750000 | \n",
" 0.500000 | \n",
" 7.075000 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 181.000000 | \n",
" 3.500000 | \n",
" 0.462500 | \n",
" 4.000000 | \n",
"
\n",
" \n",
" | 50% | \n",
" 5.100000 | \n",
" 210.000000 | \n",
" 0.600000 | \n",
" 8.000000 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 211.000000 | \n",
" 4.950000 | \n",
" 0.535000 | \n",
" 5.500000 | \n",
"
\n",
" \n",
" | 75% | \n",
" 6.125000 | \n",
" 260.250000 | \n",
" 0.600000 | \n",
" 8.400000 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 244.500000 | \n",
" 6.200000 | \n",
" 0.637500 | \n",
" 8.000000 | \n",
"
\n",
" \n",
" | max | \n",
" 8.700000 | \n",
" 336.000000 | \n",
" 0.900000 | \n",
" 8.900000 | \n",
" 0.0 | \n",
" 31.4 | \n",
" -80.867996 | \n",
" 356.000000 | \n",
" 10.000000 | \n",
" 0.960000 | \n",
" 14.000000 | \n",
"
\n",
" \n",
"
\n",
"
Help is available for the objects bound methods
"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Help on method corr in module pandas.core.frame:\n",
"\n",
"corr(method='pearson', min_periods=1) method of pandas.core.frame.DataFrame instance\n",
" Compute pairwise correlation of columns, excluding NA/null values\n",
" \n",
" Parameters\n",
" ----------\n",
" method : {'pearson', 'kendall', 'spearman'}\n",
" * pearson : standard correlation coefficient\n",
" * kendall : Kendall Tau correlation coefficient\n",
" * spearman : Spearman rank correlation\n",
" min_periods : int, optional\n",
" Minimum number of observations required per pair of columns\n",
" to have a valid result. Currently only available for pearson\n",
" and spearman correlation\n",
" \n",
" Returns\n",
" -------\n",
" y : DataFrame\n",
"\n"
]
}
],
"source": [
"help(both.corr)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Suppose we want to look at the correlation coefficient
\n",
"for one forecast grouped by day
"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"corr=both['Hs_m'][both.fcst==0].groupby(pandas.Grouper(freq='D')).corr(both['Hs_b'])"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"datetime\n",
"2018-06-01 0.799368\n",
"2018-06-02 0.617802\n",
"2018-06-03 0.985839\n",
"2018-06-04 0.846649\n",
"2018-06-05 0.310881\n",
"Freq: D, Name: Hs_m, dtype: float64"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"corr.head()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Another way, is to select a parameter and forecast to validate
"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"model=both[both.fcst==0]['Hs_m'].copy()\n",
"obs=both[both.fcst==0]['Hs_b'].copy()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"To aggregate the data by day, we need to \"group\" the data
(think of looping through and collecting your data into chunks)
"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"obsgrp=obs.groupby(pandas.Grouper(freq='D'))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Calculate some basic stats, grouped by day
"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[(Timestamp('2018-06-01 00:00:00', freq='D'), 3),\n",
" (Timestamp('2018-06-02 00:00:00', freq='D'), 7),\n",
" (Timestamp('2018-06-03 00:00:00', freq='D'), 11),\n",
" (Timestamp('2018-06-04 00:00:00', freq='D'), 15),\n",
" (Timestamp('2018-06-05 00:00:00', freq='D'), 19),\n",
" (Timestamp('2018-06-06 00:00:00', freq='D'), 23),\n",
" (Timestamp('2018-06-07 00:00:00', freq='D'), 27),\n",
" (Timestamp('2018-06-08 00:00:00', freq='D'), 31),\n",
" (Timestamp('2018-06-09 00:00:00', freq='D'), 35),\n",
" (Timestamp('2018-06-10 00:00:00', freq='D'), 39)]"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff=model-obs\n",
"diffgroup=diff.groupby(pandas.Grouper(freq='D'))\n",
"[v for v in diffgroup.groups.items()][:10]"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"diff2group=(diff**2).groupby(pandas.Grouper(freq='D'))\n",
"count=diffgroup.count()"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"datetime\n",
"2018-06-01 0.0400\n",
"2018-06-02 0.0925\n",
"2018-06-03 -0.1425\n",
"2018-06-04 -0.0300\n",
"2018-06-05 -0.0675\n",
"Freq: D, dtype: float64"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bias=diffgroup.mean()\n",
"bias.head()"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"datetime\n",
"2018-06-01 0.094163\n",
"2018-06-02 0.115866\n",
"2018-06-03 0.144827\n",
"2018-06-04 0.050498\n",
"2018-06-05 0.092060\n",
"Freq: D, dtype: float64"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rmse=diff2group.mean()**0.5\n",
"rmse.head()"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"datetime\n",
"2018-06-01 13.350430\n",
"2018-06-02 21.144840\n",
"2018-06-03 16.770537\n",
"2018-06-04 9.358023\n",
"2018-06-05 20.958924\n",
"Freq: D, dtype: float64"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter_index=100.*(rmse - bias**2)/obsgrp.mean()\n",
"scatter_index.head()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Notice these are Series
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"It would be nice to have a DataFrame of the statistics
"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"fcst=0\n",
"agg_stats=pandas.DataFrame({'Forecast':fcst,\n",
" 'Bias':bias,\n",
" 'RMSE':rmse,\n",
" 'Corr':corr,\n",
" 'Scatter_Index':scatter_index,\n",
" 'Count':count})"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Forecast | \n",
" Bias | \n",
" RMSE | \n",
" Corr | \n",
" Scatter_Index | \n",
" Count | \n",
"
\n",
" \n",
" | datetime | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 2018-06-01 | \n",
" 0 | \n",
" 0.0400 | \n",
" 0.094163 | \n",
" 0.799368 | \n",
" 13.350430 | \n",
" 3 | \n",
"
\n",
" \n",
" | 2018-06-02 | \n",
" 0 | \n",
" 0.0925 | \n",
" 0.115866 | \n",
" 0.617802 | \n",
" 21.144840 | \n",
" 4 | \n",
"
\n",
" \n",
" | 2018-06-03 | \n",
" 0 | \n",
" -0.1425 | \n",
" 0.144827 | \n",
" 0.985839 | \n",
" 16.770537 | \n",
" 4 | \n",
"
\n",
" \n",
" | 2018-06-04 | \n",
" 0 | \n",
" -0.0300 | \n",
" 0.050498 | \n",
" 0.846649 | \n",
" 9.358023 | \n",
" 4 | \n",
"
\n",
" \n",
" | 2018-06-05 | \n",
" 0 | \n",
" -0.0675 | \n",
" 0.092060 | \n",
" 0.310881 | \n",
" 20.958924 | \n",
" 4 | \n",
"
\n",
" \n",
"
\n",
"
Looks like something that might be nice to keep around?
"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"import sqlite3"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"dbfile='my_stats_table.db'\n",
"conn = sqlite3.connect(dbfile,detect_types=sqlite3.PARSE_DECLTYPES)\n",
"agg_stats.to_sql('STATS',conn,if_exists='append')\n",
"conn.close()\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"To read it back in later
"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"conn = sqlite3.connect(dbfile,detect_types=sqlite3.PARSE_DECLTYPES)\n",
"mystats=pandas.read_sql('select * from STATS',conn)\n",
"conn.close()"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" datetime | \n",
" Forecast | \n",
" Bias | \n",
" RMSE | \n",
" Corr | \n",
" Scatter_Index | \n",
" Count | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 2018-06-01 | \n",
" 0 | \n",
" 0.0400 | \n",
" 0.094163 | \n",
" 0.799368 | \n",
" 13.350430 | \n",
" 3 | \n",
"
\n",
" \n",
" | 1 | \n",
" 2018-06-02 | \n",
" 0 | \n",
" 0.0925 | \n",
" 0.115866 | \n",
" 0.617802 | \n",
" 21.144840 | \n",
" 4 | \n",
"
\n",
" \n",
" | 2 | \n",
" 2018-06-03 | \n",
" 0 | \n",
" -0.1425 | \n",
" 0.144827 | \n",
" 0.985839 | \n",
" 16.770537 | \n",
" 4 | \n",
"
\n",
" \n",
" | 3 | \n",
" 2018-06-04 | \n",
" 0 | \n",
" -0.0300 | \n",
" 0.050498 | \n",
" 0.846649 | \n",
" 9.358023 | \n",
" 4 | \n",
"
\n",
" \n",
" | 4 | \n",
" 2018-06-05 | \n",
" 0 | \n",
" -0.0675 | \n",
" 0.092060 | \n",
" 0.310881 | \n",
" 20.958924 | \n",
" 4 | \n",
"
\n",
" \n",
"
\n",
"
Thank You!
\n",
"Any Questions?
\n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}