{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "2HN7Sb7JQ3bF" }, "source": [ "Machine learning is being used on many kinds of data, being especially successful in areas with lots of data like vision. However, we see quite little of it in meteorology, where big Numerical Weather Prediction simulations run daily that try to simulate the entire atmosphere up to weeks ahead. This takes supercomputers and lots of tweaking of the parameters to achieve good results.\n", "\n", "Funny enough we know from [The Bitter Lesson](http://www.incompleteideas.net/IncIdeas/BitterLesson.html) that data + compute always wins. Currently we take the approach of putting more and more human knowledge (physics, approximations) into these NWP models, and by and large this is [a big success story](https://www.nature.com/articles/nature14956).\n", "\n", "The question is, however, not if but *when* we will hit the limit where data + compute + more general ML methods will beat these amazing numerical simulations.\n", "\n", "In this post we will do some exploratory analysis of the NOAA \"Global Surface summary Of the Day\" dataset, which is available as [BigQuery public dataset](https://console.cloud.google.com/marketplace/details/noaa-public/gsod).\n", "\n", "> You can run this notebook in [Colab](https://colab.research.google.com)! Simply download [here](https://gitlab.com/The-Fonz/the-fonz.gitlab.io/raw/master/posts/1_Exploring_global_weather_data.ipynb?inline=false) and upload to Colab.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "yuERyp9QAbUl" }, "source": [ "## Basic imports\n", "\n", "We set up plotly and gcloud libraries." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "oZrmUKJCOHiN" }, "outputs": [], "source": [ "PROJECT = '{YOUR GCLOUD PROJECT}'" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "C21uqjr_WGTh" }, "outputs": [], "source": [ "# We want plotly 4.0 as it has much better render support\n", "!pip install --upgrade plotly>=4.0.0" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "gpSalGCnQk3y" }, "outputs": [], "source": [ "from datetime import datetime\n", "\n", "import plotly.graph_objs as go\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from scipy.spatial import Delaunay\n", "\n", "from google.colab import auth\n", "auth.authenticate_user()\n", "\n", "from google.cloud.bigquery import magics\n", "\n", "magics.context.project = PROJECT\n", "\n", "from google.cloud import bigquery\n", "\n", "bqc = bigquery.Client(project=PROJECT)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "2bx26XDCAnWX" }, "source": [ "## 1. Retrieve and plot station position\n", "\n", "Here we retrieve the station positions, triangulate it (for nice plotting of mesh), then plot them on a globe. We only use stations that are active during the nine years from 2010 to 2019." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "QgFdRqF7AZJV" }, "outputs": [], "source": [ "# CONSTANTS\n", "# Must be str\n", "ds_begin = '20100101'\n", "ds_end = '20190101'" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 142 }, "colab_type": "code", "id": "-qEde4qVQk3K", "outputId": "011ac78f-8bc7-4257-de9c-db5fc35761a9" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lonlatelevbeginend
128615.68158.40652.41995101220190725
316632.73335.16720.02009070120190725
947152.500-32.1835.02001091520190724
\n", "
" ], "text/plain": [ " lon lat elev begin end\n", "1286 15.681 58.406 52.4 19951012 20190725\n", "3166 32.733 35.167 20.0 20090701 20190725\n", "947 152.500 -32.183 5.0 20010915 20190724" ] }, "execution_count": 5, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "# Use pandas directly here, could also use bigquery cell magic\n", "df = pd.io.gbq.read_gbq(f'''\n", "SELECT lon, lat, elev, begin, d.end\n", "FROM `bigquery-public-data.noaa_gsod.stations` as d\n", "WHERE begin < '{ds_begin}' AND d.end > '{ds_end}'\n", "''', project_id=PROJECT, dialect='standard')\n", "# Convert elevation to numeric\n", "df.elev = pd.to_numeric(df.elev)\n", "df.sample(3)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "OK8PlOMdSHI6", "outputId": "a725ca56-d419-47e1-e115-539616943886" }, "outputs": [ { "data": { "text/plain": [ "(23569, 3)" ] }, "execution_count": 6, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "pts = df[['lon', 'lat']].values\n", "\n", "lat = np.deg2rad(df['lat'].values)\n", "x = np.cos(np.deg2rad(df['lon'].values)) * np.cos(lat)\n", "y = np.sin(np.deg2rad(df['lon'].values)) * np.cos(lat)\n", "z = np.sin(lat)\n", "xyz = np.stack((x, y, z)).T\n", "\n", "# Do triangulation in spherical coordinates\n", "tri = Delaunay(pts)\n", "\n", "tri.simplices.shape" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "colab_type": "code", "id": "FU9hSk-J4t_e", "outputId": "291d6712-ca4c-4a7f-aace-9f59e646c06d" }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", " \n", " \n", " \n", "
\n", " \n", "
\n", "\n", "" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "layout = dict(\n", " title = f'Weather stations active during period {ds_begin} - {ds_end}
(Click and drag to rotate)',\n", " showlegend = False,\n", " width=1000,\n", " height=1000,\n", " geo = dict(\n", " showland = True,\n", " showlakes = True,\n", " showcountries = True,\n", " showocean = True,\n", " countrywidth = 0.5,\n", " landcolor = 'rgb(230, 145, 56)',\n", " lakecolor = 'rgb(0, 255, 255)',\n", " oceancolor = 'rgb(0, 255, 255)',\n", " projection = dict( \n", " type = 'orthographic',\n", " rotation = dict(\n", " lon = -100,\n", " lat = 40,\n", " roll = 0\n", " ) \n", " ),\n", " lonaxis = dict( \n", " showgrid = True,\n", " gridcolor = 'rgb(102, 102, 102)',\n", " gridwidth = 0.5\n", " ),\n", " lataxis = dict( \n", " showgrid = True,\n", " gridcolor = 'rgb(102, 102, 102)',\n", " gridwidth = 0.5\n", " )\n", " )\n", ")\n", " \n", "go.Figure(data=[\n", " go.Scattergeo(lon=pts[:,0], lat=pts[:,1], mode='markers')\n", "], layout=layout)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "UiIjUktwTRS0" }, "source": [ "## Station positions as graph\n", "\n", "So we have all these points on the globe, but how can we feed this to an ML algorithm? The grid is not regular, not even on the sphere. One idea can be to encode these positions in a graph, enabling us to connect the stations to one another." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "colab_type": "code", "id": "E6ABS9uTxvWd", "outputId": "46802e1c-875b-4c10-a651-98fefb87e152" }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", " \n", " \n", " \n", "
\n", " \n", "
\n", "\n", "" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "def plotly_trisurf(xyz, simplices, plot_edges=None):\n", " #x, y, z are lists of coordinates of the triangle vertices \n", " #simplices are the simplices that define the triangularization;\n", " #simplices is a numpy array of shape (no_triangles, 3)\n", " #insert here the type check for input data\n", "\n", " lines = []\n", "\n", " for pt1, pt2, pt3 in tri.simplices:\n", " lines.extend([xyz[pt1].tolist(),\n", " xyz[pt2].tolist(),\n", " xyz[pt3].tolist(),\n", " xyz[pt1].tolist(),\n", " [None, None, None]\n", " ])\n", " \n", " lines = np.array(lines)\n", "\n", " return go.Scatter3d(x=lines[:,0].tolist(),\n", " y=lines[:,1].tolist(),\n", " z=lines[:,2].tolist(),\n", " mode='lines',\n", "# line=dict(color= 'rgb(50,50,50)', width=1.5)\n", " )\n", "\n", "data1 = plotly_trisurf(xyz, tri.simplices)\n", "\n", "layout = go.Layout(\n", " width=1000, height=1000,\n", " title=f'Weather stations plot as mesh'\n", ")\n", "\n", "go.Figure(data=[data1], layout=layout)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "4UY0mSry9xhf" }, "source": [ "## 2. Establishing baselines\n", "\n", "When we start training crazy triple loopback LSTM-distilled manifold-free machine learning models on this data later on, we need to know if the performance we are getting is good or not. If we predict temperature for tomorrow, is a MAE of 1.5*C good, or not?\n", "\n", "To get a better \"feel\" for these numbers, we can use some (very) simple models and see how well they perform." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "3pIrfkvD_I3W" }, "source": [ "### Model 1: Persistence\n", "\n", "One of the simplest models is just taking yesterday's value, which is usually surprisingly good as the weather tomorrow has high correlation with the weather today. Let's see what MAE this gives us." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "eIsZsTZ8_edR" }, "outputs": [], "source": [ "%%bigquery df --params {\"ds_begin\": \"20100101\", \"ds_end\": \"20190101\"}\n", "-- Select station identifier and location, and the average absolute temperature difference between two consecutive days\n", "SELECT stn, lon, lat, elev, COUNT(*) as cnt,\n", " AVG(ABS(temp_c - lag_1_temp_c)) as avg_temp_diff_lag_1,\n", " AVG(ABS(temp_c - lag_2_temp_c)) as avg_temp_diff_lag_2,\n", " AVG(ABS(temp_c - lag_3_temp_c)) as avg_temp_diff_lag_3,\n", " AVG(ABS(temp_c - lag_7_temp_c)) as avg_temp_diff_lag_7,\n", " AVG(ABS(temp_c - lag_14_temp_c)) as avg_temp_diff_lag_14\n", "FROM (\n", " SELECT *,\n", " LAG(temp_c) OVER (PARTITION BY stn ORDER BY d) as lag_1_temp_c,\n", " LAG(temp_c, 2) OVER (PARTITION BY stn ORDER BY d) as lag_2_temp_c,\n", " LAG(temp_c, 3) OVER (PARTITION BY stn ORDER BY d) as lag_3_temp_c,\n", " LAG(temp_c, 7) OVER (PARTITION BY stn ORDER BY d) as lag_7_temp_c,\n", " LAG(temp_c, 14) OVER (PARTITION BY stn ORDER BY d) as lag_14_temp_c\n", " FROM (\n", " SELECT *,\n", " -- From F to C\n", " (g.temp-32)*5/9 as temp_c,\n", " -- They use a funny format where year, mo and da are strings\n", " CAST(CONCAT(g.year, '-', g.mo, '-', g.da) as DATE) as d\n", " FROM `bigquery-public-data.noaa_gsod.gsod2018` as g\n", " LEFT JOIN `bigquery-public-data.noaa_gsod.stations` as d\n", " ON g.stn = d.usaf\n", " WHERE d.begin < @ds_begin AND d.end > @ds_end\n", " )\n", ")\n", "-- Need to group by lon/lat/elev as we want to select them, but stn alone is enough really\n", "GROUP BY stn, lon, lat, elev\n", "-- Order for display purposes\n", "ORDER BY avg_temp_diff_lag_1 DESC" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 224 }, "colab_type": "code", "id": "TQq1zi9YJHmZ", "outputId": "c19489ab-3a00-4eb9-fd2d-abbf99fe484b" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stnlonlatelevcntavg_temp_diff_lag_1avg_temp_diff_lag_2avg_temp_diff_lag_3avg_temp_diff_lag_7avg_temp_diff_lag_14
042148079.46729.033+0233.0217.833333NaNNaNNaNNaN
111142013.18347.317+0647.0214.944444NaNNaNNaNNaN
242237072.36727.133+0234.0212.555556NaNNaNNaNNaN
341216152.46424.074+0015.229.277778NaNNaNNaNNaN
4999999-79.92432.775+0003.1764038.8093458.8093458.8093458.8093458.809345
\n", "
" ], "text/plain": [ " stn lon ... avg_temp_diff_lag_7 avg_temp_diff_lag_14\n", "0 421480 79.467 ... NaN NaN\n", "1 111420 13.183 ... NaN NaN\n", "2 422370 72.367 ... NaN NaN\n", "3 412161 52.464 ... NaN NaN\n", "4 999999 -79.924 ... 8.809345 8.809345\n", "\n", "[5 rows x 10 columns]" ] }, "execution_count": 10, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "biBgIfOqHd7N" }, "source": [ "We immediately see that a handful of stations stand out, with Pantnagar Airport in India having extremely high day-to-day differences. Let's take a closer look and see what's up there." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 111 }, "colab_type": "code", "id": "OEeIbSSYIRHX", "outputId": "fbe8e1aa-2635-4bf5-bae3-466bb226d9c2" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
dtemp_c
02018-08-1832.055556
12018-12-1114.222222
\n", "
" ], "text/plain": [ " d temp_c\n", "0 2018-08-18 32.055556\n", "1 2018-12-11 14.222222" ] }, "execution_count": 11, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "%%bigquery\n", "-- Select station identifier and location, and the average absolute temperature difference between two consecutive days\n", "SELECT CAST(CONCAT(g.year, '-', g.mo, '-', g.da) as DATE) as d,\n", " (g.temp-32)*5/9 as temp_c\n", "FROM `bigquery-public-data.noaa_gsod.gsod2018` as g\n", "WHERE g.stn=\"421480\"\n", "ORDER BY d" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "xL9ewutDI-i0" }, "source": [ "Aha! It has just two readings in 2018. Let's see how many readings the other stations have. I've secretly included that in the previous query, so we can graph it right away.\n", "\n", "We exclude station \"999999\" as it seems to be some kind of \"unknown station\" code that has lots of readings. (Check `df.loc[df.cnt.idxmax()]` to find out)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 281 }, "colab_type": "code", "id": "1WjDnnMPJdgp", "outputId": "d01ba65b-fbe8-46b1-ea13-ae47501a863a" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFtVJREFUeJzt3X+U5XV93/HnSxBFVtcf2I0u6EKW\nGjeSaJ1GiZrO2iQuKrYnxxi2WCUiG3PE2JYkhepRTE3UeIxBxdhtRFK1rMbEwgItonExRht1owYo\nJaJZhCWyCrq6VI2L7/5xvwPXOfPjzsyduXM/83ycc8/c7+/P+zvfec/nvr/f+/2mqpAktet+o26A\nJGl5meglqXEmeklqnIlekhpnopekxpnoJalxJnoNTZJLkrx+RNtOkvck+WaSz4yiDS1L8j+TvHjU\n7dDimOgblmRfkgNJjukb99Ike0bYrOXydOAXgOOq6mdG3ZjVJMmeJC9dwPwXJHlf/7iqOrWq/mT4\nrdNKMNG37wjglaNuxEIlOWKBizwW2FdVdw+4/iMX3qrVYZzbrtEw0bfvzcBvJnno9AlJNiWp/sTR\n3/tLcmaSv0ry1iTfSvKVJD/bjb+1+7Qw/eP8sUmuSfKdJNcmeWzfun+im3ZXkpuSvKBv2iVJ/ijJ\nVUnuBrbO0N5HJ7m8W/7mJGd3488C/hg4JcmhJK+bYdn+WO4ELkhyvySvTnJLF8t/S7K+m/9Pkpzb\nvd/Y7aeXd8M/3rXhfkmOTXJFt3/uSvKXSWb8u+rW8RvdfvxGkjf3z5vkJUlu7MpPV0/bd5Xk5Um+\nBHxphnU/MMn7ktzZteWzSTYk+V3gGcA7un3zjm7+C7vf4beT7E3yjG78NuA/Ab/Szf/FGY6Lufbb\n1DH14iRf7eJ81Uz7Qyuoqnw1+gL2AT8P/Dnw+m7cS4E93ftNQAFH9i2zB3hp9/5M4DDwq/Q+Gbwe\n+CpwEfAA4BeB7wDruvkv6YZ/rpt+IfDJbtoxwK3duo4EngR8A9jSt+xB4Gn0OiAPnCGeTwDvBB4I\nPBH4OvDMvrZ+co59MRXLK7rtHw28BLgZOBFY1+2n93bzvwTY3b3/N8CXgQ/0Tbuse/8G4F3A/bvX\nM4DM0oYCPg48HHgM8Hd9+/pfdW15fNe+VwOfmrbsNd2yR8+w7l8DdgMP6n5XTwYeMv132jf/C4FH\ndNs6F/ja1D4HLgDeN23+/uNirv22qWvrf+328U8D3wceP+q/h7X8ske/NrwGeEWSRy5i2b+vqvdU\n1T3AB4Djgd+pqu9X1UeAfwQ2981/ZVV9oqq+D7yKXi/7eOC59Eor76mqw1X1eeDPgF/uW/ayqvqr\nqvphVX2vvxHdOp4G/Meq+l5VfYFeL/5FC4jl9qp6e7f97wJnAH9QVV+pqkPA+cDp3Seca4Gndz3u\nnwN+v9s+wL/opgP8AHgU8Niq+kFV/WV1GW8Wb6qqu6rqq8AfAtu78S8D3lBVN1bVYeD3gCf29+q7\n6Xd1bZ/uB/QS9+aquqeq9lbVt2drRFW9r6ru7PbFW+j9Y37cHO3uN9d+m/K6qvpuVX0R+CK9hK8R\nMdGvAVV1PXAFcN4iFr+j7/13u/VNH7eub/jWvu0eAu4CHk2vhv6UrqzwrSTfopcwfmymZWfwaOCu\nqvpO37hbgI0LiGX6+h/draN/fUcCG6rqy8Dd9D45PIPe/rs9yeP40UT/Znq92490JZn59nF/G27p\n2gC9/XNh3765C8i0+ObaP+8FrgZ2Jbk9ye8nuf9sMyf5za5MdLDb3nrg2HnaPmXW/dY37mt97/8f\nP3qMaIWZ6NeO1wJn86OJY+rE5YP6xvUn3sU4fupNknX0Sg2300tS11bVQ/te66rq1/uWnasnfDvw\n8CQP7hv3GGD/Ato2ff2300uw/es7zH3/3K4Fng8cVVX7u+EXAw8DvgBQVd+pqnOr6kTgecB/SPIv\n52jD8X3vH9O1AXr759em7Z+jq+pTc7T/vgm9TxOvq6otwM/S+wT1opmW6+rxvw28AHhYVT2UXtks\n822nM99+0ypjol8jqupmeqWX3+gb93V6ifKFSY5I8hLgx5e4qWcneXqSo4D/DPzvqrqVXo/4nyb5\nt0nu373+eZLHD9j+W4FPAW/oTjz+FHAW8L65l5zTpcC/T3JC90/p9+jV4Q93068FzqF3bgB6depz\n6J0LuAcgyXOTbE4SesnyHuCHc2zzt5I8rCtFvZLe7wR6df7zk/xkt971SX55tpVMl2RrkpPTu1rp\n2/RKOVPtuINePX3Kg+kl5q8DRyZ5DfCQvul3AJtmO6nM/PtNq4yJfm35HXonRfudDfwWcCfwk/SS\n6VL8d3qfHu6id0LwhdDr+dI7eXs6vR7h14A30asND2o7vZN9twMfBl5bVR9dQlsvplfy+ATw98D3\n6J2snXItvaQ4leg/Se/Tzyf65jkJ+ChwCPg08M6q+vgc27wM2EvvE8GVwLsBqurD9PbHriTfBq4H\nTl1ALD8GfIhekr+xa/t7u2kXAs/vruZ5G70Sz/+idzL4li7u/rLQn3Y/70zyNzNsa779plUmc583\nkjQsSQo4qft0Ja0Ye/SS1DgTvSQ1ztKNJDXOHr0kNW5V3Bzp2GOPrU2bNi1q2bvvvptjjpl+IUm7\njLdtxtu2Yce7d+/eb1TVvN94H2miT3IacNrmzZv53Oc+t6h17Nmzh8nJyaG2azUz3rYZb9uGHW+S\nW+afa8Slm6raXVU71q9fP8pmSFLTrNFLUuNGmuiTnJZk58GDB0fZDElqmqUbSWqcpRtJapyJXpIa\nZ41ekhpnjV6SGrcqvhkrafXYdN6V977f98bnjLAlGhZr9JLUOGv0ktQ4a/SS1DhLN5LUOBO9JDXO\nRC9JjfNkrCQ1zpOxktQ4SzeS1DgTvSQ1zkQvSY0z0UtS40z0ktQ4L6+UpMZ5eaUkNc7SjSQ1zkQv\nSY0z0UtS40z0ktQ4E70kNc5EL0mNM9FLUuP8wpQkNc4vTElS4yzdSFLjTPSS1DgTvSQ1zkQvSY07\nctQNkLR6bTrvynvf73vjc0bYEi2FPXpJapyJXpIaZ6KXpMaZ6CWpcSZ6SWqciV6SGjf0RJ/k8Une\nleRDSX592OuXJC3MQIk+ycVJDiS5ftr4bUluSnJzkvMAqurGqnoZ8ALgacNvsiRpIQbt0V8CbOsf\nkeQI4CLgVGALsD3Jlm7a84ArgauG1lJJ0qKkqgabMdkEXFFVT+iGTwEuqKpndcPnA1TVG/qWubKq\nZvw6XZIdwA6ADRs2PHnXrl2LCuDQoUOsW7duUcuOI+Nt22qI97r9Mz8f4uSNw7+d+GqIdyUNO96t\nW7furaqJ+eZbyi0QNgK39g3fBjwlySTwS8ADmKNHX1U7gZ0AExMTNTk5uahG7Nmzh8UuO46Mt22r\nId4z+2570G/fGZND39ZqiHcljSreod/rpqr2AHsGmTfJacBpmzdvHnYzJC3AplmSu9qwlKtu9gPH\n9w0f140bmE+YkqTlt5RE/1ngpCQnJDkKOB24fDjNkiQNy6CXV14KfBp4XJLbkpxVVYeBc4CrgRuB\nD1bVDQvZuA8Hl6TlN1CNvqq2zzL+KpZwCWVV7QZ2T0xMnL3YdUiS5uYtECSpcSNN9JZuJGn5jTTR\ne9WNJC0/SzeS1DgTvSQ1zhq9JDXOGr0kNc7SjSQ1zkQvSY2zRi9JjbNGL0mNs3QjSY0z0UtS40z0\nktQ4T8ZKUuM8GStJjbN0I0mNG+gJU5Las+m8K0fdBK0Qe/SS1DgTvSQ1zkQvSY3z8kpJapyXV0pS\n4yzdSFLjTPSS1DgTvSQ1zi9MNaT/CzD73vicEbZE0mpij16SGmeil6TGmeglqXF+YUqSGjfSk7FV\ntRvYPTExcfYo2yGtFd6xcm3yqhupcSZ3WaOXpMaZ6CWpcSZ6SWqciV6SGmeil6TGmeglqXEmeklq\nnIlekhrnF6akBvklKfWzRy9JjVuWHn2Sfw08B3gI8O6q+shybEeSNL+Be/RJLk5yIMn108ZvS3JT\nkpuTnAdQVf+jqs4GXgb8ynCbLElaiIWUbi4BtvWPSHIEcBFwKrAF2J5kS98sr+6mS5JGJFU1+MzJ\nJuCKqnpCN3wKcEFVPasbPr+b9Y3d65qq+ugs69oB7ADYsGHDk3ft2rWoAA4dOsS6desWtew4mive\n6/bfd1//kzeuX6kmLSt/v4vTfywMy3IcU/5+l2br1q17q2pivvmWWqPfCNzaN3wb8BTgFcDPA+uT\nbK6qd01fsKp2AjsBJiYmanJyclEN2LNnD4tddhzNFe+Z/Q8HP2PmecbNTPG2/BD0YR3PZy7DVTfL\ncUz597syluVkbFW9DXjbfPMlOQ04bfPmzYve1nX7D957ULf2Ry9Jw7DUyyv3A8f3DR/XjRtIVe2u\nqh3r17dRZpCk1Wipif6zwElJTkhyFHA6cPnSmyVJGpaFXF55KfBp4HFJbktyVlUdBs4BrgZuBD5Y\nVTcsYJ0+HFySltnANfqq2j7L+KuAqxazcR8OLknLz1sgSFLjRnpTs2FcdSOpxxuZaTYj7dF71Y0k\nLT9LN5LUOO9HL40xyzUaxEh79F5eKUnLzxq9JDXO0o00Aku5MZvlGi2UJ2MlqXHW6CWpcSMt3XgL\nBGl2Ld93XyvLGr3mZcIZPevyWgoTvbRCTNYalTWR6O2RSlrLvKmZNGL9HZFLth0zwpaoVZ6M1Vjz\n05o0v2ZLN+NaDzVxSRq2ZhO92rIc/7hX4z/V6/Yf5Mwx7aRo9fKbsZLUOBO9JDXORC9JjfPySmkI\nVmO9X5ri5ZVaU1biaqxxveJL7Wrqqhv/wDQM9s7VmqYSfWtMOIvnvpPu48lYSWqcPfoxt9xfJNJ9\npu8XPyloXJjopUXyH6LGhYlei7YW6uAmc7XAGr0kNc4vTI2hUfYyh7Xt1fZpwJ67WjbSHn1V7a6q\nHevXrx9lMySpadbotaxWW89dWous0UtS4+zRqxnW2aWZmeg1dEtJuJZ6pOGzdCNJjbNHPyamerrn\nnnyY1n5tllyk5dVWxhihcS85tJxsW45NozcOf/uWbiSpcWu6Rz9XT2+Q/8yzLT8O/+GHbZBesz1r\naTTs0UtS49Zcj3419ipXY5sktWPoPfokJyZ5d5IPDXvdkqSFG6hHn+Ri4LnAgap6Qt/4bcCFwBHA\nH1fVG6vqK8BZ457oZ+tlL7TmPqz1DMtaPH8grXWD9ugvAbb1j0hyBHARcCqwBdieZMtQWydJWrJU\n1WAzJpuAK6Z69ElOAS6oqmd1w+cDVNUbuuEPVdXz51jfDmAHwIYNG568a9euRQVw4K6D3PHdRS06\nUidvvO/WzNftPzjwchuOZqB4+9ffr39bi23DSho03las5nhnO15mO9YGcejQIdatW7ekdo3CbH8v\n8+2LYce7devWvVU1Md98SzkZuxG4tW/4NuApSR4B/C7wpCTnTyX+6apqJ7ATYGJioiYnJxfViLe/\n/zLect34nVPed8bkve/PXMDJ2HNPPjxQvP3r79e/rcW2YSUNGm8rVnO8sx0vsx1rg9izZw+L/dsf\npdn+XubbF6OKd+hHVFXdCbxskHl9wpQkLb+lXHWzHzi+b/i4btzAfMKUJC2/pST6zwInJTkhyVHA\n6cDlw2mWJGlYBr288lJgEjg2yW3Aa6vq3UnOAa6md3nlxVV1w0I2bulmtPyilrQ2DJToq2r7LOOv\nAq5a7Marajewe2Ji4uzFrkOSNDfvdSNJjRvpdVxruXRj2UTjxmN2fvPto3NPPsyZ51254t9KH2mP\n3qtuJGn5WbqRpMaZ6CWpcdbo1wBrq1pJ436H1HFv/0ys0UtS4yzdSFLjTPSS1Dhr9I2yLq/VZhxr\n38vV5pXeF9boJalxlm4kqXEmeklqnIlekhrnyVhJS7IaT/zP1qb+E58LbfdqjHNQnoyVpMZZupGk\nxpnoJalxJnpJapyJXpIaZ6KXpMZ5eaWkZTPbJYlT4889+TCTK9ietcrLKyWpcZZuJKlxJnpJapyJ\nXpIaZ6KXpMaZ6CWpcSZ6SWqciV6SGucXpiSN1CD3jtfS+IUpSWqcpRtJapyJXpIaZ6KXpMaZ6CWp\ncSZ6SWqciV6SGmeil6TGmeglqXEmeklqnIlekhpnopekxg39pmZJjgHeCfwjsKeq3j/sbUiSBjdQ\njz7JxUkOJLl+2vhtSW5KcnOS87rRvwR8qKrOBp435PZKkhZo0NLNJcC2/hFJjgAuAk4FtgDbk2wB\njgNu7Wa7ZzjNlCQtVqpqsBmTTcAVVfWEbvgU4IKqelY3fH43623AN6vqiiS7qur0Wda3A9gBsGHD\nhifv2rVrUQEcuOsgd3x3UYuOpQ1HY7wNM96ZnbzxvluZX7f/4LLOs5xmire/TQu1devWvVU1Md98\nS6nRb+S+njv0EvxTgLcB70jyHGD3bAtX1U5gJ8DExERNTk4uqhFvf/9lvOW6kT4/ZUWde/Jh422Y\n8c5s3xmT974/c7YHlQxpnuU0U7z9bVouQz+iqupu4FcHmdcnTEnS8lvK5ZX7geP7ho/rxg3MJ0xJ\n0vJbSqL/LHBSkhOSHAWcDlw+nGZJkoZl0MsrLwU+DTwuyW1Jzqqqw8A5wNXAjcAHq+qGhWw8yWlJ\ndh48uPInRSRprRioRl9V22cZfxVw1WI3XlW7gd0TExNnL3YdkqS5eQsESWrcSBO9pRtJWn4jTfRe\ndSNJy2/gb8YuayOSrwO3LHLxY4FvDLE5q53xts142zbseB9bVY+cb6ZVkeiXIsnnBvkKcCuMt23G\n27ZRxevJWElqnIlekhrXQqLfOeoGrDDjbZvxtm0k8Y59jV6SNLcWevSSpDmY6CWpcWOb6Gd5Xu1Y\nm+nZvEkenuSaJF/qfj6sG58kb+vi/9sk/2x0LV+cJMcn+XiS/5PkhiSv7MY3GXOSByb5TJIvdvG+\nrht/QpK/7uL6QHc3WJI8oBu+uZu+aZTtX6wkRyT5fJIruuHW492X5LokX0jyuW7cSI/psUz0czyv\ndtxdwrRn8wLnAR+rqpOAj3XD0Iv9pO61A/ijFWrjMB0Gzq2qLcBTgZd3v8dWY/4+8Myq+mngicC2\nJE8F3gS8tao2A98EzurmP4veYzk3A2/t5htHr6R3h9sprccLsLWqnth3zfxoj+mqGrsXcApwdd/w\n+cD5o27XkGLbBFzfN3wT8Kju/aOAm7r3/wXYPtN84/oCLgN+YS3EDDwI+Bt6j9/8BnBkN/7eY5ve\nLcBP6d4f2c2XUbd9gXEeRy+xPRO4AkjL8XZt3wccO23cSI/psezRM/PzajeOqC3LbUNV/UP3/mvA\nhu59U/ug+5j+JOCvaTjmrozxBeAAcA3wZeBb1Xu+A/xoTPfG200/CDxiZVu8ZH8I/Dbww274EbQd\nL0ABH0myN8mObtxIj+m18xTiBlRVJWnuetgk64A/A/5dVX07yb3TWou5qu4BnpjkocCHgZ8YcZOW\nTZLnAgeqam+SyVG3ZwU9var2J/knwDVJ/m//xFEc0+Pao1/y82rHyB1JHgXQ/TzQjW9iHyS5P70k\n//6q+vNudNMxA1TVt4CP0ytdPDTJVKerP6Z74+2mrwfuXOGmLsXTgOcl2Qfsole+uZB24wWgqvZ3\nPw/Q+2f+M4z4mB7XRL+Wnld7OfDi7v2L6dWxp8a/qDtr/1TgYN9Hw7GQXtf93cCNVfUHfZOajDnJ\nI7uePEmOpnc+4kZ6Cf/53WzT453aD88H/qK6Qu44qKrzq+q4qtpE72/0L6rqDBqNFyDJMUkePPUe\n+EXgekZ9TI/6xMUSTng8G/g7ejXOV426PUOK6VLgH4Af0KvVnUWvRvkx4EvAR4GHd/OG3pVHXwau\nAyZG3f5FxPt0evXMvwW+0L2e3WrMwE8Bn+/ivR54TTf+ROAzwM3AnwIP6MY/sBu+uZt+4qhjWELs\nk8AVrcfbxfbF7nXDVG4a9THtLRAkqXHjWrqRJA3IRC9JjTPRS1LjTPSS1DgTvSQ1zkQvSY0z0UtS\n4/4/C2yZUuYPBkEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "plt.semilogy()\n", "plt.title(\"Number of rows per station\")\n", "df[df.stn!='999999'].cnt.hist(bins=100);" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "mci19-BALgOY" }, "source": [ "This is kind of strange, as we would expect most stations to have 365 readings (one per day) and some to have less, which is clearly the case, but there is also a small amount of stations that have more than 365 readings. So we just filter the stations with readings almost every day, which is the bulk of them anyway." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 80 }, "colab_type": "code", "id": "GNSRaAAfLg92", "outputId": "cc7fc172-cb08-4745-ecae-eb5880361e98" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
n
0147
\n", "
" ], "text/plain": [ " n\n", "0 147" ] }, "execution_count": 13, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "%%bigquery\n", "SELECT COUNT(*) as n\n", "FROM (\n", " SELECT stn, COUNT(*) as cnt, COUNT(DISTINCT dt) as cnt_days\n", " FROM (\n", " SELECT *, CAST(CONCAT(g.year, '-', g.mo, '-', g.da) as DATE) as dt\n", " FROM `bigquery-public-data.noaa_gsod.gsod2018` as g\n", " LEFT JOIN `bigquery-public-data.noaa_gsod.stations` as d\n", " ON g.stn = d.usaf\n", " WHERE d.begin < '20100101' AND d.end > '20190101' AND stn != '999999'\n", " )\n", " -- Need to group by lon/lat/elev as we want to select them, but stn alone is enough really\n", " GROUP BY stn\n", " HAVING cnt != cnt_days\n", ")" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "NWV0xMedMfxM" }, "source": [ "We see that there are just 147 stations with more readings than 365, so we can safely ignore them in further analysis. At some point we can revisit this and include stations where readings are not ambiguous." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "I129SeknHASV" }, "source": [ "We have now basically retrieved the MAE per station over the whole year of this very simple model, and we've done it for a few different days to get more feel of how our persistence model does over the span of two weeks.\n", "\n", "Let's graph the results." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 282 }, "colab_type": "code", "id": "kKiKfyTHNhRI", "outputId": "24ad4407-1188-4ad9-d850-6594cbc0e108" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXmcXEW1+L+n1+mePTOTmSSTZJKQ\nlUASEhIIOypb2ARFVBBlUxZZ9L2HCjz1CYoom+ATUfaH+JNNQSEERHYCRMieTPY9k8wks09Pb7d+\nf9zuzmSYrWe6+/ZS389nPtPd996qc6v71qmqc+ocUUqh0Wg0Gg2AzWoBNBqNRpM+aKWg0Wg0mhha\nKWg0Go0mhlYKGo1Go4mhlYJGo9FoYmiloNFoNJoYWiloYojIiSKyIwnlKhE5JNHlapKPiLwpIpdH\nXn9dRBZ1OXaMiKwXkTYROVdEKkXkbRFpFZG7rJNaMxS0UuiHyEPRKCJuq2XJBYaqQLJJAYnIYyJy\nm9VyRFFKPaWUOqXLR/8DPKCUKlBK/RW4EmgAipRS37dESM2Q0UqhD0SkBjgOUMDZSarDkYxyNelN\nKr73FNQxFljV7f1qNYgdsfo5SCOUUvqvlz/gv4H3gLuBv3f5fB5QB9i7fPZFYHnktQ34AbAR2Af8\nBRgWOVaDqWQuA7YBb0c+fyZSZjPwNnBol7LLgJeAFuBj4Dbg3S7HpwCvAfuBWuCCPu7pW8AaoBXY\nBHy7y7ETgR3AjzBHfFuAr3c5fgawOnLtTuA/uhy7AtgQkeFFYGSXYwo4JPL6TeDyLse+Gb2XyH0r\noB1oA74S+fxMYCnQBLwPHN7LvcV9feQe/xNYHrnuYaASeCVyn68Dpd2+uyuBXcDubm2QsO89UkcQ\nCETu5aXubRl5/xhwW7fv76ZImU/G036Rc78ArI3I8wDwVvT76vZdbQQMwBeR7+lu8n5+kO1xVETG\nJmAZcGIX2d4Efob5TLYCi4DyLseP7XLtduCbkc/dwK8j9ewBHgQ8Vvcv6fpnuQDp/IfZyV0NzI78\n4Cu7HNsIfKHL+2eAH0ReXw8sBqojP8jfA09HjkUfhieA/OiPE7gUKIycfy+wtEvZf478eYFpkR98\n9OHMj7z/FuAAZmF26NN6uacFwARAgBOADuCIyLETgRCmEnRHjrcDkyPHdwPHRV6Xdrnu5EidR0Su\nuz/6kEeOD0gpdD838n4WsBdTEduBSzA7cncv9xfX9ZHXizEVwajIuZ9ErssD3gB+3O27ezrS7ocB\n9cDnk/S9P0akw+/j/mLndPn+fhkpzxNP+wHlmJ3tlwAncGOkvM8ohS5t9/ne5I23PSLtvw9z8GHD\nVFD7gIouv52NwKTI+W8Cd0SOjY3I/tWI7GXAzMixezAHKsMibf0S8Aur+5d0/bNcgHT9wxx1BImM\nRDBHTzd2OX4b8EjkdSFm5zk28n4N8Lku546IlOXo8jCM76Puksg5xZEHOUikY+5Sd1QpfAV4p9v1\nvyfSkQ3gPv8KXB95He1U8rsc/wtwa+T1NuDbmGvGXct4GLizy/uCiMw1kfdDUQq/A37Wrb5a4IRe\n7ieu6/nsbOg54Hdd3n8X+GvkdfS7m9Ll+J3Aw4n+3iPvHyN+pRAA8gbTfsA3gMVd3gvmzGOwSiGu\n9sCc4TzZTaZXgUu6/HZu6XLsamBh5PUPgRd6uCfBfDYndPnsaGDzQJ6PXPzTNoXeuQRYpJRqiLz/\nU+Qzurw/L2KAPg/4RCm1NXJsLPCCiDSJSBPmwxHGHI1G2R59ISJ2EblDRDaKSAvmwwbmyK0C8yHa\n3tO1kbrmReuK1Pd1oKqnmxKR00VksYjsj5x7RqSeKI1KqfYu77cCIyOvz4+cv1VE3hKRoyOfj4yc\nB4BSqg1zhDeqJxniZCzw/W73N7qLTIm4fk+X174e3hd0K7Nr+3dtn0R+74OlXinV2eV9PO03sqt8\nyuxBt/dw3kCJqz0i53+5m6zHYiqTKHVdXndw4LsZjTmL6E4F5gz7313KXBj5XNMD2rjTAyLiAS4A\n7CIS/RG6gRIRmaGUWqaUWi0iW4HTga9hKoko24FLlVLv9VB2TeSl6vLx14BzMNdht2DOEBoxRzn1\nmKP3amBd5PzR3ep6Syn1hQHclxtzJPwN4G9KqaCI/DVST5RSEcnvohjGACsBlFIfA+eIiBO4FnMW\nMRpzfX1sl3ryMafvO3sQox3zIY3So/Lqdn+3K6Vu7+/+knR9T4zGnDmC2T67utSVqO+9+7lROvhs\n+3V1I+5+TTz3v5suvy0REQ7+rcVLvO2xHXOmcMUg65rbw+cNmIr9UKVUT79HTTf0TKFnzsUc0UwD\nZkb+pgLvYHaoUf6EuW56PKZNIcqDwO0iMhZARCpE5Jw+6isE/Jijay/w8+gBpVQYeB74iYh4RWRK\nNxn+DkwSkYtFxBn5O1JEpvZQjwtTudUDIRE5HTilh/N+KiIuETkO00j5TOT910WkWCkVxDR6G5Hz\nnwa+JSIzI4rn58CHSqktPZS9FHOG5Y24jl7W7fgeYHyX938AviMi88QkX0QWiEhhD2Un4vqBcGtE\n/kMxbTn/L/J5wr73Xu4FzPb7WmSWcRqm3acv4rn/fwCHish5EW+g6+hfafdFvO3xf8BZInJq5P7y\nxNw7Uz2Aup4CPi8iF4iIQ0TKRGSmUsrAbIN7RGR4RI5RInLqEO4rq9FKoWcuAR5VSm1TStVF/zC9\nMb7exX3uacyH8o0uy0wA92EathaJSCumsW1eH/U9gbkMsRPTu2dxt+PXYo4i64AnI/X6AZRSrZgd\n+4WYI9Y6DhgaDyJy7nWYI/xGzJHqi91Oq4sc24X5oH1HKRUdFV8MbIksdXwHc5kKpdTrwK2Ys5Dd\nmIbsC3u513sw1733AI9H6ujKT4DHI1P9C5RSSzA9mx6IyLUBc227N4Z6/UB4K1LOP4FfK6WiG7oS\n/b0/DEyL3MtfI59dD5yF6WHzdUybUK/Ec/+R3/CXgTswFdVETE+fwRJXeyiltmPOnH6EOXDZjukZ\n1m8/pZTahrm0+X1MD7ilwIzI4Zsw73tx5Lf7OjB5cLeU/UjE8KLJIETkl0CVUuqSfk/WJIzIksdm\nwKmUClkrjUaTHPRMIQMQkSkicnhk+j8Xc8nlBavl0mg02Yc2NGcGhZhLRiMxl13uAv5mqUQajSYr\n0ctHGo1Go4mhl480Go1GEyOly0fl5eWqpqYmlVVqNBpNxvPvf/+7QSmVkg13KVUKNTU1LFmyJJVV\najQaTcYT2SibEvTykUaj0WhiaKWg0Wg0mhhaKWg0Go0mhlYKGo1Go4nRr1IQkUdEZK+IrOzy2TAR\neU3MpN2viUhpcsXUaDQaTSoYyEzhMeC0bp/9APinUmoiZlCwHyRYLo1Go9FYwECiD76NGXWwK+dg\nRrgk8v/cBMul0Wg0GgsYrE2hUim1O/K6joMzKWk0Go0mQxmyoTmSsq/XAEoicqWILBGRJfX19UOt\nTqPRaFKCz+fjiiuuYNmyZVaLklIGqxT2iMgIgMj/vb2dqJR6SCk1Ryk1p6JCp0XVaDSZwebNm6mt\nreWhhx6yWpSUMlil8CIHkthfgg7jrNFospRgMGi1CCllIC6pTwMfAJNFZIeIXIaZru8LIrIeM+n4\nHckVU6PRaKxBRKwWIaX0GxBPKfXVXg59LsGyaDQaTdqRazln9I5mjUaj6YNcmylopaDRaDR9oGcK\nGo1Go4mhZwoajUajiaFnChqNRqOJoWcKGo1Go4mhZwoajUajiaFnChqNRqOJoWcKGo1Go4mhZwoa\njUajiaFnChqNRqPJWbRS0Gg0Gk0MrRQ0Go2mD7RNQaPRaDQxtE1Bo9FoNDFloGcKGo1Go4kpBT1T\n0Gg0Go2eKWg0Go3mAIZhAHqmoNFoNBr0TEGj0Wg0XdA2BY1Go9HEiC4f6ZmCRqPJaZ5//nn27t1r\ntRhpg54paDSanGXv3r3ce++93HfffVaLYjl6pqDRaHKeUCgEwPr16y2WxHq095FGo9FEyLXRcU9E\nlUKuoZWCRqP5DLk2Ou4JrRQ0Go0mgp4pQDgctloES9BKQaPRfAY9U9BKQaPRaDRd0EpBo9FoIujl\nI3jllVeA3Js1DUkpiMiNIrJKRFaKyNMikpcowTQajXXkWkfYE3V1dVaLYAmDVgoiMgq4DpijlJoO\n2IELEyWYRqOxDj1T0AHxBosD8IiIA/ACu4YukkajsYpc3bDVE7naBoNWCkqpncCvgW3AbqBZKbWo\n+3kicqWILBGRJfX19YOXVKPRJJ1c9c3vCa0U4kRESoFzgHHASCBfRC7qfp5S6iGl1Byl1JyKiorB\nS6rRaJKOVgoH0KGz4+fzwGalVL1SKgg8D8xPjFgajcYKom6YubaO3hO5qiCHohS2AUeJiFfMX9Dn\ngDWJEUuj0VhBrnaEPRGdIQQDQYslSS1DsSl8CDwLfAKsiJT1UILk0mg0FqANzQeItoU/ELBYktTi\nGMrFSqkfAz9OkCwajcZi9PLRAaKK0d/ZabEkqUXvaNZoNDFyNbRDT0TbolMrBY1Gk6topXCA6PKR\nr9OXU8tpWiloNJoY0cxrevnogIIMh8P4/X6LpUkdWiloNJoYUaWQSyPj3ug6a2pra7NQktSilYJG\no4kRDOaW+2VfhENhwJwxtba2WitMCtFKQaPRxIiOjvVMAULhEHab6aDZ3NxssTSpQysFjUYTI5Bj\nPvm94ff7UUrh0EpBo9HkMtHEMrluaG5qagLAbnMe9D4X0EpBo9HE2LNnj9UipAX79+8HwBFRCo2N\njVaKk1K0UtBoNDFisY9y3KYQVQI2mx230xNTErmAVgoajSZGLFx0jgfGi+Z+sYkdj7OAhoYGiyVK\nHUOKfaTRaLILHRDPJKoEbGLH5chj7569FkuUOvRMQaPRxIgqg1z3Qqqrq8MmdkQEr6uIuro6q0VK\nGVopaDSaGNGZQmcOhXXoiV27dmG32QHIdxfT0tpCe3u7xVKlBq0UNBpNjK45BHJ5CWnb1m3YxfQ8\nKswrBWD79u1WipQytFLQaDQxojuaDcPI2SWkpqYmmluasdtNk2tRXhkAW7dutVKslKGVgkajidE1\nHafP57NQEuvYuHEjAA6bC4CCvFLsNgcbNmywUqyUoZWCRoMZBTOXl0uidI0Mmitr6N1Zt24dcGDj\nmk1sFHsrWLt2rZVipQytFDQ5TzAY5IwzzuCJJ56wWhTLydVw0V1ZuXIlhZ5SbBFDM8Aw7whq19bG\nQotnM1opaHKe6Nr5008/bbEk1hMOhWKbl3IpXHQUpRTLl69gmHfEQZ+XF4yk098ZW1rKZrRS0Gg0\ngNkhhsJhnJH3LS0tlspjBVu3bqW5uYmKwtEHfV5eWA3A0qVLrRArpWiloNFEyHWbQkdHB0BMKeRS\nELgoy5YtA/iMUvC6CinMK9VKQaPR5A7RoG8uzHxjuagUli5ditddSIG75DPHyguqWbZ02UEeWtmI\nVgqanCf6kOd6DoF9+/YBYAcKbbacCgIXZdnSZZTlj+rxt1BRWE1bextbtmxJvWApRCsFTc6jU1Ca\nROP7OIAipdiTQ/F+wAyC17CvgbL8ET0eH5Y/EoDVq1enUqyUo6OkanKe3//+91aLkBbs2rULOKAU\ndu3caa1AKWb9+vUAlOZX9Xi8MK8Ul8Md28eQreiZgibn2bx5s9UipAU7duzAgWlPKAP27N2bU6Eu\nostCxZ7yHo+LCIV5ZVkf7kIrBU3Ok+vLRlE2b9oU8zyqAAylciYIHJhKMc/lxeXI6/WcAncp27dl\nd5topaDJebRSgFAoxNZt22JKoTLyPxc2a0XZs2cPXmdRn+d4XUXs278/q3c2a6WgyXmy3cVwIGzZ\nsoVQKIQ78r4ccIpk/fp5V+r31pPnLOjzHK+rAKWMrM7ZPCSlICIlIvKsiKwVkTUicnSiBNNoUoWe\nKRAL9uaKvLcjVAJr16yxTKZUs3//fvKc+X2eEz2ulULv3AcsVEpNAWYAufMLygI6Ojp0h4jepwCw\nZs0aPDZbbPkIoFopamtzIwhcIBCguaUZTz8zhejxbN7DMWilICLFwPHAwwBKqYBSqilRgmmSz4IF\nC3j44YetFsNy9D4FWL1qFaO6LaNVY2Zgy3ZvG4D6+noAvO5+bAqR43v27Em6TFYxlJnCOKAeeFRE\nPhWRP4rIZ+ZeInKliCwRkSXRhtekB+FwWIeL5sBMIVeVQkdHB5s3b6a62+fR6D8rV65MtUgpJ+pl\n1VN4i664HV6cDjc7duxIhViWMBSl4ACOAH6nlJoFtAM/6H6SUuohpdQcpdScioqKIVSn0SSH2EyB\n3FQKa9euxVCK0d0+LwXybbas38ELsGnTJqD3PQpRRISivLKs9soailLYAexQSn0Yef8sppLQZAC5\nOiruieiaea56Ia1atQrgMzMFQag2DFYsX556oVLMqlWrKMgr6XOPQpRSbyVr16zNWlvLoJWCUqoO\n2C4ikyMffQ7I/iFFlpCrHWBPRB/uYCBIZ2enxdKknmXLljHcZsPLZw3tY4EdO3fS1JS95kLDMFi2\nbBnl+aMGdH5FYTWd/s6sddcdqvfRd4GnRGQ5MBP4+dBF0qQCrRQO0HXEl82uhj0RCoVYsWwZNb38\nHmoi/z/55JOUyZRqamtraWlpYXjR2AGdX1E4BhA++uij5ApmEUNSCkqppRF7weFKqXOVUrkXgD1D\neeCBB6wWIW0IBAIxe0KuOUOsXLkSn9/PhF6OjwTybDY+/vjjVIqVUt5//31EhBHF4wZ0fp7TS3nB\nCN57970kS2YNekdzjrJhwwarRUgbAoFALF5wNrsa9sS7776LXYTxvRy3IxxiGLz/7rsxg3y28dZb\nb1NeMAq30zvga0aUHELtutqs/L1opaDJaTo7O83OLrJrqy6HcggYhsGbb/yLCUqR14M9Ico0oLG5\nmeVZaHDevn07W7ZsZlTJpLiuqy6dCMDbb7+dDLEsRSuFHEV7H5lEcwjgAJvXxs4cyiHw6aefsreh\nnhn9nDcZcIvwyiuvpEKslBLt1EdFOvmBUpg3jGJvuVYKmuxBG5pNtm3bZr6wQzg/nBO7d6O89NJL\neMTG1H7OcyEcphT/euMNWltbUyJbqnj77XcYll9Ffj87mXtiZPFEli9fkXWeWVop5Ch6pmAS3bSE\nA1SRYtPmTTmhMBsaGnjrzTeZpQycfSwdRTkSM+RFNs0WGhsbWbt2DSNKejOz982o0kNQyuDDDz/s\n/+QMQiuFHCUXOr6BsLZ2rZmpXoBS6PR1ZnUIgygvvvgihmEwd4Dnj0QYI8Jzzz6bNQbnxYsXo5Ri\nZPHglEKptxKPq4APPvggwZJZi1YKOUq2PNhDIRwOs3z5cpTTnDWpMvN/NhpUuxIIBPjrCy8wCSgb\nwCwhytFKsbuuLms6wQ8//BCPK58S7/BBXS8iVBbV8NGHH2XV86SVQo6SrVv042HNmjV0tHfEPI8o\nBPFK1i0HdOfNN9+kqbmZo+K8bipQZLPxwvPPJ0OslGIYBkuWLGF44dghhUyvKqqhrb2N2traBEpn\nLVop5ChaKcA777xjPgHRzDIC4aowiz9cnNXhLl78298os9l63ZvQG3aEOYbBx0uWHPDaylA2bNhA\nS0sLlUU1QypneNEYAP79738nQKr0QCuFHCUYDMZe+3w+CyWxhnA4zGuvv4Yarg56CtRohb/Tz3vv\nZedu1V27drF8xQpmGQa2OJaOoszCNL+8+uqrCZctlUQ78coBhrbojTxnPqX5w1myZEkixEoLtFLI\nUfx+f+z17t27LZTEGj766CMa6hswxnUzuFeA5AsvvviiNYIlmahf/eGDvL4EYSzw1ptvJkokS/j0\n008p8gzD4+o709pAKC8YzcqVqw56pjIZrRRyFL/fj01Mw2o0wUgu8cwzzyAeMYP7dEUgPC7Mp59+\nesBdNYt4/733qBKhdBCzhChTgE2bN2dsiIdwOGxGRS3onkHiAJ9ue4Omjr00dezlX2v/zKfb3uj1\n3OGFYwgGA6zJknzWWinkIOFwmM7OTgqdCgE2b95stUgppba2liVLlhCeEO7xCVDjFeIQnnrqqdQL\nl0SCwSCrV69m/BD3qERtEcuWLRu6UBawefNmfD4fFQW9h8pu6thLMOwnGPZT37qdpo69vZ5bXmCO\nLLIlQ51WCjnItm3bUErhdSiq8lXWjHAGyqOPPoq4BHVIL52jG8Ljw7z++uts2bIlpbIlk02bNhEI\nBj+TYS1eKgGXSMb+bqKZ5MoKuk8TB4fb6aXIMyxrMtRppZCDREd4HodiYlGAlSuWZ5WfdV8sXbqU\n999/n/CkA0HwekJNUeCAB3//YOqESzLR5bDKIZZjQxgObM7Q5bUNGzbgcrjJ7ycfczwUe4azYX12\nRB7WSiEHWbx4MQ4buGyK6WVBWtvaM3bUFw+hUIi77r4LyRfUpH6WUNwQnhLm/ffez5rNWtu3b8cG\nDEtAWeVKsS1D40Rt2rSJorzyIe1P6E6xp5y6PXV0dHQkrEyr0Eohx2htbeXjjz6k0Gl63cwoC+Gw\nwb/+9S+LJUs+Tz31FFu3bCU0I2SGtugHNUkhRcKvfv0r2traki9gktm5cyelNhv2IRiZowwDGvbv\nz0iPmx3bd1CQV5rQMgvzTFWb6fs3QCuFnGPRokUEQ2GKXaZSyHcqZpUHeHXhKxn5gA+U2tpaHn3s\nUYzRBgwsFS/YIDQnRENDA/fdd19S5UsFO7ZtY1iCYl6VRf5nWifo9/vZ37ifggQuHQGx8rIh9LpW\nCjlEKBTimb/8PyYUG+TZDyyffL7aT0trG4sWLbJQuuTR3t7Orf99K8qtUEfE6XlTBsYUg1dffTWj\n28cwDHbs2BHrzIdKtJxY6PEMYe9e04vI64o/VHZfeF2FQHZk7tNKIYd4/fXX2bW7jjPHHrzuOa00\nxPhigyefeNxMTZlFKKX45S9/SV1dHaF5oQMhLeIpY5qCCvjVr36Vsd5IdXV1+Pz+IRuZo1Rg7mze\nuHFjgkpMDdFOO9qJJwqXw4PD7owpnUxGK4Ucwefz8cc/PMT4IoPZFcGDjonAl8e3U7dnL88995xF\nEiaHF154gTfffBNjugHlgyzEBuF5YQIS4JZbb8nIuEhr164FYEQf57yMYjewG3gYxcv0PqtyIVTY\nbNRGys0UYkphEEl1+kJE8LqL9ExBkzk88cQT7K1v4GsT27H1YGc8rCzEzPIgjz36SFaMdsAcxd7/\nwP2oEQo1eYhJhTwQmhti29Zt3H///YkRMIWsWLECpwhVfZyzG/BH/rZE3vfFaMNgxYoVGZWbo66u\nDhHB40zsTAHA6yhk967MDxmjlUIOUFtby9NPP83xI/xMKe09Ouo3JnVghAL8+te/yvjMbKFQiNtu\nvw3DYWAcaZAAhxuoBGOywUsvvcTHH3+cgAJTx7+XLGGMUgnxPIpSA7S1t7NhQ+b45+/cuZN8dxF2\n2wDcz+Ik312ScYb3ntBKIcvx+Xz8z09/QrHL4OuT+o6GOtxrcMH4dhYv/pC//e1vKZIwObz88sts\n3LCR0KwQuBNXrjpUIYXCvffdmzEb/hoaGtiydWvcobL7I5qvLJMU5M6dO/G6ipNSdkFeCW3tbRmf\nx1orhSxGKcVdd93Fjh07+c60VvKd/Y/+vzDaz2FlIR64/37WrVuXAikTTzgc5rHHHzNtCAN1Px0o\ndghND7F923bezJBIodGkQRMTXG4hQpUIixcvTnDJyWPHjp0UuBO7RyFK1C0109O5aqWQxTz//PMs\nWrSIL473ceiwgSXVsQlcdWgbBY4gt9z8I5qampIsZeJZunQpDfUNhA8JJ2bZqDujzPDaryzMjCT2\nH374IYU2W5/2hMEyUSlWrlhBe3t7EkpPLK2trbS2tiR8j0KU6Ia4TN+rkHNK4Wc/+xlvvfWW1WIk\nnY8//pj777+fWeVBzh0Xn7dMkUtx/fQW9tXv5Zabb844N9VPPvnEVAZ9udoMBYHwiDBLly5N+yWk\ncDjMko8+YqJhIEnQkBOBsGGYbZ7mRPOGFLiTs3yUH9n7kOn5SXJOKbz22mvceuutVouRVDZt2sSt\nt9zCqPwQV01v69HbqD8mFIe5Yloby1es4M4778wow/OWLVuQIgFHEispgYA/kPYuiLW1tbR1dMTW\n/xPNaMyIqZmQeSzqVedJ8Ma1KA67izynN+O993JOKWQ7DQ0N/Nd//gcu5eM/ZrTgHULHOL8qyJfG\n+1i0aBGPPvpo4oRMMh0dHShHcpWYithn0j2VaTTtZKKNzFEcCGOVYslHHyWphsTR2NgIgMeZn7Q6\n3E4v+/fvT1r5qUArhSzC7/fzox/+gObGfXz/8BbK8obeMZ4zrpPjR/h57LHHeO211xIgZfJxuVyI\nkQxjwgEkLLG60plPP/mEKrFRkBTjisl4YPvOnTQ0NCStjkTQ0tICmLuPk4XL7qG5uTlp5aeCISsF\nEbGLyKci8vdECJRMfvOb31gtQlK56667qK1dx9XTWhlXlJi1bhG4dGoHU0pD/PKOOzLCJ72srAzx\nDawTlKUCTUAT2N60me8HQmSCMGxYIgJRJ4dQKMSKFSsYq5K7uWxc5H+6Z2Lz+XyICHZb8tYVHTZn\n2s8e+yMRM4XrgYwIxr9+/XqrRUgab7zxBgsXLuTscT5mDw/2f0EcOGzw3cPa8NqD/OTH/532huex\nY8didBowAPu6NAkSjPzVC9I0QKXQDMPKhpGfn7yliKGybt06/IEANUmupwpwi6S9UggGg9jEPqA8\nCsGwH4/Hw5e+9CU8Hg/B8MAiCNtsdgKBxD5/qWZISkFEqoEFwB8TI05yySRjaTz4fD7u/819jCsy\nOC9OT6OBUuxSXDGllW3bd/CXv/wlKXUkiilTppgvkri0a2+0M23qtORVkABWrFgBwNgk12NHqFaK\nFcuXJ7mmoTPQxDrBkJ8FCxZw3XXXsWDBAoKhgYeVT2DuHksY6kzhXuC/gF7npyJypYgsEZEl9fX1\nQ6xuaGSrUli0aBH79jfy9Ynt2JNoJZpRHmJGWZBn/t+f03q2MGXKFOx2O9KQpKezE1Sr4rDDDktO\n+Qli5cqVDLPZKEyiPSHKGGDT5s1pnXnM4XAQDocG1A84HW7+8Y9/8Jvf/IZ//OMfOB0D2xZvGAYO\nRzLd3pLPoLsQETkT2KuU+nczm+0aAAAgAElEQVRf5ymlHlJKzVFKzamoqBhsdQkhkwJ3xcObb77J\nqAKDySUD26A2FE6u9tPY3BIbhaYjbrebyVMmY2tIkoaM2FMPP/zw5JSfINasWsWoFP3mqzEHXem8\nC97j8aBQhFX/z4nT7sbn8/Hss8/i8/lw2gemFEJGIK2XFAfCUJ6aY4CzRWQL8GfgZBH5v4RIlSTS\nfaPRYNmwfh0Ti4IpmbZOKjYfqHQ3OM84fAbSKJCEr1z2CU6nk0mTJiW+8ATR3NzM3oYGRqaovmg9\n6awUCgvNyKiBUPJCnwcNPwUFBUkrPxUMWikopX6olKpWStUAFwJvKKUuSphkSSAUSv5I2gr8fj+e\nJPvlR4nWk+45BSZNmoQyFLQkvmxpEsZPGI/T6Ux84QkimgwoUUl1+qMAId9mY/PmzSmqMX5KS80w\nFP5g8pa4/MH2tPZIGwg5tU+haw7ibFIQpSUlNPpT81VG60n3H/7o0aPNF0kIyWPvsDN2TLLNt0Oj\nrq4OgOSEfuuZUqVi9aYj0eXrjkASRgpA2AjhC3Rg9TL5UElIT6KUelMpdWYiykomXf2H03lEEy+T\npkxlfYuLgdrRn6z1sLXVztZWO7ctKeDJ2oFv5lnXZBrRJk5MdMzNxBIdFUog8Wtqyq8oKUlOULVE\nEd1Vm8qFjEKl2J/GG9hGjDCDYbX7k7O5zFQ2KlZPppJTM4WW1lZUJLlGdPt/NjBv3jz2+WBz68AS\nh2xtteML2/CFbaxtcrJ1gNcBLNnrpLSkOO2VQjAY8RVPgp1FRA6Un6ZEvcPiWeDqhIN88+NdIHRA\nWrdLcXEx+fn5tHY2JqX81k5TEY8aleh47aklZ5TC5s2bCfj9KGc+Kr8sY2LhD4QTTjgBp9PB27uS\nG3KhJSB8us/F579wCnZ74jNXJZLt27cDoAoSb2sx8o1Y+elK9PuJx/eoEw7yzY9XKRiAzZa+XYqI\nUFNTQ0tncmYzzb59ANTU1CSl/FSRvt9ggnn55ZcBUM48gsPGs3r1arZu3WqxVImhsLCQk046mffq\nPHQm0VTy1i4XIQPOOuus5FWSIN5//33ELpCEVR6jzGDpsqVpnWErurwVj0klDw7yzc+Ls852oDTN\nbU2HHHIIzb76pOxZamrfw/DhlbnrfZRJtLe389Lf/45y5IHYCFZMBJudZ555xmrREsa5556LL6R4\nry45swVDwRu7vMycOSPtR0ItLS0sfHUh4ZHh+NZPBogaqwgFQ/z97+kb7quqykypE8+m7jw4yDc/\nXqXQaLNRlebr6ZMnTyYQ8tPmT3zyqKbOvUyZMjnh5aaanFAKTz/9NB3t7RiuyKYSp5dg+UT+8Y9/\npP0ywEA59NBDGT+uhrd2JycC5JpGB/UdcPbZ5ySl/ETyxBNP4OvwoaYmyU13GFAFTzz5RNpmppsw\nwcygkKp0L+0omg2D8eOTFaQ7MUybZoYm2deW2OxoncF2Wn2NHHrooQkt1wqyXimsXbuWp556ilDZ\nBLAfGDYGRx2BIXZ+/otfZIV7qohwyqmnsanZxr7OxFtXP97rJM/t5rjjjkt42Ylk3bp1PPvssxjj\nDUhOgi0AwoeHae9o57e//W3yKhkCpaWlVA0fTqqGPNsi/6OdbrpSU1NDfn4+9a2JVQoNESUzffr0\nhJZrBVmtFHbv3s1NP/ghYYcH/9ijDzqmXF58Y49h1cqV3HnnnVkRAuOoo44CYE1j4tdM1jS5mXXE\nLNzugW33t4JQKMQv7/wlyq1QhyV5M18xGJMMXn311bT1ZJsxaxZbbTYUyd/YuBUztlAsGGGaYrPZ\nmDFjBg3tiVWXe1u243a70/7+B0LWKoWtW7dyzTXX0tTaRsekU8D52RXScPkEAqNmsXDhQn6RBTOG\nMWPGYLfb2NWe2K/VULC7XZgw4ZCElptoXn75ZdavW094RhhSkPtGTVNIoXD3PXenZQiVWbNm0WYY\npCI55CYRph96aFoPGqLMnj2bVl9jQvcr1Ldt4/DDD0/rXe4DJSuVwkcffcSV3/4O+1ra6Zh8Bsrb\nu0dEcNQRBKpn8+qrr/K9730/bdeIB4LD4cDpcBBMcNYxQ5l/6fzAB4NBHn70YSgHVZ2iaLh2CE0P\nsX3bdl5//fXU1BkHRxxxBACbklxPB4o6pZg9Z06Sa0oMRx55JAB1LVsSUl6Hv4XmjoZYuZlOVikF\nwzB48skn+c///E984qZ92lkY+WV9XyRCcNQs/ONPYNny5Vx2+eWsXbs2NQInmMbGRjr9AYa5E7sU\nZhcocgu7d6fKbBk/ixcvpnFfI+HJ4aRsWOuVUSBFwosvvZjCSgdGVVUVI6uq2JjkejYDigNKKN0Z\nO3YsFeUV1DUlRl3ubjGjI8ydOzch5VlN1iiFpqYm/uumm/jDH/5AcNh42qeehXIXDvj6UMVEOqae\nSX2Lj6uuvprnnnsu4/IvvPvuuwBMSnAIbRGYVOTng/ffS9sltvfeew9xi5kGLB6CB+/iJd4NuQLh\n6jArlq9Iy30Lc+bOZYsI4STaFTYCnrw8pk6dmrQ6EomIMP+Y+exp3UbYGPrveXfTJoYPr2TcuHH9\nn5wBZIVSWLJkCZdc8k0+/ngJ/pr5+CeceJCn0UAxCipon3YOgYIR3HffffzwRz/KmOWkUCjEn5/+\nE9WFivEJys/clRNGBmhsambhwoUJLzsRrF27FqPUiP8XHTx4F2/cSgFQ5WaHm45ho4888kj8SiXN\nC0mhWG+zccTs2RmVXGb+/PmEwgH2tg6tZUJGkL2tWznmmPkDzuqW7mS0Uujs7OQ3v/kN3/ve92j0\nKzqmnU2octrQ8uE58+icdAr+MUfxwQeLufgb3+C9995LnNBJ4plnnmH7jp18eVx7UvIqzCwPMrEk\nzB9+/yDNzckJKDYU9uzdg8ofxGjYefAu3kFtdotsf9mzZ88gLk4uRx55JA6Hg2QtiNYBTYbBMccc\nk6QaksMRRxyB2+1mV+PQ8oLsbdlGKBzMuPvvi4xVCqtXr+Zb37qUZ599lmDlNNqnndO//WCgiBAa\nMZ32Q8+mKWjnhz/8IXfccQft7UmIw5wAtmzZwsN//COzK4IcUZGcgGQi8K3J7bS0tHDPPXen1dKa\nYRi0t7XDYOzgzoN38Q5KKUTqTcflI6/Xy5zZs1lps2EkYQlpJaab5/z58xNedjJxu93MnTuX3S0b\nh/Rb3tm4AY/Hw6xZsxIonbVknFIwDIPHH3+cq6++mp37mvFNOYNAzXywJ37qqrxldEw7m8DIGbz8\nyit885vfYvXq1QmvZygEAgF++pOfkGcL8q0pyZklRBlTGOa88T7eeONfLFq0KHkVxUksT4ZVMfoi\nP710TTx06mmn0WwYJDpYvIFimc3G3COPTPv8Gj1x7LHH0uFvpbFjcDM8pRR1LRs56qijssIVNUpG\nKQWfz8ePbr6Zhx9+mEDpONoP/SJGcZITDtrsBEcfiW/qmexpbueaa69Nq3X1xx57jI2bNnH5lDZK\n3MkfvZ9V08mkkjD33HN32iyXxIzfVikFMf+i4arTjWOPPZbiwkIWJ7jcNUCzYXBmBgRI7Imjjz4a\nERu7mga3hLS/fTe+QDvHHntsgiWzloxRCn6/n/+66Sbef/99/GOPNo3JjhTsUIpgFFbSfug5BPKH\n8/Of/zwtgqFt3LiRP/3pTxw/wp+0ZaPu2AS+M62NcMC056QDsY1jVtr5hLTdFe92uzn73HOpBfYl\ncAnpAxEqKyoydj29pKSEadOmUtf82TlUiXc4Trsbp91NReFoSrzDP3POrqZNiNiYN29eKsRNGRmj\nFO6//36WLV1K5/gTCFUdOjRj8mBx5NE56VTCxdX8+q67WLVqVepl6MLv/vd/8ToMvjrR1//JCWS4\n1+Ccmnbeeecdli5dmtK6eyLm9WJVn6zMutM5x8R5552Hw+Hg3QSVtxXFVqW48GtfS+v77o/58+ez\nv72OzuDB9sJZY06mxDucEu9wTppyIbPGnPyZa/e0bObQQ6dRVFSUKnFTQkYohdraWl588UWCVYcR\nLrc41ILNTuchJ2M4PNx99z2WGVzXrFnDRx9/zJljOih0xSeDLyQH+eb7QvEr2NNG+ynNg8cefTTu\naxONx+Mxk7v4+z83KUQmaencOZSVlXHGggV8KkJrAmYL7wDFhYWmG28GM3v2bADq43RNDYQ6aezY\nw5wM2cUdDxmhFF577TXE5iAwKk0s/A4X/hGHs379OssS9bzwwgvkOYSTq+PvCTtCcpBvfscglILL\nDl8Y5eOTTz9l27Zt/V+QROx2OyNHjUSaLVo/injojh492pr6B8hXvvIVFPDBEMvZi6IWOP/LXyYv\nL96sC+nFpEmT8Hg8ce9XaGjbiVIqq7yOomSEUqivr0e584dkQ3Bt/QBbxz5sHfvIW/13XFuH9mgY\nHjOz1d69qQg3djB+v5+333qTecM78Q7C6crrUAf55nsdgxs5HjfSj0BaxP2ZM3sOtnobWLDhWnYK\ndoc97cMmV1dXc8yxx/KJzUaoh9nCCEzvWjdQE3nfEx8BToeDc889N1mipoxoZNd4PZD2t+9GxJYV\nUVG7kxFKYfz48eBrQXyD311sa9+HhINIOIi9tQ5b+74hyWRv2nFAthSzdOlSOnydzB0+OG8Xj0Md\n5JvvGaRSKHUrJpWEeO/ddwZ1fSI57bTTUCGFbEjxbMEP9q12jj3mWAoLBx5WxSrOPfdc2g2DNT0c\nOwNhBKYyuAzhjB4s90EUy0Q48aSTYik/M50pU6bQ1LEXQw3cKNXYvpcxY8aYoVGyjIxQCmeffTYe\nrxfPprcgAbFKhoqttQ5X3UpOPfVUysvLU17/J598gsMGU0utb4vDy4Ks37DR8l3O06ZNY+68udhr\n7dCRunplpUAQLr300tRVOgSOOOIIhpWUMFgXiQ1Ap1KceuqpiRTLUmpqajCMMO1xpOhsC+xj/Pjs\niHXUnYxQCqWlpdxy84+QtnrcG98CC3fTiq8Z7/rXGTlyBNddd50lMqxdu4axhWFcaeD0cUixqZhq\na2stlgRuvOFGXDYX9iV2UpBXBnaBbZONCy+8MGOCodntdo49/ng2DDJI3jrAm2U7eMeMGQNAa+fA\nMlqHjRCtvmbGjh2bTLEsIyOUAsBxxx3HNddcg2P/Zpw7llgjRMiPd/0iCjwufv2rX1m2XLB7504q\nPdbPEgCqvOYegbq6OoslgVGjRnH9ddfDHpBVSV5GagfHxw4mHDKByy67LLl1JZjZs2fjV4pdg7h2\ni83GjJkzs2oH78iR5gbYdn/LgM7vCLQCihEjerO6ZDYZoxQALrjgAhYsWIBr13JsranvhFzbPsTm\nb+UXP/851dXVKa8/SjgcxpEmARkdkV9QumQeW7BgAaeffjq2NTYG1esNhDA4PnDgcXq4/bbbcblS\nt4kyEUQN4vFmKe5E0WAYaW9Qj5eSkhJcLteAl486IhnbqqrijdOeGWSUUhARrrvuOkpKS3DtWpba\nuv3tOBvW86UvfYnDDz88pXV3p7yigr2dabB2BOzpMH9CVthWekJE+N73vseEQybg+NiRFPuCfCqo\nRsWP//vHsVFmJlFeXk5RQQHxBimJnn/IIemdljVeRIThwysjM4D+aY+cV1lZmUyxLCOjlAKYG5Xm\nH300Dt/QvIfixdaxD5TihBNOSGm9PXHE7Dmsb3LQ0Gn9dGHxHhd2u43DDjvMalFiuN1ubvvZbbjt\nbuwfJdi+sANsm218/etf5+ijj05gwalDRBhTU0NDnNdFz8/GtfSqqkp8wYEphY5ACyJCRUVFkqWy\nhoxTCqFQiNVr1mA481Nar3KZ9a1Z05MzX2o555xzwGbjLxusdYfb1W7jX7s8fOELp6Sde+KoUaO4\n8YYboR5kU4KUZxAcSzPTjtCd6upqGm3xPf77AbvNxvDhn40DlOlUVVXRERioTaGFYaXDssqu0pWM\nUgqhUIg777yTLZs3469K7cjU8A4jXDyK3//+IRYvTnS8yfioqqri4ou/wft1bt7eZc16dmcY/ndV\nIXkeL1deeaUlMvTHaaedxoyZM7CvsveaUU2VKJQz8lehUCW9TytkraA6FTf9100ZlWWsJ6qqqmgx\njB43sfVGE1BRXp7x994TVVVV+ALthML9B5Zs9zdTlaVGZhiCUhCR0SLyLxFZLSKrROT6RArWnd27\nd3Ptd7/LwoULCYyaRbgsxZvGROiccCIBVyE33XQTjz76qKX5ir/xjW9wxKxZPLw2n5X7UvuQhgz4\n7YoCtrXaufmWW9PGntAdEeGq71yF8qteZwtqpoISoASMEw3zfU8Ewb7RzkknnZQVu1grKytRQDxp\ngZohazvDmAdSoP/9Nh3BZqqrRyVbJMsYykwhBHxfKTUNOAq4RkSmJUasAxiGwfPPP883vnEJa2rX\n0znhJILVsxNdzcBweuiYuoDgsAk8+uijXPntb7Nhw9DS+Q0Wh8PBz267jbFja7h7eSFrG1OjGAwF\nD67K59MGJzfceGPaZ9yaNm0a0w+bjn3z0GwLsk1QQcWFX7kwccJZSNRIGk+MgGabjcos9bgZNcrs\n5Ns6G/s8L2yEaO9szUgHg4EyaKWglNqtlPok8roVM+dGQtXnrl27+O53r+Pee++lI6+M9unnES6f\nkMgq4sfuwn/IiXRO/Bwbt+7k8iuu4JFHHrFk1lBYWMjd99xL1chqfr2siA3NyfVIMhQ8tNrL4j0u\nrrrqqoyJfXP6aaejWlUscN1gkB1C9ejqrJglwAGlMNAmCaNoMYys9biJupi39eOWah5XlrqkJ5uE\n2BREpAaYBXzYw7ErRWSJiCypr68fcJn//Oc/ueSb32Tlmlr844+nc/JpKHdBIsRNCOFh42ibfh6B\nknE89thjXH31NZZkIhs2bBj33Hsfwyoq+dWyYna1J89M9Kd1Ht7d7ebSSy/lq1/9atLqSTRHHXUU\nALJ3kAbnMNgabMw/ej5iRR6PJBA1Fvc9Lj5AK2a6imw0MoM5wCooKOx3phDdy6BnCn0gIgXAc8AN\nSqnPmO+VUg8ppeYopeYMxIXLMAx+97vf8dOf/hSfs5j26V8kVDHJmqQ6/eHMM2cNh5xM7YaNXHrZ\nZSxfvjzlYlRUVHD3Pffi8hbx62XFdPQzaRlbGMZjN/DYDaaUBBlb2P/Gs3/tdLFwex7nn38+l1xy\nSYIkTw0VFRVUDK8w3WcGQzMoQ2XVpi23201pScmAl4+iXWW27uIFGDliBG3+vudObZ3m8ehyUzYy\nJKUgIk5MhfCUUur5oQoTCoX42c9+xtNPP01w+FR8UxYkbnYQDhyUWIZw4vLphsvG0z7tHFqDNm64\n8UbefTdR+a0GzsiRI7n957+godPGY2u9fZ578WQfYwvDjC0Mc8ucNi6e3Hfmtt3tNp5YV8CcObO5\n9tprM3K0PGniJOwtg1tei+ZpmDhxYiJFspxR1dUD1pPR87J5hFw1oorOUN+m945AMy6Xi+Li4hRJ\nlXqG4n0kwMPAGqXU3UMVJBQK8dOf/pR//vOfBEYfSaBmPsTpR90XEgoclFhGQolNsq48xbRPPZOg\nu5Rbbr2Vd95JfTjp6dOnc/HFF/N+nZt1TYmzL/xpvRenK4+bb74lY1MvjhkzxrQrDMbY3Goa9rMt\nrEE8exX2Ye5RyFabApgzyg5/30rBF2ijorwiIwdGA2Uove4xwMXAySKyNPJ3xmALe+CBB3jrrbfw\nj5lHcOSMhC8XKYfroMQyaggJe3rFmUfH5NMIe8r48U9+wurVqxNfRz987Wtfo7iokFe2JiYj1q52\nG582OPnq175GWVlZQsq0glGjRpmL4oMIeyFtQmVVZcYqxN6orq6m2TAIDEBT7sNcXsnGPQpRysrK\nCIYDhPsIz98Z6mBY2bAUSpV6huJ99K5SSpRShyulZkb+Xh5MWUuXLuX5558nWHUooRFJ2pRmdx2U\nWAZ7kjZ9OVx0TDqFkD2P226/PeVeSR6Ph9NOP4NPGlz92hYGwgd1LmwinHXWWUMvzEJia8DtfZ/X\nE7YOG6Or0zvV5mCIhoweSLiLfTYb1ZHzs5Vo1ONAqLPXc0KGPyOSKQ2FtNjR/PLLLyPOPAKjj7Ra\nlMTgzKNz1Bx2bN9uSViMY445hrCCNfuHvg1/+X4XU6dOZdiwzB4dRdfCpS3OGagC2rJzLT2aU7q/\nKGIGiv1KxZRIthLNohYyet/VHDKCWZltrStpoRRCoRBIWoiSOCLLX1bsX5g6dSpOh531zUOb6gfC\nsKXVzoyZMxMkmXUMHz7cXPpoi/NCP6hgdvqlRxVdf8bmViCoVFZ73ACx5UFD9e6Np5SR1UtokCZK\nYcGCBaigD/f6f0KCDcBWYGutw7P1fUZVV1sSPdTtdjN+/Hg2tQ7tx7u9zU7YMJVMpmO32xk5aiTS\nEudMIeJknY2jZI/HQ1FhYb8b2KLHs9nIDKAiGR2lr24xiw3MUdJCKcyePZvvf+97uFp2kr/qBez7\nt1iacnPQhAI4t32IZ80/qCwv5Vd33mnZqGLS5ClsbXUOqRk3t5ojp0mTJiVIKmuZMH4C9tb4jMVR\nJZIp6TbjZVhpab9mlujkKtOXEPsjGDSXjey23n8jdrHj9/tTJZIlpIVSADMc9AMPPMDYyjLy1r+O\nZ+3L2FpTv0N4UBhhHHUrKVjxLK66lSw44wwefeQRS5ccpkyZQntQscc3+K94c4uDosKCrHHFnDBh\nAkab0WvE1B5pgvzC/LQN+jdUCouK6N2sahLtArPdwNrWZqo/Rx9OKA6bi46OJGRuSiPSanFs+vTp\nPProI7z44os8+tjjNK9+iXDJaALVczDy09AdUhk46tfj3vUp+NuYMXMWV131nbRYbonG6NnUYqfK\nawyqjE2tLqZMnZY1PtmxjGHNwAD7eFuLjYkTJmZNG3THZrf365AaPW5L4L6hdKSlpQUQnHZ3r+c4\n7Xk0Nw0hiFYGkFZKAcxNQueddx6nn346zz33HE899SfaV/6VYMVE0zvJmR6Wf1vzLjzbFkPHfiZP\nmcK3r7ySOXPmWC1WjHHjxuF2OdnY7GB+VTxDY5POEOxoEz43LeGBby0jugQkzYIqH8C6mjKXj8aP\nT3GY9hQSDAT6XS6IHo8ur2QrTU1N5Lk82PpwenE7PDQ2xpuzLrNIO6UQxePxcNFFF3HOOefw5JNP\n8swzz+Bq2o6v5hjCw2riLs/ILzNTagKGt2zwMw8jhGvrYpx711JZWcU1P/gfTjjhhLQbSTocDiZN\nnsym7cuAvsNY9MTmVgdKZYeROUplpbkBzWgf4Mwpiz2PorS2tFDUzznRYVhrazzZFzKPxsZG3I6+\nQ8S4HV527Y8n4HjmkfbzwcLCQq6++moeeeQRJo4bQ97613HuXBp3OYGxR5vKwFtG57QzCYwdRH7d\nYCfeNS/jrK/lK1/5Cv/3f09y4oknpp1CiDJ58hS2tTowBmFs3hoxyE6ePDnBUlmHzWajuLT4wCJ5\nf0TOy2YDa2NjI/1FF4smvm1qyu7OsK2tDae970gALoeHYDBAIJD5XpK9kfZKIcq4ceP43//9LZ//\n/Odx7ViCo25lagUwQnjXvYrT38jtt93GNddcg9vd+9pjOjB+/Hj8YUXDIIzNO9rslBQVZl2H6HK6\nzHAXAyFyXrbm4g0EArR1dPSrFKLm5X37+tvmltm0t7XjsPX9XTsjRuhsNjan7fJRT7hcLm6++WZ8\nvk7ee/99jIJKjIL+w3EnpO5tHyFt9fz09ts57rjjUlLnUImGOW7otDE8TmPzvk4bI7JwF6+v0wcD\nDQsV8Uzs7OzPPyczMQ2r0PeCyYHlo+bm7Dawiq3/Gb+KmN3TdXUgEWTMTCGK3W7nRz/6ISUlxeRt\n+yAl+xls7ftw7lnN+eefnzEKAQ64EHaE4v8Bd4RtFBZlV3hgwzBobWmFgU7wIudl67KJz2famvpr\nDheCdDk/W3G73X2GuABiwfJcriTFTksDMk4pgNnZXX7ZZUjrXmwtu5Nen3P3cjweL5deemnS60ok\nURfCwahNpSTrRkNtbW0YYWPgMwUXINmrFKJhHfpLsWRExsfZFiW2O5WVlfiCn8kTdhDt/iYKC4uy\nOv5RRioFgFNPPRWPNx9Hw/rkVhQO4GzcwmmnnZpxm3fCYfNxH8yXLKIwjMHtb0hX2tsje3cHumgq\nYHPZYpuaso2SkhLAjG3UF9HjpaWlSZXHasaMGUN7ZwuBUO+eCC2+fYwZk30Rc7uSsUrB7XYz/+ij\ncLXsTOoSkr1lN8oIc/zxxyetjmQRnSkMpmtXZN9mpdiUv//soyYKVEiRl5eY3BTphtfrpbysjP7i\nBkSPjx07NtkiWcrhhx8OwN7WrT0eD4b97G/fzcwsCBDZFxn91M+ePRsV6EB8yZve21t243A6MzI/\nb1GR6YHe5I//a24O2GPXZwvFxcW489yxIHf90g4qrLI6L/HsOXPYaLMR7mORcT3gdjqZlkUbGXti\n+vTpeD1edjVt7PF4XfMWDGUwb968FEuWWjJaKRx5pJl/wdHYs2YfMkrhbN7OzBkz0979tCeGDx9O\nSXERa5viczLb22Fjny97AuFFcTgczJwxE3udfUDTJ9ll2lRmzZqVZMms49hjj6XDMOi5G4QwilU2\nG0fOnZu1M6YoDoeD4084np1N6wmFP2tw3rpvFWXDyiyJfJxKMlopVFZWcviMGbjr10I48XkL7I1b\nwdfMaaedmvCyU4GI8LnPf4El9S4aOgduNF643Y3NZuOkk05KonTWcPbZZ6PaFbKjn/YIg32jnemH\nTc/qZZP58+dTWlzMh70cXw20GgZnnX12KsWyjNNPP51gyM+OxnUHfe4LtlHXvJlTTj0l6w3uGa0U\nAC6/7DKUvw3X1sS6p0qgHc+2Dxg7toaTTz45YeWmmgsvvBCn08Xja/MH1DxbW+28viOP008/neHD\nhydfwBQzf/58xk8Yj32lHfoYR8h6QbUpvvXNb6VOOAtwOp2ce955rAP2dFtCUijeE2HUiBHMnTvX\nGgFTzIwZMxg5ciSbGurnARUAAA6NSURBVJYd9Pnm+pUYyuDMM8+0SLLUkfFKYebMmVx00UU462tx\nbfswIYpB/G14axfikjA//vF/Z3SmpcrKSi6/4ko+bXDy2va+l8B8IXhgZSGlpaV85zvfSZGEqcVu\nt3PD9TeYs4XVvcwW2sC+xs78Y+bHliizmfPOOw+3y8W73T7fBOxUiq9ddFHWj46j2Gw2zjnnHBpa\nd8aWkJRSbNm3glkzZ8VSmGYzGa8UAC6//HK++MUv4qxbSd66RRAc/A5UW/Mu8le/iEf5+dWddx4I\nt5zBfOlLX2L+/KN5ar2X2qaeH25Dwe9X5bPHZ+fHP/kpxcXZtXGtKzNnzmTBggXY1tmgu4+CAvsn\ndtxONzfecKMl8qWa4uJizjzrLJZz8OTpXYTS4mJOOeUUq0SzhFNOOQWbzUZn0HRhbmjbQVtnEwvO\nXGCxZKkhK5SCzWbjhhtu4Prrr8fVuouCVS9gb9oRXyFGyMyatvYVRlWW8eCDv8sa1zObzcbNN99C\n1YgR3L+iiGb/Z0fIL291s6TexVVXXZU1990XV111FcXFxdg/sR+0u092COyBb1/57axPP9mVL3/5\nyyAS25NQj2IDivO//OWMdLIYCmVlZcydOxd/qAOUYtu+NeS58zIqmsFQyAqlAKZR9fzzz+fBBx+k\nurKMvNqFuDa/NyADtK19H/mrX8S1ewVnn30Wf/zDH7Iu/WJhYSG33f5zOgwnD64+OATaxmY7f9no\n5YQTTuCCCy6wSMLUUlRUxNVXXQ37OBA1NQz2lXbGTxjPueeea6V4KWfkyJHMnTuXNkwd+Qlgt9lY\nsCA3RsfdOf744zFUmJARZHfLJubOm5vVu5i7kjVKIcrkyZN55OGHueCCC3DWryV/9d+QjsaeT1YK\nR91qvKtfpNQFv/zlL/mP//gPvN7+QoRlJhMmTODqa65hxT4HzQHzqw8b8Ie1ZjTUm266KetCW/TF\nKaecQk1NDdIhZkKdbaZx+arvXJUza+hdOf2MMwhj6shVNhtz582jrCwNMx6mgKOOOgqAkBGiw9/K\n0UcPItR+hpJ1SgHM3c7XXnstd991F8UuRf6alz67nKQMXFvew731febNPZLHH38sJ774c889l6lT\nJlPfaaZhfHu3ix2twvU33EhBQX9BlLMLu93OV7/6VXOHcxDsm+yMGz8uZzxtujNv3jxEhE6g0TBy\nZrmkJ8rLy6mqrKK10wwXnu17E7qSlUohypw5c3j4j39k3JjReNa/BqFIYgylcG1+F+fetVx00UXc\ncccdsTgw2Y7NZuPyK64kZEBLwMYr271MmjQxI8N4JIJYkqQOUPsVZ5x+Rk7Nlrri9Xo5dNq0mO39\niCOOsFQeq5k6zcw66PF4sjr7XneyWimAuav3N7+5j+rqUdg7m0AZOOprcdav4+KLL+bKK6/Muhg/\n/TFnzhxcLid7fTZ2tQnnnHNuznaEHo+HgsICJGjefy7MFvtiSiT9aoHXm9XhPQZCVBGMGDEip/qI\nnLjToqIibr/tNlAG0tlK3o6PmTFzJpdddpnVolmCiFBUVExYmR3h/PnzLZbIWvK9ZsJJu8OeE37o\nfRG9//zCwpwdKESJKsVszbzXGzmhFABqamoYNmwYtpAPFfRz5RVX5JT2705+fn7sda4aE6NEY/rY\nbfac7wijbri53g5wILR4rrVFTvWKXTdkZWLU00SSzZmj4iXqfpzLg4Qo2bxpMV6yLUrwQBnSUyAi\np4lIrYhsEJEfJEqoZNHVzzjXtH93Jk6caLUIacPll18O6N8EHJhB6rYgFt4m19pi0EpBROzAb4HT\ngWnAV0UkrQOuT5482WoR0oZrrrnGahHShlzZlDQQojNIlYLc55lCrrXFUGYKc4ENSqlNSqkA8Gfg\nnMSIlRxuuOEGq0VIG/Ty0QF0W3yWTA4CmWhyrS2GohRGAdu7vN8R+ewgRORKEVkiIkvq6+uHUN3Q\nybVpYF/k2g+9L3Jx93JvlJaWMnz4cK644gqrRbGcsWPHMm7cuJzzUkx6z6CUegh4CGDOnDm5NQ9L\nY7RR9QB6sHAAj8fDs88+a7UYaUF+fj6PP/641WKknKH0DDuBrk7d1ZHPNBmA7ggPkGtrxhpNXwxF\nKXwMTBSRcSLiAi4EXkyMWMkl1zdraQ4mqiCzLSe1RjMYBr18pJQKici1wKuAHXhEKbUqYZIliUce\neSSn4uT3x1lnnWW1CJbj9Xq57rrrcj7EhUYDIKmcOs+ZM0ctWbIkZfVp+mbXrl2UlZXlXBIVjSbT\nEJF/K6XmpKIu7YKSw4wcOdJqETQaTZqhXVA0Go1GE0MrBY1Go9HE0EpBo9FoNDG0UtBoNBpNDK0U\nNBqNRhNDKwWNRqPRxNBKQaPRaDQxUrp5TUTqga0pq7BnyoEGi2VIF3RbHEC3xQF0WxwgXdpirFKq\nIhUVpVQppAMisiRVOwPTHd0WB9BtcQDdFgfIxbbQy0cajUajiaGVgkaj0Whi5KJSeMhqAdII3RYH\n0G1xAN0WB8i5tsg5m4JGo9FoeicXZwoajUaj6QWtFDQajUYTQysFjUaj0cRIC6UgIt8UEcszvojI\niSLy98jrs0XkB5HXFSLyoYh8KiLHiciXRWSNiPyrv3IGIUOPbSEit4vIdhFpG0y5g5AjLdtCRLwi\n8g8RWSsiq0TkjsGUHaccadkWkc8XisiySFs8KCL2wZQfhxwpb4uB9A8iUigiS7v8NYjIvfHeXz91\nWP476EO2PvsHETlfRJSI9L/nQill+R/wJjAnDeQ4Efh7D59fCPyxy/uFwLHxljOUtgCOAkYAbbnc\nFoAXOCny2gW8A5yei20R+bwo8l+A54ALs60tBtM/AP8Gjs/0e49Dtl77B6AQeBtY/P/bO9cQq6oo\njv+Wj2YcfGQmZqRNWCiZFWYPU0JMspelEKmR4aeIgvqQGPSGLNMgtYdlkRVYmT0sK0uiNHAMTXR0\nxhdjKmklE+RjVHyEqw9r3Znj9Z4754zO3GHu/sPh7HPO3nvt9d9rr3X23od7k/CYVvBXTvYm4EHg\nIeCVyPPJwBuefgbYBqwEPgGmxNR5D3DI81YCnYBrgF9c1jKgd8Q4ZgFrgS3AtcCXQA0wzfOUA1uB\njzzP50BZHp1u9fzrgNcynZXRBbga+AP4x9v3nLe3Dtibiwvv9Eovf52XPwocBL4FpmDObRGwGVgM\nrAamJuDicOCiwS6AOZ6n2LnoA3zjcgvFxTHgCFDreu0F1sT4h/dcp5Xe/l3AemAV0N/zlPn1Sef1\nMDAslx1kjYkq4ESB7WAbEd8YFxTcDn6N0T3bDhp36LmDwmzgDhIG17RB4Tw/dwKqgV7A9sjz74Hh\n3hmVQCkWpWqICQqRzhzi6Y5OTk+/Hg/Mj+Sb4enHgL8wx1AC7AF6eKcrMMzzzY+T7e3bDVyGvWkt\nyu707HSkHSPjuPBOr3UuRkS4GOMGM8WPeZ7/CuA/YEgCLk4ELuq5WADswN6AipmLDZhD/pjC2sW7\nEd0PA7cQ7x+2A0sx//A7MNXzjAK+8PQUYJ7rNN65uJ5G/APwLPBzC+uebQexzpdTg0JXoEOc7tl2\nkDYoAIMjdeZtV+ZIu6fwqIhswAZhH+ASYIeI3CAiPYABQAUWzb9W1aOqWoe9wSRFfyfhRxGpBJ4G\nLoo8X+LnKmCTqv6tqscw59DHn+1W1QpPL8AMMRcGADtVtUaNtQUp2jkhjgusoztjXIzAjHIt8BJw\n3MsPBxYCqGo1sDGHjFxcSOR5sXNxK/b2dpTi5qIdtoxWAnSncFxcFdG9BLiQeP+wEjjp/uEnYJKI\nVGNv+gOzucACx0bMqTfmHyYAn7aw7k1FN+CzfLrnsYO8EJF2wKvA42nKdUghYAQWyYaq6hERWYFF\n0oXAvdgUa7GqqojE1pNEFNaZQ2OeH/PzyUg6c53RR7PKZF+fKc7Fpn1xXJwE9joXtwO7VPU2ESnH\n3uqS4jQusjaSipYLEZmPvRXNFpGxFDEX9TdFHgCmUxguBmDO+VLXfT8WGJL4h2FAjaqOdS5WNCIr\nn3/oi+m4jcLZQRq8ACxX1XEJdU+DLlgAXeGcXwAsEZG7VHVtXKE0M4VuwD7v8AHYxgbYetfdwEQa\nonoFMEZESkWkM3BnI3XXuQJgndlTRDKDv6OIDIwtmRt9M+WB+7C3klzYCpSLSD+/npiw/vbAwTxc\n3IxNXcGm0f1FpBRbYy7z+xXYYEFELgcG+f28XJD+i7E2x4WITMMccNqfIGhrXIwUkd4+Rq7E1o2P\nJGzr2eaiE3AsontXvx/nH24E2rl/6Asc8GeTI3VmuKjDnNsg7Kf38/mH0dgeZj40tx2kQTfgT09P\njtyPs4PEUNUDqnq+qparajk2g8sbECCdg/kB6CAiW4CXXQCqug/bsLlYVdf4vd+wadxGbB2xioZO\nz4UPgLd9Otge23ye4VPRSsyA0mAb8Ii3tTvwVq5MqnoUG5Dficg6bL03Cf4F2ufhohew3/M+iXX8\nfuB+bJngADAXM+7NwDRsg/AA8VzUYgO+VET2YNPoJGhrXMwCnsIc4FJ/1jthW9saFy9iyyp12Dir\npSHoNIazzUUV5uQzuh/M0j3bP6wCbvJ2rwVGi8h6Tl29mAv0xDZzX8fW1Q+R3z+MovGg0Nx2kAYz\ngelxuuewg5wQkZnuF8pEZI+IPN/kFjW26dDUA+js5zKs0wc3l6wsueVAdUvIOhMusIFd6vf7ATuB\ncwIXgYti4CKJf2gOLlqD7gnb2ex2EHck3lNoAt7xaU8p8KGqrmtGWa0dp3EhIl2A5b4kJMDDqno8\nby1tA4GLBhQzF0n8QxnFwUUuFEz3lv47zjexTaUo5qjq+y0gezH2JUgUT6jqsjOoczT2CVwUO1V1\nXIKyObnAvmsvybo/SVWrmtrOHLIDFw2yAxcNsgvKRVL/ICKrOctcFFr3FHWedd1Pk9GSQSEgICAg\noHWjVfz2UUBAQEBA60AICgEBAQEB9QhBISAgICCgHiEoBAQEBATU43848teYVNw/vAAAAABJRU5E\nrkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "# Choose stations with lots of samples (this does not decrease our dataset that much)\n", "df_subset = df[(350 < df.cnt) & (df.cnt <= 365)]\n", "plt.title(\"Average absolute temperature difference\")\n", "sns.violinplot(data=df_subset[[\n", " \"avg_temp_diff_lag_1\",\n", " \"avg_temp_diff_lag_2\",\n", " \"avg_temp_diff_lag_3\",\n", " \"avg_temp_diff_lag_7\",\n", " \"avg_temp_diff_lag_14\"\n", "]]);" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "4HAcGpA5PMJ7" }, "source": [ "We can see that for one day ahead this is approx. 2\\*C on average, increasing to about 4\\*C if we would forecast 14 days ahead using persistence. The lumps in the violin plot are probably due to discretization (and not at rounded numbers as we converted from F to C).\n", "\n", "We are starting to get a \"feel\" for these numbers, and can now intuïtively say that a model that has MAE of 1\\*C for dayahead forecast on the same subset is doing quite well compared to this \"persistence\" baseline model." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "joW3UKMdRjiD" }, "source": [ "Let's have some more fun and find the station where the temperature tomorrow is most similar to today's, and where the temperature in 2 weeks is most similar to today's." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 131 }, "colab_type": "code", "id": "EGxctXZUPYU_", "outputId": "904e99ee-3b10-4026-e5c4-b7f22f53d598" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stnlonlatelevcntavg_temp_diff_lag_1avg_temp_diff_lag_2avg_temp_diff_lag_3avg_temp_diff_lag_7avg_temp_diff_lag_14
1072889578077.02-76.420+2830.03580.0949270.1562110.2266040.4979420.949451
10709914900-157.351.986+0001.53580.3431370.4422600.4625980.5210510.541990
\n", "
" ], "text/plain": [ " stn lon ... avg_temp_diff_lag_7 avg_temp_diff_lag_14\n", "10728 895780 77.02 ... 0.497942 0.949451\n", "10709 914900 -157.35 ... 0.521051 0.541990\n", "\n", "[2 rows x 10 columns]" ] }, "execution_count": 15, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "df_subset.loc[[\n", " df_subset.avg_temp_diff_lag_1.idxmin(),\n", " df_subset.avg_temp_diff_lag_14.idxmin()\n", "]]" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "HotfEFE-SjJ2" }, "source": [ "The first is a [station in the antarctic](https://www.google.com/maps/search/?api=1&query=-76.42,77.02), the second [Kiribati island](https://www.google.com/maps/search/?api=1&query=1.986,-157.35) in the pacific near the equator. This again makes sense, as the antarctic probably has less weather systems and energy influx to change the temperature drastically overnight, while on an island at the equator surrounded by sea we have the more stable long-term temperature." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "odbRBN5SNZDn" }, "source": [ "So now we've done some good basic analysis on this dataset, and are ready to dive in deeper. Next post will be about exporting this dataset and doing some epic ML on it." ] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "1. Exploring global weather data.ipynb", "provenance": [], "version": "0.3.2" }, "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.7.3" }, "nikola": { "date": "2019-07-28 12:00:00 UTC", "slug": "global-meteo", "title": "Exploring global weather data" }, "pycharm": { "stem_cell": { "cell_type": "raw", "source": [], "metadata": { "collapsed": false } } } }, "nbformat": 4, "nbformat_minor": 1 }