{ "metadata": { "name": "", "signature": "sha256:1dec014908b350d0c54505ef3415b91eabd06f7a0dacf1eba412ea9fed14afea" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Weather Changes: an Introduction to Data Analysis with Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a basic introduction to real-world data analysis and visualization with Python. To see the finished product and a bit more analysis, read my [blog](http://corbt.com/posts/2014/05/10/climate-change-since-1950.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step in data analysis is finding the right data. For this project we'll be using the Global Historical Climatology Network (GHCN) monthly data compiled by the US government. In a [previous post](http://corbt.com/posts/2014/05/01/world-temperatures.html) I explain where this data comes from.\n", "\n", "A professor of mine, and an experienced data scientist, once mentioned that \"80% of the effort in machine learning is getting the source data into a manageable format.\" The following section simply does exactly that. :) The source files are in a custom fixed-width-column format specified on the [NOAA website](ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v3/README).\n", "\n", "The first block here parses three \\*.dat files, one each for minimum, maximum and average temperature measurements. Each of these files consists of one line per station per year. A line includes the station ID, year of observation, and measurements for all 12 months.\n", "\n", "These files are simply read in, split at the correct locations, and written out to a CSV file. This will allow us to import the data into pandas directly using pandas's `read_csv` functionality." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import division\n", "from glob import glob\n", "import csv,re\n", "\n", "for name in glob('data/raw/*.dat'):\n", " label = name[15:19]\n", " print \"Parsing data:\",label\n", " with open(name) as in_file:\n", " with open(\"data/csv/\"+label+\".csv\",'w') as out_file:\n", " writer = csv.writer(out_file)\n", " header = [\"id\",\"year\"]\n", " header += [str(month+1)+\"_\"+label for month in xrange(12)]\n", " writer.writerow(header)\n", "\n", " for line in in_file:\n", " # if int(line[11:15]) < 2000:\n", " # continue\n", " components = [line[0:11],line[11:15]]\n", " components += [int(line[19+x*8:24+x*8])/100 for x in range(12)]\n", " writer.writerow(components)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Parsing data: tavg\n", "Parsing data:" ] }, { "output_type": "stream", "stream": "stdout", "text": [ " tmax\n", "Parsing data:" ] }, { "output_type": "stream", "stream": "stdout", "text": [ " tmin\n" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "A `country-codes` file is provided, which correlates each station ID with a country. Because this file is small we can read it straight into a python dictionary of the form `{[country_code]: [country_name]}`, which we will use in the next step to associate countries with station IDs." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def title(name):\n", " \"\"\"Title case with correct handling of internal apostrophes\"\"\"\n", " return re.sub(\"([a-z])'([A-Z])\", lambda m: m.group(0).lower(), name.title())\n", "\n", "cc_file = open('data/raw/country-codes')\n", "ccodes = {l[:3]:title(l[4:].strip()) for l in cc_file}" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Index files, ending in \\*.inv, correlate the station IDs in the .dat files with station metadata, such as latitude, longitude and station name. The following lines open each .inv file and convert them to a CSV, as above. It also uses the ccodes index from above to look up the country name for each file." ] }, { "cell_type": "code", "collapsed": false, "input": [ "for name in glob('data/raw/*.inv'):\n", " label = name[15:19]\n", " print \"Parsing index:\",label\n", " with open(name) as in_file:\n", " with open(\"data/csv/\"+label+\"-inv.csv\",'w') as out_file:\n", " writer = csv.writer(out_file)\n", " writer.writerow([\"id\",\"lat\",\"long\",\"name\",\"country\"])\n", " for line in in_file:\n", " # if not line[74:79].strip().isdigit() or int(line[74:79]) < 100:\n", " # continue\n", " name = title(line[38:68])\n", " name = name[:name.find(\" \")].strip()\n", " components = [line[0:11],line[12:20],line[21:30],name,ccodes[line[:3]]]\n", " writer.writerow(components)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Parsing index: tavg\n", "Parsing index:" ] }, { "output_type": "stream", "stream": "stdout", "text": [ " tmin\n", "Parsing index:" ] }, { "output_type": "stream", "stream": "stdout", "text": [ " tmax\n" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aggregating the Data\n", "\n", "Now that all the relevant data is stored as easily-parsed CSV files, we can begin processing it!\n", "\n", "The main tool we'll be using for this is a pandas `DataFrame`. Think of a dataframe as much like a table in SQL or even a spreadsheet in Excel with lots of data indexing, retreival and aggregation tools built right in. For more information on DataFrames, see the pandas [introduction](http://pandas.pydata.org/pandas-docs/stable/dsintro.html). Pandas lets you build a dataframe from a dictionary, but in this case we'll be reading one right from a CSV file on disk. It looks something like this:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import pandas as pd\n", "\n", "tavg_dataframe = pd.read_csv(\"data/csv/tavg.csv\")\n", "tavg_dataframe.head()" ], "language": "python", "metadata": {}, "outputs": [ { "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idyear1_tavg2_tavg3_tavg4_tavg5_tavg6_tavg7_tavg8_tavg9_tavg10_tavg11_tavg12_tavg
0 10160355000 1878 9.60 10.2 11.80 16.80 20.5 23.10 25.60 27.50 23.9-99.99 14.40 12.20
1 10160355000 1879 12.50 12.4 12.90 16.20 16.3 23.40 24.70 25.80 23.1 18.20 15.20 9.70
2 10160355000 1880 10.30 11.8 13.10 15.90 17.8 21.40 26.50 26.40 23.6 21.30 16.40 13.60
3 10160355000 1931-99.99 10.4-99.99-99.99 19.2 24.60-99.99 26.70 22.3 20.00 16.20 11.30
4 10160355000 1932 10.80 10.5-99.99 14.90 19.1-99.99 23.60-99.99 25.1-99.99-99.99-99.99
\n", "

5 rows \u00d7 14 columns

\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 4, "text": [ " id year 1_tavg 2_tavg 3_tavg 4_tavg 5_tavg 6_tavg 7_tavg \\\n", "0 10160355000 1878 9.60 10.2 11.80 16.80 20.5 23.10 25.60 \n", "1 10160355000 1879 12.50 12.4 12.90 16.20 16.3 23.40 24.70 \n", "2 10160355000 1880 10.30 11.8 13.10 15.90 17.8 21.40 26.50 \n", "3 10160355000 1931 -99.99 10.4 -99.99 -99.99 19.2 24.60 -99.99 \n", "4 10160355000 1932 10.80 10.5 -99.99 14.90 19.1 -99.99 23.60 \n", "\n", " 8_tavg 9_tavg 10_tavg 11_tavg 12_tavg \n", "0 27.50 23.9 -99.99 14.40 12.20 \n", "1 25.80 23.1 18.20 15.20 9.70 \n", "2 26.40 23.6 21.30 16.40 13.60 \n", "3 26.70 22.3 20.00 16.20 11.30 \n", "4 -99.99 25.1 -99.99 -99.99 -99.99 \n", "\n", "[5 rows x 14 columns]" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we're able to open our CSV files and create dataframes, we can do some basic number crunching with them. For this project we're not interested in the month-by-month breakdown of temperatures, we just want the average for a year. So for each of our dataframes, we need to take the mean of the months. That's done with the line \n", "\n", "```python\n", "typ_data[typ] = typ_data.drop(['id','year'],axis=1).mean(axis=1,skipna=True)\n", "```\n", "\n", "typ_data is a dataframe of the schema we can see above. For each row, we want to average together all the columns *except* year and id. That's why we call `drop` -- that removes those two columns for the rest of the calculation. We then take the mean of the rest of the columns, and store it back in the dataframe as `typ` (a stand-in for tavg, tmin or tmax). The argument `axis=1` means that we want to operate across columns, instead of rows which is the default.\n", "\n", "We then remove all columns except those we're interested in, with the syntax `typ_data = typ_data[['id','year',typ]]`. This assigns a new dataframe to `typ_data` composed only of the columns from the original included in the double brackets.\n", "\n", "Next, we merge all three dataframes together one by one. This works just like a join in SQL, and gives you similar options. We need to make sure that all three readings are from the same year are matched together, so we join on `id` and `year`. By default, this performs an inner join, so only stations and years with all three types of readings are included in the final dataframe.\n", "\n", "Finally, we merge the completed dataframe with the metadata index, giving us a latitude, longitude, station name and country for each data reading." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "\n", "data = None\n", "\n", "for typ in [\"tavg\",\"tmin\",\"tmax\"]:\n", " # -99.99 is used in this dataset to represent missing data, so we replace it with NaN (not a number)\n", " # This keeps it from being used in calculations.\n", " typ_data = pd.read_csv(\"data/csv/\"+typ+\".csv\").replace(-99.99,np.NaN)\n", " typ_data[typ] = typ_data.drop(['id','year'],axis=1).mean(axis=1,skipna=True)\n", " typ_data = typ_data[['id','year',typ]]\n", " if data is None:\n", " data = typ_data\n", " else:\n", " data = pd.merge(data, typ_data, on=['year','id'])\n", "\n", "index = pd.read_csv(\"data/csv/tavg-inv.csv\")\n", "data = pd.merge(data, index, on=['id'])\n", "data = data.dropna()\n", "\n", "data.head()" ], "language": "python", "metadata": {}, "outputs": [ { "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", "
idyeartavgtmintmaxlatlongnamecountry
0 10160355000 1996 19.370000 15.860000 24.200000 36.93 6.95 Skikda Algeria
1 10160355000 1997 18.790000 14.942857 21.828571 36.93 6.95 Skikda Algeria
2 10160355000 1998 18.550000 14.730000 22.460000 36.93 6.95 Skikda Algeria
3 10160355000 1999 19.536364 12.914286 20.600000 36.93 6.95 Skikda Algeria
4 10160355000 2000 19.070000 16.111111 23.933333 36.93 6.95 Skikda Algeria
\n", "

5 rows \u00d7 9 columns

\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": [ " id year tavg tmin tmax lat long name \\\n", "0 10160355000 1996 19.370000 15.860000 24.200000 36.93 6.95 Skikda \n", "1 10160355000 1997 18.790000 14.942857 21.828571 36.93 6.95 Skikda \n", "2 10160355000 1998 18.550000 14.730000 22.460000 36.93 6.95 Skikda \n", "3 10160355000 1999 19.536364 12.914286 20.600000 36.93 6.95 Skikda \n", "4 10160355000 2000 19.070000 16.111111 23.933333 36.93 6.95 Skikda \n", "\n", " country \n", "0 Algeria \n", "1 Algeria \n", "2 Algeria \n", "3 Algeria \n", "4 Algeria \n", "\n", "[5 rows x 9 columns]" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How Many Stations are There, Anyway?\n", "\n", "Just for fun, let's take a quick look at the stations we've got.\n", "\n", "This plot shows the number of stations with reported readings by year. That big drop in the 90's is interesting!\n", "\n", "To produce this, we just need to combine our dataframe with some simple plotting code from [matplotlib](http://matplotlib.org/). All we need for a simple plot like this are two lists: \"x\" and \"y\" values. For the \"x\" values, we simply enumerate all of the years with available data. This is done with `data['year'].min()` and `data['year'].max()` -- the brackets in this case select the `year` column, and `.min()` and `.max()` perform the expected aggregation functions.\n", "\n", "We then count the number of data points for each year, by calling `data[data.year == year].count()`. This selects all of the rows whose `year` matches the year in question, and then counts them. A simple query. :)\n", "\n", "The matplotlib code is fairly self-explanatory -- we just pass in the two lists, label the axes, and show the plot." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "years = xrange(data['year'].min(),data['year'].max())\n", "counts = [data[data.year == year].count() for year in years]\n", "\n", "plt.plot(years,counts,color='red')\n", "plt.xlabel('Year'); plt.ylabel('Stations Online')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEPCAYAAABlZDIgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1cVGXeP/DPQVBDodSNgWa0MRnEkcfUMVMTX4gPmKwP\n3RSVYOJd6Za6dlfW766w3Q162l2spZf3RkV2b2q3ie6m5GqisZuohGmSOe5iMcNAJpIoGE/X74/j\nDKIiT2fmzMDn/XrxOoczc875znCY71zXda7rkoQQAkRERN3kpXYARETUMzChEBGRIphQiIhIEUwo\nRESkCCYUIiJSBBMKEREpwukJpampCdHR0ZgzZw4AoKqqCnFxcQgJCcH06dNRXV3teG56ejoMBgNC\nQ0Oxc+dOx/aioiKEh4fDYDBgxYoVzg6ZiIi6wOkJJTMzE0ajEZIkAQAyMjIQFxeHEydOIDY2FhkZ\nGQCAkpISbNy4ESUlJcjLy8OyZctg7yKzdOlSZGdnw2w2w2w2Iy8vz9lhExFRJzk1oVgsFmzfvh1L\nlixxJIdt27YhJSUFAJCSkoLc3FwAwNatW5GUlAQfHx/o9XoEBwejsLAQNpsNNTU1MJlMAIDk5GTH\nPkRE5D6cmlB+/etf49VXX4WXV8tpKisrodFoAAAajQaVlZUAgPLycuh0OsfzdDodrFbrVdu1Wi2s\nVqszwyYioi5wWkL529/+hoCAAERHR6Ot0V0kSXJUhRERkWfzdtaB//nPf2Lbtm3Yvn07Ll68iHPn\nzmHhwoXQaDSoqKhAYGAgbDYbAgICAMglj7KyMsf+FosFOp0OWq0WFoul1XatVnvNcwYHB+Nf//qX\ns14SEVGPM2LECJw8eVKZgwkXyM/PF3fffbcQQognn3xSZGRkCCGESE9PF08//bQQQohjx46JyMhI\n8fPPP4t///vf4rbbbhPNzc1CCCFMJpPYv3+/aG5uFrNmzRI7duy45nlc9HI65YUXXlA7hKswpo5x\nx5iEcM+4GFPHuGNMSn5uOq2EciV71dbq1auRmJiI7Oxs6PV6bNq0CQBgNBqRmJgIo9EIb29vZGVl\nOfbJysrCokWLUFdXh/j4eMycOdNVYRMRUQe5JKFMmTIFU6ZMAQAMHjwYu3btuubznn32WTz77LNX\nbR8zZgyOHj3q1BiJiKh72FPeyWJiYtQO4SqMqWPcMSbAPeNiTB3jjjEpSbpUh9YjSJLU5h1lROTG\nMjMBjoKhCiU/N1lCISLXqasDoqIALy/glluAggJAkoCVK+Vlv35X7yNJLT+/+pXrY6YOY0IhIteo\nqwN8fYGvvgKEAGw2YPJk+bFVq4AbbgDq64HU1JZ95s+Xlz4+8jIrC9i40bVxU4exyouIXMPLS04k\nR44A4eFygpk6FdizR04mgFwKAeTnXfm7PSEBgMUCtNEfjTqHVV5E5FnuuENOChMmyMkEkJPI/v0t\nyQQA7rlHXr74IjBpkry+Zk3L848ckdcNBtfETZ3CEgoROd+VJY+OPNe+3tzc+nFvb6CpqWPHonax\nhEJEnqOuTl6GhXXs+cnJcvVYevrVyQQAtm2Tl9OmKRMfKYYlFCJyrjvvBL74AqitbV291R2dKfHQ\ndbGEQkSeY/9+ealUMgGAu+6Sl5IE3HijcselbmFCISLnEqLltl+l7N0LREbK6+fOyYmF8ySpjgmF\niJynoEBe/vd/K3/sw4flZPXUU/LvOh1waUpxUgfbUIjIOcLDga+/lted/X958mTLrcQhIcC33zr3\nfD0I21CIyL35+7ckk+Bg558vOFhOWn36ACdOyFVgnObC5ZhQiEh5NTUtPePNZtedt7ER+Phjef3T\nT4Ht2113bmJCISKF2cfieu89dc4/b558izIAzJ6tTgy9FNtQiEhZ7tJH5K23gGXL5PVx44ADB9SN\nx02xDYWI3FNVlbzU61UNAwCwdCnw6qvy+sGDcqJ76y11Y+rhWEIhIuUMHgycPat+6eRKL74IvPBC\ny+833wz88IN68bgRjyihXLx4EePHj0dUVBSMRiOeeeYZAEBaWhp0Oh2io6MRHR2NHTt2OPZJT0+H\nwWBAaGgodu7c6dheVFSE8PBwGAwGrOCsbkTu6+xZtSO4tuefl5PcmDHy76dPy3eEkaK8nXXg/v37\nY8+ePfD19UVjYyMmTZqEgoICSJKEVatWYdWqVa2eX1JSgo0bN6KkpARWqxXTpk2D2WyGJElYunQp\nsrOzYTKZEB8fj7y8PMzkLYFE7mX3bnl5773qxnE9hw7JyzvuAAoL5TvRrjUAJXWJU9tQfC9NhlNf\nX4+mpiYMGjQIAK5ZvNq6dSuSkpLg4+MDvV6P4OBgFBYWwmazoaamBiaTCQCQnJyM3NxcZ4ZNRF0R\nHy8vN2xQN46O2L8fmDhRLrVIErBli9oR9QhOTSjNzc2IioqCRqPB1KlTMXr0aADAG2+8gcjISKSm\npqK6uhoAUF5eDp1O59hXp9PBarVetV2r1cLKMXuI3E99vfyN31MUFADvvy+v26capm5x6l/fy8sL\nhw8fhsViwb59+5Cfn4+lS5eitLQUhw8fRlBQEJ544glnhkBErmDvQPjyy+rG0VkLFwKzZsnrJ0+q\nG0sP4LQ2lMvdeOONmD17Ng4dOoSYmBjH9iVLlmDOnDkA5JJHWVmZ4zGLxQKdTgetVguLxdJqu/Y6\nc0mnpaU51mNiYlqdj4icJClJXv7Xf6kbR1ds3y5Xe40fD5w5o3Y0Tpefn4/8/HznHFw4yenTp8XZ\ns2eFEELU1taKyZMni127dgmbzeZ4zu9//3uRlJQkhBDi2LFjIjIyUvz888/i3//+t7jttttEc3Oz\nEEIIk8kk9u/fL5qbm8WsWbPEjh07rnlOJ74cIroeuTVC7Si6ztPj7wYlPzedVkKx2WxISUlBc3Mz\nmpubsXDhQsTGxiI5ORmHDx+GJEkYPnw41q1bBwAwGo1ITEyE0WiEt7c3srKyIF3qcZuVlYVFixah\nrq4O8fHxvMOLyB1duunGIxmNQEmJPF2xkhOB9TLs2EhE3WMf4mTXLiA2Vu1ouqauDvD1BYYNA777\nTu1oXErJz00mFCLqHl9f+QPZ0//37GOQAfJw+K4cJVlFTChtYEIhUoG7DAapBKtVnvkR6BmvpwM8\nYugVIuoF7IMt3nqrunEoRasFoqLk9bo6dWPxQCyhEFHX9aTSiZ29PUWvB0pL1Y7G6VhCISL19e0r\nLz1hqJXOsN/ldeqUqmF4IiYUIuq8LVuAhgYgKMi9B4PsqsBAtSPwSKzyIqLO64lVXZerqgKGDJGn\nE7bPUd9DscqLiNSTkSEvx49XNw5nGjxYXnLa4E5hCYWIOqenl07sJAkYOBCoqVE7EqdiCYWIXGfS\npJbZDR94QF7ef7968bjSxYtqR+BRWEIhotaOHgUiIuT1Pn2ApiZ5XZJaSiW94f9MkuSfHj6jI0so\nROQ80dEt601Ncme/GTNaksiaNerEpYbekDgVxBIKEbWw393k7w/89FPrx0JDge+/B2pr1YnN1XpJ\nWxHH8moDEwpRN/XrJ0/lW1vLYdz79JGru3r4ZwqrvIjIOerr5R7wvT2ZAIC3Sya07VGYUIhIdued\n8vLYMXXjcBcDBqgdgcdhlRcRyXpJm0GHjRwJnDjR498PVnkRkbLsQ7UHBakbhzsZOVLtCDwOEwoR\nAaNHy8t//UvdONzJL3+pdgQex2kJ5eLFixg/fjyioqJgNBrxzDPPAACqqqoQFxeHkJAQTJ8+HdXV\n1Y590tPTYTAYEBoaip07dzq2FxUVITw8HAaDAStWrHBWyES9l33eDzbGt7CPBnDwoLpxeBCnJZT+\n/ftjz549OHz4MI4cOYI9e/agoKAAGRkZiIuLw4kTJxAbG4uMSwPNlZSUYOPGjSgpKUFeXh6WLVvm\nqNdbunQpsrOzYTabYTabkZeX56ywiXqvgAC1I3Av9uT63nuqhuFJnFrl5evrCwCor69HU1MTBg0a\nhG3btiElJQUAkJKSgtzcXADA1q1bkZSUBB8fH+j1egQHB6OwsBA2mw01NTUwmUwAgOTkZMc+RKSA\n3bvl5bvvqhuHu2IJpcOcmlCam5sRFRUFjUaDqVOnYvTo0aisrIRGowEAaDQaVFZWAgDKy8uh0+kc\n++p0Olit1qu2a7VaWK1WZ4ZN1LssXCgv4+PVjcNdlZWpHYHHcGrPHS8vLxw+fBg//fQTZsyYgT17\n9rR6XJIkSPZbFRWSlpbmWI+JiUFMTIyixyfqcWw2tSNwb1cOQePh8vPzkZ+f75Rju6Qr6I033ojZ\ns2ejqKgIGo0GFRUVCAwMhM1mQ8ClelutVouyy74JWCwW6HQ6aLVaWCyWVtu1Wm2b57o8oRBRB3nx\nhs821derHYGirvyivUbBwT6ddhX9+OOPjju46urq8Pe//x3R0dFISEhATk4OACAnJwdz584FACQk\nJGDDhg2or69HaWkpzGYzTCYTAgMD4e/vj8LCQgghsH79esc+RKQQo1HtCNxXDx++XklOK6HYbDak\npKSgubkZzc3NWLhwIWJjYxEdHY3ExERkZ2dDr9dj06ZNAACj0YjExEQYjUZ4e3sjKyvLUR2WlZWF\nRYsWoa6uDvHx8Zg5c6azwibqXezT+fLOybb18J7ySuLQK0S9mb+/PMUt/2+urRcMR8OhV4hIGT18\nvvRuY9tSp/DdIurt+vZVOwL35eOjdgQehQmFqLeyDwg5Y4a6cbgzDmHfKUwoRL1VcrK83LhR3Tjc\nWWCg2hF4FCYUot5q2zZ5yQEh2xYernYEHoUJhai36mEd9pxi3jy1I/AoTChEvdnAgWpH4N4SEuQl\nB4jsECYUot7IPsDqr3+tbhzuzl4duHevunF4CCYUot7o7rvl5YsvqhuHpyguVjsCj8CEQtQbHTmi\ndgSe5bvv1I7AIzChEPVGzc3sBd4Zp0+rHYFH4BVF1NucPCkvH3xQ3Tg8SQ+bE8VZODgkUW+j0QA/\n/NCjBzxUlCQBfn7AuXNqR+IUSn5uMqEQ9Ta9YARdRUkS0K8fcPGi2pE4BUcbJqLuGTxY7Qg8S1OT\n2hF4hHYTSnNzM9avX48XL91e+P333+PAgQNOD4yInOCRR+Slk+YU77E4a2OHtFvl9eijj8LLywuf\nffYZjh8/jqqqKkyfPh2HDh1yVYwdxiovonb4+ACNjazu6oweXkWo5Odmu1MAFxYWori4GNHR0QCA\nwYMHo6GhQZGTE5GLNTaqHQH1YO1WefXt2xdNl9Ufnj59Gl68f53Ic7H9hJyk3czw+OOPY968efjh\nhx/w7LPPYuLEiXjmmWc6dPCysjJMnToVo0ePRlhYGNauXQsASEtLg06nQ3R0NKKjo7Fjxw7HPunp\n6TAYDAgNDcXOnTsd24uKihAeHg6DwYAVK1Z09nUSUUGBvMzKUjcO6rE6dNvwN998g927dwMAYmNj\nMWrUqA4dvKKiAhUVFYiKisL58+cxZswY5ObmYtOmTfDz88OqVataPb+kpAT3338/Dh48CKvVimnT\npsFsNkOSJJhMJrz55pswmUyIj4/H8uXLMXPmzNYvhm0oRG0LCgIqKnpsW4DTeHnJ75mz37eoKOCL\nL1w+P41L21AAICQkBP7+/mhsbIQkSfj+++8xbNiwdvcLDAxE4KUZzwYOHIhRo0bBemmU02u9gK1b\ntyIpKQk+Pj7Q6/UIDg5GYWEhbr31VtTU1MBkMgEAkpOTkZube1VCIaLrqKhQOwLP5OXl/NuG/fyA\n8+eB0aOBf//buedyonarvN544w1oNBrExcXh7rvvxuzZszF79uxOn+jUqVMoLi7GHXfc4ThuZGQk\nUlNTUV1dDQAoLy+HTqdz7KPT6WC1Wq/artVqHYmJiDqhTx+1I/A8zn7PMjLkZAIApaXOPZeTtVtC\n+eMf/4hvv/0WQ4YM6fJJzp8/j3vuuQeZmZkYOHAgli5diueffx4A8Nxzz+GJJ55AdnZ2l49/ubS0\nNMd6TEwMYmJiFDkuUY8wcaLaEXgeb2/nzW5ZUADY26THjAGKiuS5arRa55wPQH5+PvKd1A+p3YQy\nbNgw+Pv7d/kEDQ0NWLBgAR588EHMnTsXABAQEOB4fMmSJZgzZw4AueRRVlbmeMxisUCn00Gr1cJi\nsbTarm3jDb88oRDRJQ88IC/z8tSNwxP17w/U1ip/3JkzgU8/ldc3bADuvVfu82I0OnUwyiu/aK9Z\ns0axY7ebUIYPH46pU6di9uzZ6Nu3LwC5EefKBvVrEUIgNTUVRqMRK1eudGy32WwICgoCAGzZsgXh\n4eEAgISEBNx///1YtWoVrFYrzGYzTCYTJEmCv78/CgsLYTKZsH79eixfvrxLL5ioV9q0SV66uMG3\nR/DzA6qqlD+uPZnU1rb+u5w719KZsk8f4Fe/AjIzlT+/E3SohDJs2DDU19ejvr4eQghI9hfbjn/8\n4x/44IMPEBER4egY+dJLL+HDDz/E4cOHIUkShg8fjnXr1gEAjEYjEhMTYTQa4e3tjaysLMe5srKy\nsGjRItTV1SE+Pp4N8kSdwQ6NXTdokPMm2Bo5snUySU4G3n+/5femJmDtWvmnb9+Wqrfly90yyXC0\nYaLeQJIAnQ64rEqZOig+HtixQ9nbhuvqAF9f4NVXgf/6r+s/d/58YMuWth9/+20gNbXLobhk+PoV\nK1YgMzPT0b5xZQDbtm1TJAAlMaEQXUNGhtzwe+QIcKl6mTrhiSeA3/9e2YTy/PPAb37T9WMePCgn\non375C8L3Ri80iUJ5dChQxg7dmybdwO4491TTChE1+DrK38j5v9G12zcCNx3n7Lvn14vV6N195j+\n/kBNzdXtMJ3ACbbawIRCdA09fLRcp6uqAoYMAc6cUW4ctH795PaQ7v5N7LENG9bldh6XJJTw6xSN\nJUnCkSNHFAlASUwoRNcgSfI3Wc6L3nWSJDeWL1yo3PEAZZJ8N4/lkqFX/vrXvypyAiJS0dCh8vJ/\n/kfdOHqCQ4eUSyhKuvLOMBWxyouop3r6aeCVV3r0fOguI0nArFnA9u3KHa9/f7ltq7vsd4x18W4v\nl84pv3nzZhgMBvj7+8PPzw9+fn7d6jlPRE5WUCAnkVdekX9nMlFGeXnX9z158upten3Xj3c5e2P8\npf58amo3oTz11FPYtm0bzp07h5qaGtTU1ODcuXOuiI2IOuvoUWDyZLnB18sLMJvVjqjnOHu28/vs\n3i2XRgyGlm32UonSo318+62yx+uCdhNKYGBgh+c/ISKVRUTIy9pauZd1cLC68fQkNTWde/6LLwLT\nprX8bq+O+u//lpdLlyoTl92FC8oerwvabUNZsWIFKioqMHfu3FZjec2fP98lAXYG21CoV7v1VuD7\n74F77gE++kjtaHoWSZLbKTr6ob17d0syqa2V97XPq2L/Oyn5WdWNO71c2g9l0aJFjpNe7t1331Uk\nACUxoVCvZO+4CMhDrTc0qBtPTyRJ8lhaP//c/nPtjeRAS98Vb285mQihXB+UK+MD3D+heBImFOpV\n7J3aAPkD7KuvWMXlLJIkj/zbkUE2Bw0CqquBXbuA2Fh5m334lgMHgEszzyqaULqRpFx2l9f27dtx\n1113YciQIRgyZAimTJmCTz75RJETE1EX1NXJgwV6ebUkk8WL5aoYJhPn6uh4WdXVcvKxJxMAeP11\neWlPJkr1uLdT+nhd1GbHxj//+c9Yt24dXnnlFYwZMwYAUFRUhNWrV8NiseCRRx5xWZBEvZa9KsPX\n9+pJngYNkmf34xwnrtGRb/G/+IW8LC5u+zkBAUBlpTIx2U2YcP0RiV2kzSqvUaNGoaCg4Kqpf8+c\nOYOJEyfi+PHjLgmwM1jlRT3K5W0jduPGAXv3Mom4WkfbKCSppfH9Slu2AH/7G6DQdOetHD0q3+Fn\nNne6pOqyKq9rzSM/ZMiQDk+wRURd9OKLcjIZPFj+ELP/HDjAZOKunn9eXrY1bNW8ec5JJkDLtAQv\nvOCc43dQm1Ve/v7+OHz4MKKiolpt/+qrr+Dn5+f0wIh6rdDQlk5qZ86oGwt13Esvycv4ePVi2LdP\nvXPjOgnl9ddfxy9/+Us89NBDGDNmDIQQKCoqwnvvvYcPPvjAlTES9R7220slibMruhMvr/Yb5Zua\n5Oep6YcfVD19m69+0qRJKCwsRFNTE9577z3k5OSgubkZhYWFmDx5sitjJOodVqyQP5SMRvnDS6tV\nOyKyay9RHDwoLx991PmxXI/KfZCc2g+lrKwMycnJ+OGHHyBJEh5++GEsX74cVVVVuPfee/Hdd99B\nr9dj06ZNuOmmmwAA6enpeOedd9CnTx+sXbsW06dPByDfYbZo0SJcvHgR8fHxyMzMvPrFsFGePBkn\nwnJfN9wgD7LZ1t/GPnOimn+7Ll4/Lh1tuDt8fHzwhz/8AceOHcP+/fvxpz/9Cd988w0yMjIQFxeH\nEydOIDY2FhkZGQCAkpISbNy4ESUlJcjLy8OyZcscL3Tp0qXIzs6G2WyG2WxGXl6eM0Mncr7LvxT1\n6SMv335bnVjo+nx8rv94Z8f5cga1q9vg5IQSGBjoaNQfOHAgRo0aBavVim3btiElJQUAkJKSgtzc\nXADA1q1bkZSUBB8fH+j1egQHB6OwsBA2mw01NTUwXeoUlJyc7NiHyOMkJMjfJleulJeSJFdxhYZ2\naT4LcgH7UCrXM3Cg8+O4ngED1D0/OplQmpqaujx0/alTp1BcXIzx48ejsrISGo0GAKDRaFB5qZNP\neXk5dDqdYx+dTger1XrVdq1WC6vV2qU4iFTj7y8nD/ttpW+/3fJB9fHHwDffqBcbXd+lKvnruvFG\n58dxPZePbKySNu/ysktKSsK6devQp08fjBs3Dj/99BNWrFiBp556qsMnOX/+PBYsWIDMzMyrbjmW\nJEnRfi1paWmO9ZiYGMTExCh2bKIu699fHliwTx/g3XdbppJlicQzaDTtzzcybJhrYmnLxx936Gn5\n+fnIz893SgjtJpSSkhL4+/vjf//3fzFr1ixkZGTg9ttv73BCaWhowIIFC7Bw4ULMnTsXgFwqqaio\nQGBgIGw2GwICAgDIJY+yy26VtFgs0Ol00Gq1sFgsrbZr27gD5vKEQuQWhgyRk4leD5SWqh0NdUV2\ntjxJ1q9+BfzpT9d+zrhxro2pi678or1mzRrFjt1ulVdjYyMaGhqQm5uLOXPmwMfHp8MlCiEEUlNT\nYTQasXLlSsf2hIQE5OTkAABycnIciSYhIQEbNmxAfX09SktLYTabYTKZEBgYCH9/fxQWFkIIgfXr\n1zv2IXJrb70ljwo8cCCTiSezD2fyP/9z9WP24XGWLHFdPO5KtCMzM1PccsstYubMmaKpqUmUlpaK\nSZMmtbebEEKIzz//XEiSJCIjI0VUVJSIiooSO3bsEGfOnBGxsbHCYDCIuLg4cfbsWcc+v/vd78SI\nESPEyJEjRV5enmP7oUOHRFhYmBgxYoR4/PHHr3m+DrwcIteyD5pCnq+tv+Xbb3v031jJz81O90MR\nQqCpqQne3u3Wlrkc+6GQW7H3TUhPB1avVjsa6q7Ro4GSkqv7ecyfLw/86KGfPS6dYOvixYvYvHkz\nTp06hcZLk8tIkoTn7QOhuREmFHILBw+2zHsxaJBc5UWezz4T46xZwPbtLdtHjQKOH2dCQQca5X/5\ny1/ipptuwpgxY9C/f39FTkrUo9TVybeV1te33s653XsW+yjPn37aervK42e5k3ZLKGFhYfj6669d\nFU+3sIRCqrDPWyJJ8pwUf/6zx9zxQ53k5dUylYCd/ZZwD/3scenQK3feeSeOHDmiyMmIepyCAjmZ\nDBwo93Y/fJjJpCcbOvTqbSoPyOhO2i2hjBo1CidPnsTw4cPRr18/eSdJcsskwxIKuRwHdOxdUlOB\nd95p/ff28GvApY3yp06dcpwUgOPEer1ekQCUxIRCLjN4MHD2rLweEwPs2aNqOOQiVVVyR9Vdu4DY\nWHkbE0rLsTpy2/Dhw4fx+eefQ5IkTJ48GZGRkYqcXGlMKOQSPj5AY6M8jMqnn7Z8sFDvIEnAXXcB\ne/e2/A4woaADbSiZmZl48MEHcfr0aVRWVuLBBx/E2rVrFTk5kcfx9paTSViYvGQy6Z2OHm39u336\ngV6u3RJKeHg49u/fjwGXhka+cOEC7rjjDhy98g11AyyhkFPZv4lGRQHFxerGQuqRJDmBXOqXB0kC\n+vWTJ+DyQC6fYMvrsolbvNxgEhcilxo9uiWZ3HMPkwnJUzVfzg3mInEH7XZsfOihhzB+/HjMnz8f\nQgjk5uZi8eLFroiNSF2X93gH5PlLONw8SdLV7SWBgerE4mY61ChfVFSEgoICR6N8dHS0K2LrNFZ5\nkWIyMoBnnpHXORYXXe7KjoyS5NGjIrjkLq9z587B398fVZfGIbI/zX778ODBgxUJQElMKKQI+5hN\nAHDmjHyLMJHd8OHAqVOtE8r777dMmuZhXDKWV1JSEj755BPcfvvt15z/pJRzO1BPM3SoPEJwSYn8\n++efM5nQ1RISAPudrgUF8vKee9SLx410evh6d8YSCnWZnx9w/nzL7+ysSG2xl2A//lju4JiV5bF9\nUAAXjzYcGxuL3bt3t7uNyONkZ8sjBO/fLyeTwYPlKi6i67GPOvzyy/L4beTQZkKpq6tDbW0tTp8+\n7WhHAeS2FavV6pLgiJzCapWrt678VsZkQp3xzTe8XfgKbSaUdevWITMzE+Xl5RgzZoxju5+fHx57\n7DGXBEekqKoqICioZd6SsDAgLw9YuVJuVCXqjPPnWzo3EoAOtKGsXbsWy5cvd1U83cI2FLou+80l\nvr7Ajz+2VF0QdZb9WrIPxePBnzsu7Sm/fPlyfP3119i0aRPef/99x09HLF68GBqNBuHh4Y5taWlp\n0Ol0iI6ORnR0NHbs2OF4LD09HQaDAaGhodi5c6dje1FREcLDw2EwGLBixYrOvD4imX2Eh1dfBS5c\nYDKh7rGP3cUSSmuiHS+88IKIiYkRN998s1i0aJHQaDRiwYIF7e0mhBBi37594ssvvxRhYWGObWlp\naeL111+/6rnHjh0TkZGRor6+XpSWlooRI0aI5uZmIYQQ48aNE4WFhUIIIWbNmiV27NhxzfN14OVQ\nbzR4sDwN0wlzAAAWAklEQVTH3oQJakdCPYXFYp+3Uf7xYEp+brZbQvm///s/7Nq1C0FBQXj33Xfx\n1Vdfobq6ukPJavLkyRg0aNC1kthV27Zu3YqkpCT4+PhAr9cjODgYhYWFsNlsqKmpgenSEBjJycnI\nzc3t0PmJEB8vt534+QH//Kfa0VBPodUCGzaoHYXbaTeh3HDDDejTpw+8vb3x008/ISAgAGVlZd06\n6RtvvIHIyEikpqY6klN5eTl0Op3jOTqdDlar9artWq2Wd5lRx2RmAvYq1XPn1I2Fep577wXMZsAN\nZ69VS7v9UMaOHYuzZ8/iP//zPzF27FgMGDAAd955Z5dPuHTpUjz//PMAgOeeew5PPPEEsrOzu3y8\nK6WlpTnWY2JiEBMTo9ixyYO89ZZ89xYA1NaqGwv1XMHBakfQafn5+cjPz3fKsdtNKG+99RYA4NFH\nH8WMGTNw7ty5bs3YGBAQ4FhfsmQJ5syZA0AueVxe8rFYLNDpdNBqtbBYLK22a7XaNo9/eUKhXqqq\nCli2TF6vrWUDPNFlrvyivWbNGsWO3W6VV+xlM9INHz4ckZGRrbZ1ls1mc6xv2bLFcQdYQkICNmzY\ngPr6epSWlsJsNsNkMiEwMBD+/v4oLCyEEALr16/H3Llzu3x+6gU0Gnl55gyTCZELObWnfFJSEvbu\n3Ysff/wRQ4cOxZo1a5Cfn4/Dhw9DkiQMHz4c69atAwAYjUYkJibCaDTC29sbWVlZjkEps7KysGjR\nItTV1SE+Ph4zZ87szmumnuzkSflWTj8/DuxI5GJtdmz84x//6Ogpf8sttzi2+/n54eGHH3bL3vLs\n2Ejo00ceX4lVXUQd4pL5UOzYU548xqBBQHW1PLxKebna0RB5BJcklIMHD0Kn0yEoKAgAkJOTg82b\nN0Ov1yMtLY0TbJF78faW5/keNEhulCeiDnHJ0CsPP/ww+vXrBwDYt28fVq9ejZSUFPj7++Phhx9W\n5OREiggNlZPJmDFMJkQqarNRvrm52VEK2bhxIx555BEsWLAACxYs6NZtw0SK2rgR+PZbue3k0CG1\noyHq1dosoTQ1NaGhoQEAsGvXLkydOtXxWCMHRCN3cd998pLXJJHqrjun/JQpU/CLX/wCvr6+mDx5\nMgDAbDbjpptuclmARG2y3z5u78RIRKq67l1eX3zxBSoqKjB9+nQMuDQz2YkTJ3D+/HncfvvtLguy\no9go38vY56Tg35yoy1w2p/yECROu2hYSEqLIiYm6ZdQoefn22+rGQUQO7fZD8SQsofQiLJ0QKcKl\nMzYSuZ1HHpGX6enqxkFErbCEQp6HpRMixbCEQr1XQYG8HDNG3TiI6CosoZDnqKsDfH3ldf6diRTB\nEgr1TvZkwju7iNwSEwp5Bnu7yapVQGqqurEQ0TWxyovcn4+PPLTKnDnAtm1qR0PUo7h0PhRPwoTS\nA3l5ye0lYWHA0aNqR0PU47ANhXoHSZKTycSJTCZEHoAJhdyTfdrpVatabhUmIrfm1ISyePFiaDQa\nhIeHO7ZVVVUhLi4OISEhmD59Oqqrqx2Ppaenw2AwIDQ0FDt37nRsLyoqQnh4OAwGA1asWOHMkMkd\nnDwJ2GzyHCevv652NETUQU5NKA899BDy8vJabcvIyEBcXBxOnDiB2NhYZGRkAABKSkqwceNGlJSU\nIC8vD8uWLXPU6y1duhTZ2dkwm80wm81XHZN6GINBXtbUqBsHEXWKUxPK5MmTMWjQoFbbtm3bhpSU\nFABASkoKcnNzAQBbt25FUlISfHx8oNfrERwcjMLCQthsNtTU1MBkMgEAkpOTHftQD2QfRXjcOOCG\nG9SNhYg6xeVtKJWVldBoNAAAjUaDyspKAEB5eTl0Op3jeTqdDlar9artWq0WVqvVtUGTa2RnA8eP\ny1VdBw6oHQ0RddJ150NxNkmSINk7rCkkLS3NsR4TE4OYmBhFj09OsGULsGBBy3AqnM6XyGny8/OR\nn5/vlGO7PKFoNBpUVFQgMDAQNpsNAQEBAOSSR1lZmeN5FosFOp0OWq0WFoul1XatVtvm8S9PKOQB\njh4F5s+X1wcPBi77WxOR8q78or1mzRrFju3yKq+EhATk5OQAAHJycjB37lzH9g0bNqC+vh6lpaUw\nm80wmUwIDAyEv78/CgsLIYTA+vXrHfuQh6urAyIi5HWzGThzhu0mRB7MqSWUpKQk7N27Fz/++COG\nDh2KF198EatXr0ZiYiKys7Oh1+uxadMmAIDRaERiYiKMRiO8vb2RlZXlqA7LysrCokWLUFdXh/j4\neMycOdOZYZOr2Ad73LABCA5WNxYi6jYOvULqGDwYOHsWuOce4KOP1I6GqNfi0Cvk2VJS5GTi58dk\nQtSDsIRCrnX0aEu7Cf9WRKrjaMNtYELxAPbbxM+ckau9iEhVrPIiz2TvBZ+ezmRC1AOxhEKuYy+d\n8G9E5DZYQiHPM3asvHz/fXXjICKnYQmFXIOlEyK3xBIKeRZ7B8a331Y3DiJyKiYUcq5HHpGHWAkI\nAFJT1Y6GiJyIVV7kXKzqInJrrPIiz5CQIC/ZEE/UK7CEQs7D0gmR22MJhdzf7t3ycuJEdeMgIpdh\nCYWcw9sbaGpi6YTIzbGEQu6vqQno21ftKIjIhZhQSHn224O3b1c3DiJyKVZ5kfLYGE/kMVjlRe7P\n3jueiHoNJhRSlr3vyf796sZBRC6nWkLR6/WIiIhAdHQ0TCYTAKCqqgpxcXEICQnB9OnTUV1d7Xh+\neno6DAYDQkNDsXPnTrXCpvb87W/yMjxc3TiIyOVUSyiSJCE/Px/FxcU4cOAAACAjIwNxcXE4ceIE\nYmNjkZGRAQAoKSnBxo0bUVJSgry8PCxbtgzNzc1qhU7XI0RLGwoR9SqqVnld2RC0bds2pKSkAABS\nUlKQm5sLANi6dSuSkpLg4+MDvV6P4OBgRxIiN6TXqx0BEalA1RLKtGnTMHbsWPz5z38GAFRWVkKj\n0QAANBoNKisrAQDl5eXQ6XSOfXU6HaxWq+uDpuvbskVecuwuol7JW60T/+Mf/0BQUBBOnz6NuLg4\nhIaGtnpckiRI16k6aeuxtLQ0x3pMTAxiYmKUCJc64le/kpeTJqkbBxG1KT8/H/n5+U45tmoJJSgo\nCABw8803Y968eThw4AA0Gg0qKioQGBgIm82GgIAAAIBWq0VZWZljX4vFAq1We83jXp5QyMUqKtSO\ngIjaceUX7TVr1ih2bFWqvGpra1FTUwMAuHDhAnbu3Inw8HAkJCQgJycHAJCTk4O5c+cCABISErBh\nwwbU19ejtLQUZrPZcWcYuRE2yBP1aqqUUCorKzFv3jwAQGNjIx544AFMnz4dY8eORWJiIrKzs6HX\n67Fp0yYAgNFoRGJiIoxGI7y9vZGVlXXd6jBSURslRyLq+Tj0Cilj925g2jRg1y4gNlbtaIiog5T8\n3GRCIWX84hfAmTMcv4vIw3AsL3IvGzfKycSLlxNRb8YSCnVPXV3LQJC1tcANN6gbDxF1Ckso5D5u\nvVVebtjAZELUy7GEQt3DuU+IPBpLKOQenn9eXirYMYqIPBdLKNR1LJ0QeTyWUEh9GzfKy5tvVjcO\nInIbLKFQ5913X0tC4Z1dRB6NJRRSzwMPMJkQ0TWpNtoweaD164G//EVeZ0mQiK7AhEIdEx0NHD4s\nr9fWqhsLEbklJhS6vst7wgOA2cxqLiK6JrahUNsuTybLlsnVXMHB6sZERG6LCYWutns34O3dkkye\new7405/UjYmI3B6rvEiWnQ1s3gz8/e9AY6O8rW9fYPt2zm9CRB3ChNIbnTwJRETIVVrXMmgQUFXl\n2piIyOOxyqs3yMwEbrkF6NNHHi7FYGidTAID5cZ2IeQfJhMi6gKPSih5eXkIDQ2FwWDAyy+/rHY4\n7m3KFDl5SBKwciVgswHNzfJjISEtyUMI+TE2thNRN3lMQmlqasJjjz2GvLw8lJSU4MMPP8Q333yj\ndljtys/Pd/5JsrOB+Hhg1Ch51kRJAvbtkx8LCAAsllYJJH/dOufH1EkueZ86yR1jAtwzLsbUMe4Y\nk5I8JqEcOHAAwcHB0Ov18PHxwX333YetW7eqHVa7FL+A6uqAp58GBgxoKYEsWQLs2AEcPy4nDW9v\nYNcueb2yEtBqnRuTAhhTx7ljXIypY9wxJiV5TKO81WrF0KFDHb/rdDoUFhaqGFEXWK3yh/6JE8Cp\nU8D58/IdVcePyz/V1UBDQ+eGNbn5ZuCTT4Bx45wWNhFRR3hMQpHsc2+0/0TnBtIV3Z2Ayl4SGTIE\nWLgQ+O1v2VudiNyP8BBffPGFmDFjhuP3l156SWRkZLR6zogRIwQA/vCHP/zhTwd/RowYodjntMfM\nh9LY2IiRI0di9+7duOWWW2AymfDhhx9i1KhRaodGRETwoCovb29vvPnmm5gxYwaampqQmprKZEJE\n5EY8poRCRETuza1vG168eDE0Gg3Cw8Md2w4cOACTyYTo6GiMGzcOBw8ebLXP999/j4EDB+L11193\nbCsqKkJ4eDgMBgNWrFjh0piOHDmCCRMmICwsDBEREaivr1c1posXLyIpKQkREREwGo3IyMhw7OPs\nmL766itMmDABERERSEhIQE1NjeOx9PR0GAwGhIaGYufOnarH9Pe//x1jx45FREQExo4diz179jgl\nps7GZafGdX69mNS6ztuKyVXXeVlZGaZOnYrRo0cjLCwMa9euBQBUVVUhLi4OISEhmD59Oqqrqx37\nOPta72xMil7rirXGOMG+ffvEl19+KcLCwhzbpkyZIvLy8oQQQmzfvl3ExMS02mfBggUiMTFRvPba\na45t48aNE4WFhUIIIWbNmiV27NjhkpgaGhpERESEOHLkiBBCiKqqKtHU1KRqTO+++6647777hBBC\n1NbWCr1eL7777juXxDR27Fixb98+IYQQ77zzjnjuueeEEEIcO3ZMREZGivr6elFaWipGjBghmpub\nVY2puLhY2Gw2IYQQX3/9tdBqtY59lIyps3HZqXGdtxWTmtd5WzG56jq32WyiuLhYCCFETU2NCAkJ\nESUlJeLJJ58UL7/8shBCiIyMDPH0008LIVxzrXc2JiWvdbcuoUyePBmDBg1qtS0oKAg//fQTAKC6\nuhrayzrt5ebm4rbbboPRaHRss9lsqKmpgclkAgAkJycjNzfXJTHt3LkTERERjm9UgwYNgpeXl6ox\nBQUF4cKFC2hqasKFCxfQt29f+Pv7uyQms9mMyZMnAwCmTZuGzZs3AwC2bt2KpKQk+Pj4QK/XIzg4\nGIWFharGFBUVhcDAQACA0WhEXV0dGhoaFI+ps3EB6l3nbcWk5nXeVkyuus4DAwMRFRUFABg4cCBG\njRoFq9WKbdu2ISUlBQCQkpLiOIcrrvXOxqTkte7WCeVaMjIy8MQTT2DYsGF48skn8dJLLwEAzp8/\nj1deeQVpaWmtnm+1WqHT6Ry/a7VaWK1Wp8aUnp4OQL7YJUnCzJkzMWbMGLz66quqxWR/n2bMmAF/\nf38EBQVBr9fjySefxE033eSSmEaPHu0Y3eCjjz5CWVkZAKC8vLzVuXU6HaxW61XbXRnT5TZv3owx\nY8bAx8fHJe/T9eJS8zpvK6YTJ06odp23FZMa1/mpU6dQXFyM8ePHo7KyEhqNBgCg0WhQWVkJwPXX\nekdiulx3r3WPSyipqalYu3Ytvv/+e/zhD39AamoqACAtLQ2//vWv4evrC+Hi+wyujGnx4sUAgIaG\nBhQUFOAvf/kLCgoKsGXLFnz22Wcd76SpYEz29+mDDz5AXV0dbDYbSktL8dprr6G0tNTp8QDAO++8\ng6ysLIwdOxbnz59H3759XXLe7sR07NgxrF69GutcPP5ZW3GpeZ23FVNjY6Nq13lbMbn6Oj9//jwW\nLFiAzMxM+Pn5tXpMkiSXvBfdjUmJa91jbhu2O3DgAHbt2gUAuOeee7BkyRLH9s2bN+Opp55CdXU1\nvLy8cMMNN2D+/PmwWCyO/S0WS6tqMmfGNHToUNx1110YPHgwACA+Ph5ffvklHnzwQdVi+uc//4l5\n8+ahT58+uPnmmzFx4kQUFRVh0qRJTo9p5MiR+PTTTwHI32o/+eQTAPI3n8tLBhaLBTqdDlqtVrWY\n7OebP38+1q9fj+HDhztidXZM14pr+/btANS9ztt6r9S8ztt6n1x5nTc0NGDBggVYuHAh5s6dC0Au\nAVRUVCAwMBA2mw0BAQEAXHetdyYm+/mUuNY9roQSHByMvXv3AgA+++wzhISEAAD27duH0tJSlJaW\nYuXKlfh//+//YdmyZQgMDIS/vz8KCwshhMD69esdb7CzY5o+fTqOHj2Kuro6NDY2Yu/evRg9erSq\nMYWGhuKzzz4DAFy4cAH79+9HaGioS2I6ffo0AKC5uRm//e1vsXTpUgBAQkICNmzYgPr6epSWlsJs\nNsNkMqkaU3V1NWbPno2XX34ZEyZMcDw/KCjI6TFdK65HH30UgLrXeVvv1YwZM1S7ztt6n1x1nQsh\nkJqaCqPRiJUrVzq2JyQkICcnBwCQk5PjOIcrrvXOxqTotd6l2whc5L777hNBQUHCx8dH6HQ68c47\n74iDBw8Kk8kkIiMjxR133CG+/PLLq/ZLS0sTr7/+uuP3Q4cOibCwMDFixAjx+OOPuzSmDz74QIwe\nPVqEhYU57qpQM6aLFy+KBx54QISFhQmj0djqLiFnxpSdnS0yMzNFSEiICAkJEc8880yr5//ud78T\nI0aMECNHjnTcnaZmTL/5zW/EgAEDRFRUlOPn9OnTisfU2bgu58rrvL2Y1LjOrxeTq67zzz//XEiS\nJCIjIx3XyY4dO8SZM2dEbGysMBgMIi4uTpw9e9axj7Ov9c7GpOS1zo6NRESkCI+r8iIiIvfEhEJE\nRIpgQiEiIkUwoRARkSKYUIiISBFMKEREpAgmFKIuEkJg8uTJyMvLc2z76KOPMGvWLBWjIlIP+6EQ\ndcOxY8fwH//xHyguLkZDQwNuv/12fPrpp47hKzqjsbER3t4eNxoSkQMTClE3Pf300/D19cWFCxcw\ncOBAfPfdd/j666/R0NCAtLQ0JCQk4NSpU0hOTsaFCxcAAG+++SYmTJiA/Px8PPfccxg8eDCOHz+O\nb7/9VuVXQ9R1TChE3VRbW4vbb78dffv2xd13343Ro0fjgQceQHV1NcaPH4/i4mJIkgQvLy/069cP\nZrMZ999/Pw4ePIj8/HzcfffdOHbsGG699Va1XwpRt7B8TdRNvr6+uPfeezFw4EBs2rQJf/3rX/Ha\na68BAH7++WeUlZUhMDAQjz32GL766iv06dMHZrPZsb/JZGIyoR6BCYVIAV5eXvDy8oIQAh9//DEM\nBkOrx9PS0hAUFIT169ejqakJ/fv3dzw2YMAAV4dL5BS8y4tIQTNmzMDatWsdvxcXFwMAzp0755hm\n9f3330dTU5Mq8RE5ExMKkUIkScJzzz2HhoYGREREICwsDC+88AIAYNmyZcjJyUFUVBS+/fZbDBw4\nsNV+RD0BG+WJiEgRLKEQEZEimFCIiEgRTChERKQIJhQiIlIEEwoRESmCCYWIiBTBhEJERIpgQiEi\nIkX8f/FMQqBFR/RKAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding Averages\n", "\n", "The next section finds the average temperature for each station for the periods of 1950-1954 and 2008-2012. Most of these operations we've seen before, but we do introduce two new ones -- queries with multiple selectors and `groupby`.\n", "\n", "```python\n", "early_years = data[(data.year >= 1950) & (data.year < 1955)][['id','tavg']]\n", "```\n", "This line uses simple ampersand syntax (`&`) to show that we want to make *two* selections on the data -- the year should be greater than or equal to 1950, and less than 1955. We also select only the `id` and `tavg` columns from the result.\n", "\n", "```python\n", "early_years = early_years.groupby('id').mean()\n", "```\n", "We then use the powerful `groupby` syntax to find the mean of all the years in the collection that have the same station ID. This gives us the average temperature over the whole range for each station.\n", "\n", "Finally, we use a simple subtraction operator to find the difference in averages. `dropna` is called to remove stations that do not have data available in both ranges." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from datetime import datetime\n", "\n", "data = data[(data.year >= 1950) & (data.year <= 2013)]\n", "data['date'] = data['year'].apply(lambda x: datetime(x,1,1))\n", "\n", "early_years = data[(data.year >= 1950) & (data.year < 1955)][['id','tavg']]\n", "early_years = early_years.groupby('id').mean()\n", "\n", "late_years = data[(data.year >= 2008) & (data.year < 2013)][['id','tavg']]\n", "late_years = late_years.groupby('id').mean()\n", "\n", "deltas = (late_years['tavg']-early_years['tavg']).dropna()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We combine the evaluated deltas with the index using the same square brackets to assign a column.\n", "\n", "The last line groups by country, and counts the occurances of each. 1219 of the 1529 stations represented are in the USA." ] }, { "cell_type": "code", "collapsed": false, "input": [ "change = index.copy().set_index('id')\n", "change['delta'] = deltas\n", "change = change.dropna()\n", "print \"Total stations represented:\",len(change)\n", "change.groupby('country').count().sort('country',ascending=False)[['country']]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Total stations represented: 1686\n" ] }, { "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", " \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", " \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", "
country
country
United States Of America 1222
Russian Federation (Asian Sector) 71
China 52
Australia 49
Japan 46
Canada 34
Russian Federation (European Sector) 19
Pakistan 17
Sudan 15
Italy 15
Philippines 14
Thailand 13
Malaysia 11
Ukraine 10
Poland 10
South Africa 9
Republic Of Korea 9
Algeria 6
Turkmenistan 6
Uzbekistan 6
Kazakhstan 5
Zambia 3
Belarus 3
Egypt 3
Federated States Of Micronesia 3
Bangladesh 3
Zimbabwe 2
Libya 2
Lithuania 2
Tanzania 2
Nigeria 2
Cocos Islands (Australia) 1
Antarctica 1
Armenia 1
Wake Island (U.S.A.) 1
Austria 1
Saudi Arabia 1
Belau 1
Bermuda (U.K.) 1
Portugal 1
Northern Mariana Islands (U.S.A.) 1
Latvia 1
Kyrgyzstan 1
Cuba 1
Czech Republic 1
Norfolk Island (Australia) 1
Estonia 1
Ghana 1
Iceland 1
Indonesia 1
Moldova 1
Marshall Islands 1
Coral Sea Islands (Australia) 1
\n", "

53 rows \u00d7 1 columns

\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ " country\n", "country \n", "United States Of America 1222\n", "Russian Federation (Asian Sector) 71\n", "China 52\n", "Australia 49\n", "Japan 46\n", "Canada 34\n", "Russian Federation (European Sector) 19\n", "Pakistan 17\n", "Sudan 15\n", "Italy 15\n", "Philippines 14\n", "Thailand 13\n", "Malaysia 11\n", "Ukraine 10\n", "Poland 10\n", "South Africa 9\n", "Republic Of Korea 9\n", "Algeria 6\n", "Turkmenistan 6\n", "Uzbekistan 6\n", "Kazakhstan 5\n", "Zambia 3\n", "Belarus 3\n", "Egypt 3\n", "Federated States Of Micronesia 3\n", "Bangladesh 3\n", "Zimbabwe 2\n", "Libya 2\n", "Lithuania 2\n", "Tanzania 2\n", "Nigeria 2\n", "Cocos Islands (Australia) 1\n", "Antarctica 1\n", "Armenia 1\n", "Wake Island (U.S.A.) 1\n", "Austria 1\n", "Saudi Arabia 1\n", "Belau 1\n", "Bermuda (U.K.) 1\n", "Portugal 1\n", "Northern Mariana Islands (U.S.A.) 1\n", "Latvia 1\n", "Kyrgyzstan 1\n", "Cuba 1\n", "Czech Republic 1\n", "Norfolk Island (Australia) 1\n", "Estonia 1\n", "Ghana 1\n", "Iceland 1\n", "Indonesia 1\n", "Moldova 1\n", "Marshall Islands 1\n", "Coral Sea Islands (Australia) 1\n", "\n", "[53 rows x 1 columns]" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "`describe` gives you standard dataset summary data." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print change['delta'].describe()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "count 1686.000000\n", "mean 0.570878\n", "std 1.711249\n", "min -14.261167\n", "25% 0.051986\n", "50% 0.661217\n", "75% 1.250096\n", "max 9.831667\n", "Name: delta, dtype: float64\n" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "matplotlib is great for quick analysis and visualizations. Let's plot a histogram of our temperature deltas." ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.hist(list(change['delta']),50)\n", "plt.xlim([-5,5])\n", "plt.title('Average Temperature Change, 1950-2010')\n", "plt.xlabel('Delta Temperature (C)')\n", "plt.ylabel('Count of Stations')\n", "plt.savefig('viz/temp_change.png')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEZCAYAAACaWyIJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVGXfP/DPIKghICKxyKCjouIosnmjmelwI6gZ5haK\nG6bmnabPrzRzuTOhRelxuzX1VgsNM/dScM11LM0e3Jew0ARlFxcEBGS7fn+QJ0eWIWBmwPm8Xy9e\nzpz1O8c55zvXcq4jE0IIEBGR0TIxdABERGRYTAREREaOiYCIyMgxERARGTkmAiIiI8dEQERk5JgI\niOoYlUqFiIgIQ4dBRoSJQA9UKhVsbGxQUFBg6FBq5O2334alpSUsLS3RqFEjNGzYUHo/YMAAQ4dX\nLaGhoRgzZoxe91lQUIDQ0FC0b98eFhYWaN26NSZMmIBbt24BAGQyGWQymV5jqqm0tDQMHDgQTk5O\nMDExwe3btzXmJycn4/XXX0fz5s3h7OyMtWvXasw3MTGBhYWF9H2aNGmSxvxly5bB0dERTZs2xYQJ\nEyo9lyIjI9G1a1c0bdoUzs7OmDVrFoqLi6X59+/fx+DBg2FhYQGFQoEtW7ZI8woLCzFs2DC0bt0a\nJiYmOHHiRJntz5o1C7a2trC1tcXs2bP/1nGqq5gIdCwhIQExMTGws7NDdHR0rW+/qKio1rdZkTVr\n1iA7OxvZ2dmYO3cuRowYIb3ft2+f3uKoKn0cm+rsY9iwYdi7dy+2bNmCrKwsXLp0CV27dsWxY8d0\nEKF+mJiY4NVXX8V3331X7vzRo0ejbdu2uHPnDvbt24e5c+dCrVZrLHPlyhXp+7Ru3Tpp+g8//IDP\nP/8cx44dw61bt3Dz5k3Mnz+/wljy8vKwfPly3Lt3D//3f/+Ho0ePYvHixdL8d955B40bN8adO3fw\n7bffYvLkyYiNjZXm9+rVC5s2bYKDg0OZhLx27VpERUXh8uXLuHz5Mvbs2VMmqdVLgnQqLCxMBAYG\nik8//VS89tprQggh8vPzRdOmTcXVq1el5e7cuSNeeOEFkZGRIYQQYs+ePcLd3V1YW1uLHj16iMuX\nL0vLtmrVSnz++efCzc1NNG7cWBQVFYmFCxeKtm3bCktLS6FUKsWuXbuk5YuLi8X06dOFra2taN26\ntfjiiy+ETCYTxcXFQgghMjMzxfjx44Wjo6NwcnISH374oTSvIvPnzxejR4+W3p8+fVq89NJLwtra\nWri7uwu1Wi3N6927t/jwww9Fjx49hIWFhQgMDBQZGRli5MiRwsrKSvzjH/8QCQkJ0vIymUysWLFC\ntGnTRtja2oqZM2eKkpISaX5ERITo2LGjaNasmejbt6+4deuWxrqrVq0SLi4uok2bNkIIIf7nf/5H\nODs7CysrK+Ht7S1++uknIYQQBw4cEA0bNhRmZmbCwsJCeHh4SMf3yJEj5X7W+Ph4IZPJREREhGjZ\nsqXo3bu31piedvjwYfHCCy+IpKSkCo+tSqUS8+bNEy+//LKwtLQUAQEB4u7du9L8YcOGCQcHB9G0\naVPRq1cv8euvv0rzQkJCxJQpU8SAAQOEpaWl6Natm/jjjz+k+T/88INo3769aNq0qZgyZYro1auX\n+Oqrr6p0bKuisLBQyGQyjfWys7OFTCaTvttCCDFp0iQxZswY6b1MJhM3btwod5vBwcHi3//+t/T+\n2LFjwsHBocoxLV26VAQGBgohhMjJyRENGzYU169fl+aPHTtWzJ49u8x6crlcnDhxQmPaSy+9JL78\n8kvp/fr160X37t2rHEtdxUSgY23bthWbNm0ScXFxwszMTNy5c0cIIcT48eM1vtwrV64U/fv3F0II\ncf78eWFnZydiYmJESUmJiIyMFAqFQhQUFAghSi9Unp6eIikpSeTn5wshhNixY4dITU0VQgixbds2\n0aRJE5GWliaEEOK///2vUCqVIjk5WTx48ED4+fkJExMT6WI/aNAg8fbbb4vc3Fxx584d4ePjI9au\nXVvp53r64piUlCSaN28uDhw4IIQovdg1b95cunj17t1btGvXTty8eVM8fPhQKJVK4eLiIo4ePSqK\niorE2LFjxZtvviltWyaTiX/+85/iwYMH4vbt26J9+/bSxWr37t3CxcVF/Pbbb6K4uFh8+umnokeP\nHhrrBgQEiAcPHkjHZtOmTeL+/fuiuLhYLFmyRDg4OIjHjx8LIYQIDQ3VuCAJIYRCoRBHjx6V3oeG\nhpZJBCEhISI3N1fk5eVpjelps2bNEiqVqtJj27t3b9G2bVtx/fp1kZeXJ1QqlcaFasOGDSInJ0cU\nFBSId999V0pgQpQmgubNm4szZ86IoqIiMWrUKDFixAghhBAZGRnCyspK7Nq1SxQXF4vly5cLMzMz\nERERUaVjWxXlJYKsrCwhk8mk774QQkycOFF4enpK72UymWjRooVwcHAQQ4YM0fhh4O7uLrZv3y69\nv3v3rpDJZOL+/ftViun1118Xc+bMEUKUnlvm5uYa85csWSIliqeVlwiaNm0qYmJipPdnz54VlpaW\nVYqjLmMi0KGffvpJNG7cWGRlZQkhSr/Qy5YtE0IIceTIEdG2bVtp2R49eohvvvlGCCHE22+/LebN\nm6exrQ4dOogff/xRCFF6odqwYUOl+/bw8BDR0dFCCCF8fX3FunXrpHlHjhyRSgRpaWmiUaNGIi8v\nT5q/efNm4evrW+n2n04E4eHhZS6mffv2FZGRkUKI0l+4CxYskObNmDFDvPrqq9L7PXv2aFzMZDKZ\n+OGHH6T3q1evFn5+fkIIIfr16ydduIQoLe2Ym5uL27dvS+seP3680tibNWsmlbCeLdkIUTYRlFci\niI+Pl+Zri+lpEydOlC7MFVGpVOKzzz6T3q9evVr069ev3GUfPHggZDKZ9B0bN26ceOutt6T5+/fv\nF66urkIIISIjI8tc2J2dnaXY/87nqEh5iUAIIXr27CmmTZsm8vPzxblz54SNjY0UlxCl50phYaHI\nzMwUU6dOFZ07d5Z+qLRt21bj+1BQUFDuPsoTEREhnJ2dxb1794QQQvz4449lShPr1q0rNzmXlwga\nNGggfv/9d+l9XFyckMlkWuOo69hGoEORkZEICAiApaUlAOCNN95AZGQkgNIG5NzcXMTExCAhIQGX\nLl3C4MGDAQC3bt3CkiVL0KxZM+kvKSkJKSkp0radnZ019rVx40Z4enpKy1+9ehV3794FAKSmpmos\nL5fLpde3bt1CYWEhHB0dpXXffvttZGRkVPlz3rp1Czt27NCI99SpU0hLS5OWsbe3l143btwYdnZ2\nGu9zcnI0tvl0vC1btpQ++61bt/D//t//k/bTvHlzAKWNkRUdm8WLF0OpVMLa2hrNmjXDw4cPpWNT\nXU/voyoxPWFra4vU1FSt23dwcJBev/DCC9LxKS4uxuzZs+Hi4oKmTZuidevWAKDxeZ4+1k+vm5KS\novF/D5T9LlT1c/xd3377LeLj4+Hs7Ix33nkHo0ePhpOTkzS/Z8+eMDU1RdOmTbF8+XIkJCTg2rVr\nAAALCwtkZWVJyz58+BAAYGlpiW+//bbCDgu7d+/G3LlzceDAAdjY2JS7rSfbe3KOalNeLBYWFn/j\nSNRNpoYO4HmVl5eH7du3o6SkBI6OjgCAx48fIzMzE5cvX0aXLl0QFBSELVu2wM7ODoGBgWjSpAmA\n0gvfv//9b8ydO7fC7T/diHXr1i1MmjQJx44dw0svvQSZTAZPT0+IPweWdXR0RGJiorT806+dnZ3R\nqFEj3Lt3DyYm1ftd0LJlS4wZM0ajga8yVekRc/v2bXTs2FF6/eSi0bJlS8ybNw/BwcFV2v5PP/2E\nRYsW4dixY+jUqRMAwMbGRjo25cXSpEkTPHr0SHr/dEIrbx9ViemJPn36YPny5UhOTta4EFbV5s2b\nER0djaNHj6JVq1bIzMzU+DyVadGiBfbs2SO9F0IgKSmpWp/j72rZsqXGvkeOHIlu3bqVu+yTz/Lk\n306dOuHixYsYNmwYAODSpUuwt7dHs2bNMGrUKIwaNarMNg4ePIhJkyZh//790v87ALRv3x5FRUW4\nceMGXFxcpO117ty5Sp/jSSxdu3b92+vWZSwR6Mju3bthamqKa9eu4dKlS7h06RKuXbuGV155BRs3\nbgRQejJs3boVmzdvxsiRI6V133rrLaxZswYxMTEQQuDRo0fYt29fmV/NTzx69AgymQy2trYoKSnB\nhg0bcPXqVWl+UFAQli9fjpSUFGRmZuLzzz+XLmSOjo4ICAjA9OnTkZ2djZKSEvzxxx/48ccfq/xZ\nR48ejT179uDQoUMoLi5Gfn4+1Gq1xi/Jpy9UVbloLV68GJmZmUhMTMSKFSswfPhwAKVdWBcsWCD1\n8nj48CF27NhR4Xays7NhamoKW1tbFBQU4OOPP9b4Refg4ICEhASNmDw8PLB161YUFRXh7Nmz+O67\n7ypNXn8nJj8/P/j7+2Pw4ME4f/48ioqKkJ2djTVr1mDDhg1aj1FOTg4aNWoEGxsbPHr0qMyPhcqO\n7auvvoorV64gKioKRUVFWLVqlUaS0/Y5VCoVwsLCKtx+fn4+8vPzy7wGgN9++w3Z2dkoKCjApk2b\ncPjwYUyfPh0AEBsbi4sXL6K4uBg5OTmYPn065HK59ENg7NixiIiIwLVr1/DgwQN88sknePPNNyuM\n49ixYxg1ahS+//576YL9RJMmTTBkyBB89NFHyM3NxcmTJ7Fnzx6NLsSPHz+WYn/69ZNYli5dipSU\nFCQnJ2Pp0qUYN25chbHUF0wEOrJx40aMHz8ecrkcdnZ2sLOzg729PaZOnYrNmzejpKQEPj4+sLCw\nQGpqKvr37y+t6+3tjS+//BJTp06FjY0N2rVrh40bN1Z4MVIqlZgxYwZeeuklODg44OrVq+jZs6c0\n/6233kJAQAC6dOkCb29vDBgwAA0aNJBKABs3bkRBQQGUSiVsbGzwxhtvlPsr+GlP93WXy+WIiorC\nggULYGdnh5YtW2LJkiUaF6WnYy+vn/yz719//XV4e3vD09MTr732GsaPHw8AGDRoEGbNmoURI0ag\nadOmcHNzww8//FDhdvr164d+/fqhffv2UCgUeOGFF9CyZUtp/htvvAEAaN68uXTR+OSTT/DHH3+g\nWbNmCA0NLfOL89l9aIvpWTt37sSrr76K4cOHw9raGm5ubjh//jz8/f21Hq+xY8eiVatWcHJyQufO\nnaUSYFWOra2tLXbs2IEPPvgAtra2uHbtGrp27YpGjRpV6XMkJSVpfK+eZW5uDisrK8hkMri6ukol\nXKC0C2jbtm1hY2ODdevW4YcffpCqntLT06V9tm3bFomJidi7dy8aNGgAAOjbty8++OAD+Pr6QqFQ\noG3btpUmpE8//RTZ2dno379/udVGq1evRl5eHuzs7DB69GisWbNGSjoA0KFDB5ibmyMlJQV9+/ZF\nkyZNpPsi/vWvfyEwMBBubm7o0qULAgMDy9zzUC/puhGiqKhIeHh4SF0n7927J/r06SPatWsn/P39\nxYMHD6RlFyxYIFxcXESHDh00Goeodu3fv1+0atXK0GFUSCaTaXR5JN0oLi4WLVq00OjqW5HExETx\n8ssv6yEqMgSdlwiWL18OpVIp/SoJDw+Hv78/4uLi4Ofnh/DwcAClxcNt27YhNjYWBw8exJQpU1BS\nUqLr8IxCfn4+9u/fj6KiIiQnJyMsLAxDhgwxdFhkAIcOHUJmZiYeP36MBQsWAAC6d++udT25XI6T\nJ0/qOjwyEJ0mgqSkJOzfvx8TJ06Uqgmio6MREhICAAgJCcHu3bsBAFFRUQgODoaZmRkUCgVcXFwQ\nExOjy/CMhhACoaGhsLGxgZeXFzp16oSPP/7Y0GFVqL4Nr1CfnD59Gi4uLnjxxRexb98+7N69W6oa\nIuOl015D7733HhYtWqTROJeeni51b7O3t0d6ejqA0q5tT/8ykcvltdJtjUq7ENanpPr0uDBUu+bP\nn1/p8AxknHRWIti7dy/s7Ow0ujE+S9vgWvxlSESkezorEfz888+Ijo7G/v37kZ+fj6ysLIwZMwb2\n9vZIS0uDg4MDUlNTpRuLnJycNPq3JyUlldvP2sPDA5cuXdJV2EREzyV3d3dcvHix/Jn6aJFWq9VS\nr6GZM2eK8PBwIYQQCxcuFLNmzRJCCPHrr78Kd3d38fjxY3Hz5k3Rpk0bjYHGntBTyFrNnz/f0CHU\nGTwWf+GxKMXj8Je6ciwqu3bq7c7iJ9U8s2fPRlBQECIiIqBQKLB9+3YApX3hg4KCoFQqYWpqitWr\nV7NqiIhID/SSCHr37o3evXsDKL29/8iRI+UuN3fu3EqHVSAiotrHO4urSaVSGTqEOoPH4i88FqV4\nHP5SH46F7M+6o3pDJpNVaawaIiL6S2XXTpYIiIiMHBMBEZGRYyIgIjJyTAREREaOiYCIyMgxERAR\nGTkmAiIiI8dEQERk5JgIiIiMHBMBEZGRYyIgIjJyTAREREaOiYCIyMgxERARGTkmAiIiI8dEQERk\n5JgIiIiMnM4SQX5+Prp16wYPDw8olUrMmTMHABAaGgq5XA5PT094enriwIED0joLFy5Eu3bt4Orq\nikOHDukqNCIieopOH1WZm5sLc3NzFBUVoWfPnli8eDGOHj0KS0tLTJ8+XWPZ2NhYjBw5EmfOnEFy\ncjL69OmDuLg4mJho5io+qpKMkZWVDbKzH1R7fUvLZsjKul+LEVF9Y7BHVZqbmwMACgoKUFxcjGbN\nmgFAucFERUUhODgYZmZmUCgUcHFxQUxMjC7DI6o3SpOAqPZfTZIIPf90mghKSkrg4eEBe3t7+Pr6\nolOnTgCAL774Au7u7pgwYQIyMzMBACkpKZDL5dK6crkcycnJugyPiIig40RgYmKCixcvIikpCT/+\n+CPUajUmT56M+Ph4XLx4EY6OjpgxY0aF68tkMl2GR0REAEz1sZOmTZtiwIABOHv2LFQqlTR94sSJ\nCAwMBAA4OTkhMTFRmpeUlAQnJ6dytxcaGiq9VqlUGtskIiJArVZDrVZXaVmdNRbfvXsXpqamsLa2\nRl5eHvr27Yv58+ejU6dOcHBwAAAsW7YMZ86cwebNm6XG4piYGKmx+MaNG2VKBWwsJmNUeh7U5HvP\n88bYVXbt1FmJIDU1FSEhISgpKUFJSQnGjBkDPz8/jB07FhcvXoRMJkPr1q2xdu1aAIBSqURQUBCU\nSiVMTU2xevVqVg0REemBTruP6gJLBGSMWCKgmjJY91EiIqr79NJYTGTsanpDGJEusWqISA9qo2qH\nVUNUE6waIiKiCjEREBEZOSYCIiIjx0RARGTkmAiIiIwcEwERkZFjIiAiMnJMBERERo6JgIjIyDER\nEBEZOSYCIiIjx0RARGTkmAiIiIwcEwERkZFjIiAiMnJMBERERk5niSA/Px/dunWDh4cHlEol5syZ\nAwC4f/8+/P390b59ewQEBCAzM1NaZ+HChWjXrh1cXV1x6NAhXYVGRERP0ekTynJzc2Fubo6ioiL0\n7NkTixcvRnR0NGxtbfHBBx/g888/x4MHDxAeHo7Y2FiMHDkSZ86cQXJyMvr06YO4uDiYmGjmKj6h\njOojPqGMDM1gTygzNzcHABQUFKC4uBjNmjVDdHQ0QkJCAAAhISHYvXs3ACAqKgrBwcEwMzODQqGA\ni4sLYmJidBkeERFBx4mgpKQEHh4esLe3h6+vLzp16oT09HTY29sDAOzt7ZGeng4ASElJgVwul9aV\ny+VITk7WZXhERATAVJcbNzExwcWLF/Hw4UP07dsXx48f15gvk8n+LDKXr6J5oaGh0muVSgWVSlUb\n4RIRPTfUajXUanWVltVpIniiadOmGDBgAM6dOwd7e3ukpaXBwcEBqampsLOzAwA4OTkhMTFRWicp\nKQlOTk7lbu/pREBERGU9+yM5LCyswmV1VjV09+5dqUdQXl4eDh8+DE9PTwwcOBCRkZEAgMjISAwa\nNAgAMHDgQGzduhUFBQWIj4/H9evX4ePjo6vwiIjoTzorEaSmpiIkJAQlJSUoKSnBmDFj4OfnB09P\nTwQFBSEiIgIKhQLbt28HACiVSgQFBUGpVMLU1BSrV6+utNqIiIhqh067j+oCu49SfcTuo2RoBus+\nSkREdR8TARGRkWMiICIyckwERERGjomAiMjIMREQERk5JgIiIiPHREBEZOSYCIiIjBwTARGRkdOa\nCG7cuIH8/HwAwPHjx7FixQqNx0sSEVH9pjURDB06FKamprhx4wb+9a9/ITExESNHjtRHbEREpAda\nE4GJiQlMTU3x/fffY9q0aVi0aBFSU1P1ERsREemB1kTQsGFDbN68GRs3bsRrr70GACgsLNR5YERE\npB9aE8H69etx+vRp/Pvf/0br1q1x8+ZNjB49Wh+xERGRHvB5BER6wOcRkKFVdu3U+oSykydPIiws\nDAkJCSgqKpI2ePPmzdqNkoiIDEJriaBDhw74z3/+Ay8vLzRo0ECabmtrq/PgysMSAdVHLBGQodWo\nRGBtbY3+/fvXelBERFQ3aG0s9vX1xcyZM3H69GmcP39e+quKxMRE+Pr6olOnTujcuTNWrFgBAAgN\nDYVcLoenpyc8PT1x4MABaZ2FCxeiXbt2cHV1xaFDh6r5sYiIqKq0Vg2pVKo/i7Wajh8/rnXjaWlp\nSEtLg4eHB3JycuDt7Y3du3dj+/btsLS0xPTp0zWWj42NxciRI3HmzBkkJyejT58+iIuLg4nJX/mK\nVUNUH7FqiAytRlVDarW62jt2cHCAg4MDAMDCwgIdO3ZEcnIyAJQbUFRUFIKDg2FmZgaFQgEXFxfE\nxMSge/fu1Y6BiIgqp7VqKDMzE++99x68vb3h7e2NGTNm4OHDh397RwkJCbhw4YJ0Uf/iiy/g7u6O\nCRMmSGMXpaSkQC6XS+vI5XIpcRARkW5oLRGMHz8ebm5u2LFjB4QQ+Oabb/Dmm2/i+++/r/JOcnJy\nMGzYMCxfvhwWFhaYPHkyPvroIwDAvHnzMGPGDERERJS7bnnVUqGhodJrlUoFlUpV5ViIiIyBWq2u\nco2O1jYCd3d3XLp0Seu0ihQWFuK1115D//798e6775aZn5CQgMDAQFy5cgXh4eEAgNmzZwMA+vXr\nh7CwMHTr1u2vgNlGQPUQ2wjI0Cq7dmqtGnrhhRfw008/Se9PnjwJc3PzKu1YCIEJEyZAqVRqJIGn\nB63btWsX3NzcAAADBw7E1q1bUVBQgPj4eFy/fh0+Pj5V2hcREVWP1qqhNWvWYOzYsVK7QLNmzRAZ\nGVmljZ86dQqbNm1Cly5d4OnpCQBYsGABtmzZgosXL0Imk6F169ZYu3YtAECpVCIoKAhKpRKmpqZY\nvXp1uVVDRERUe6o81lBWVhYAwMrKSqcBacOqIaqPWDVEhlat7qPffPMNxowZgyVLlmj8KhdCQCaT\nlbkHgIiI6qcKE0Fubi4AIDs7m9UzZPSsrGyQnf3A0GEQ6YTWqqGTJ0+iZ8+eWqfpC6uGyBDqQtUO\nq4aoJiq7dmpNBJ6enrhw4YLGNC8vryqPN1TbmAjIEJgIqL6rVhvB6dOn8fPPPyMjIwNLly6VNpCd\nnY3i4mLdREpERHpXYSIoKCiQLvrZ2dnSdCsrK+zcuVMvwRERke5prRpKSEiAQqHQUzjasWqIDIFV\nQ1Tf1Wj0UXNzc7z//vuIjY1FXl6etMFjx47VbpRERGQQWoeYGDVqFFxdXXHz5k2EhoZCoVCga9eu\n+oiNiIj0QGvV0JMeQl26dMHly5cBAF27dsXZs2f1EuCzWDVEhsCqIarvalQ11LBhQwClD5nZu3cv\nWrRogQcPeGMNEdHzQmsi+PDDD5GZmYklS5Zg2rRpyMrKwrJly/QRGxER6YHWRGBtbS39PXnIwcmT\nJ3UdFxER6Um17iwub5q+sI2ADIFtBFTf1fqdxSUlJbqJlIiI9I53FhMZBdMajSJsadkMWVn3azEe\nqkv+1p3F9+/fh7W1NUxMtN5+oDOsGiJDeB6qhli1ZNyq9czisLAwXLt2DQqFAo8fP4avry9cXFzg\n4OCAw4cP6yxYIiLSrwoTwbZt2+Dq6goAiIyMhBACGRkZOHHiBObOnau3AImISLcqTASNGjWS6hQP\nHjyIESNGoEGDBujYsSOKioqqtPHExET4+vqiU6dO6Ny5M1asWAGgtIrJ398f7du3R0BAADIzM6V1\nFi5ciHbt2sHV1RWHDh2qyWcjIqIqqDQRXLlyBRkZGVCr1QgICJDmPXmMpTZmZmZYtmwZfv31V/zy\nyy9YtWoVrl27hvDwcPj7+yMuLg5+fn4IDw8HAMTGxmLbtm2IjY3FwYMHMWXKFPZQIiLSsQoTwX/+\n8x8MGzYMHTp0wHvvvYc2bdoAAPbt2wcvL68qbdzBwQEeHh4AAAsLC3Ts2BHJycmIjo5GSEgIACAk\nJAS7d+8GAERFRSE4OBhmZmZQKBRwcXFBTExMjT4gERFVrsLuo927d8fvv/9eZvqAAQMwYMCAv72j\nhIQEXLhwAd26dUN6ejrs7e0BAPb29khPTwcApKSkoHv37tI6crkcycnJf3tfRERUdVqHmKgNOTk5\nGDp0KJYvXw5LS0uNeTKZrNL+zeXNCw0NlV6rVCqoVKraCpWI6LmgVqulYYG00XkiKCwsxNChQzFm\nzBgMGjQIQGkpIC0tDQ4ODkhNTYWdnR0AwMnJCYmJidK6SUlJcHJyKrPNpxMBERGV9eyP5LCwsAqX\nrbCNYMeOHQCAmzdvVjsQIQQmTJgApVKJd999V5o+cOBAREZGAijtmvokQQwcOBBbt25FQUEB4uPj\ncf36dfj4+FR7/0REpF2FdxY/GViuJgPMnTx5Er169UKXLl2kKp6FCxfCx8cHQUFBuH37NhQKBbZv\n3w5ra2sAwIIFC7B+/XqYmppi+fLl6Nu3r2bAvLOYDIB3FvO8q+8qu3ZWmAj69OkDmUyGM2fO4JVX\nXimzwejo6NqPtAqYCMgQmAh43tV31UoEBQUFOH/+PEaPHo2IiAiNDchkMvTu3Vs30WrBRECGwETA\n866+q1YieCIjIwMvvvgicnJyAJTeD2BITARkCEwEPO/qu2oNOvdEWloaPD09oVQqoVQq4e3tjatX\nr9Z6kET6J2S0AAAU5ElEQVREZBhaE8GkSZOwdOlS3L59G7dv38aSJUswadIkfcRGRER6oDUR5Obm\nwtfXV3qvUqnw6NEjnQZFRET6o/WGstatW+OTTz7BmDFjIITAt99+K407RERE9Z/WEsH69etx584d\nDBkyBEOHDkVGRgbWr1+vj9iIiEgPtPYaqmvYa4gMgb2GeN7VdzXqNURERM83JgIiIiOnNRGcPHmy\nzLRTp07pJBgiItI/rW0E5Q06V5OB6GqKbQRkCGwj4HlX31V27ayw++jp06fx888/IyMjA0uXLpU2\nkJ2dzecIExE9RypMBAUFBcjOzkZxcTGys7Ol6VZWVti5c6degiMiIt3TWjWUkJAAhUKhp3C0Y9UQ\nGQKrhnje1XfVqhp64vHjx3jrrbeQkJCAoqIiaYPHjh2r3SiJiMggtJYIunTpgsmTJ8PLywsNGjQo\nXUkmg7e3t14CfBZLBGQILBHwvKvvalQiMDMzw+TJk2s9KCIiqhu03kcQGBiIVatWITU1Fffv35f+\niIjo+aC1akihUEgPnn9afHy81o2PHz8e+/btg52dHa5cuQIACA0NxVdffYUXX3wRQOnD6vv37w+g\n9MH269evR4MGDbBixQoEBASUDZhVQ2QArBrieVff1ehRlTXx008/wcLCAmPHjpUSQVhYGCwtLTF9\n+nSNZWNjYzFy5EicOXMGycnJ6NOnD+Li4mBiolloYSIgQ2Ai4HlX39WojSAyMrLcEsHYsWO17viV\nV15BQkJCmenlBRMVFYXg4GCYmZlBoVDAxcUFMTEx6N69u9b9EBFR9WlNBGfOnJESQV5eHo4dOwYv\nL68qJYKKfPHFF9i4cSO6du2KJUuWwNraGikpKRoXfblcjuTk5Grvg4iIqkZrIli5cqXG+8zMTAwf\nPrzaO5w8eTI++ugjAMC8efMwY8YMRERElLtseSURoLSd4QmVSgWVSlXteIiInkdqtRpqtbpKy2pN\nBM8yNzevUkNxRezs7KTXEydORGBgIADAyckJiYmJ0rykpCQ4OTmVu42nEwEREZX17I/ksLCwCpfV\nmgieXKgBoKSkBLGxsQgKCqp2cKmpqXB0dAQA7Nq1C25ubgCAgQMHYuTIkZg+fTqSk5Nx/fp1+Pj4\nVHs/RERUNVoTwYwZMwCUVtOYmpqiZcuWcHZ2rtLGg4ODceLECdy9exfOzs4ICwuDWq3GxYsXIZPJ\n0Lp1a6xduxYAoFQqERQUBKVSCVNTU6xevbrCqiEiIqo9Veo+mpaWJjUa+/j4aFTv6Bu7j5IhsPso\nz7v6rkbPLN6+fTu6deuGHTt2YPv27fDx8cGOHTtqPUgiIjKMKg06d+TIEakUkJGRAT8/P1y+fFkv\nAT6LJQIyBJYIeN7VdzUqEQghpOEgAKB58+b8QhARPUe0Nhb369cPffv2xciRIyGEwLZt26SxgYiI\nqP6rUmPxd999h1OnTgEoHTZi8ODBOg+sIqwaIkNg1RDPu/quWoPOXb9+Henp6ejZs6fG9JMnT8LR\n0RFt27at/UirgImADIGJgOddfVetNoJ3330XVlZWZaZbWVnh3Xffrb3oiIjIoCpMBOnp6ejSpUuZ\n6V26dKnREBNERFS3VJgIMjMzK1wpPz9fJ8EQEZH+VZgIunbtinXr1pWZ/uWXXxrswfVERFT7Kmws\nTktLw+DBg9GwYUPpwn/u3Dk8fvwYu3btkgaO0zc2FpMhsLGY5119V+1HVQohcPz4cVy9ehUymQyd\nOnXCP//5T50FWhVMBGQITAQ87+o7gz2zWBeYCMgQmAh43tV3NRpigoiInm9MBERERo6JgIjIyDER\nEBEZOSYCIiIjx0RARGTkdJoIxo8fD3t7e7i5uUnT7t+/D39/f7Rv3x4BAQEaQ1ksXLgQ7dq1g6ur\nKw4dOqTL0IiI6E86TQRvvvkmDh48qDEtPDwc/v7+iIuLg5+fH8LDwwEAsbGx2LZtG2JjY3Hw4EFM\nmTIFJSUlugyPiIig40TwyiuvoFmzZhrToqOjERISAgAICQnB7t27AQBRUVEIDg6GmZkZFAoFXFxc\nEBMTo8vwiIgIBmgjSE9Ph729PQDA3t4e6enpAICUlBTI5XJpOblcjuTkZH2HR0TlMoVMJqvRn5WV\njaE/BFVA6zOLdenJF6Sy+eUJDQ2VXqtUKqhUqlqOjIg0FaFmQ1QA2dkVn+tU+9RqNdRqdZWW1Xsi\nsLe3R1paGhwcHJCamgo7OzsAgJOTExITE6XlkpKS4OTkVO42nk4ERERU1rM/ksPCwipcVu9VQwMH\nDkRkZCQAIDIyEoMGDZKmb926FQUFBYiPj8f169fh4+Oj7/CIiIyOTksEwcHBOHHiBO7evQtnZ2d8\n/PHHmD17NoKCghAREQGFQoHt27cDAJRKJYKCgqBUKmFqaorVq1dXWm1ERES1g8NQE1UBh6Gu6fql\n2+C5azgchpqIiCrEREBEZOSYCIiIjJxB7yMg0hcrKxtkZz8wdBhEdRIbi8kosLHX0OuXboPnruGw\nsZiIiCrEREBEZOSYCIiIjBwTARGRkWMiICIyckwERERGjomAiMjIMREQERk5JgIiIiPHREBEZOSY\nCIiIjBwTARGRkWMiICIycgYbhlqhUMDKygoNGjSAmZkZYmJicP/+fQwfPhy3bt2SnmdsbW1tqBCJ\niIyCwUoEMpkMarUaFy5cQExMDAAgPDwc/v7+iIuLg5+fH8LDww0VHhGR0TBo1dCzY2NHR0cjJCQE\nABASEoLdu3cbIiyqg6ysbCCTyar9R0QVM2iJoE+fPujatSu+/PJLAEB6ejrs7e0BAPb29khPTzdU\neFTHlD5dTNTgjwzPtEbJ3MrKxtAf4LllsDaCU6dOwdHRERkZGfD394erq6vGfP6SI3reFKEmSTk7\nm9cDXTFYInB0dAQAvPjiixg8eDBiYmJgb2+PtLQ0ODg4IDU1FXZ2duWuGxoaKr1WqVRQqVR6iJiI\nqP5Qq9VQq9VVWtYgzyzOzc1FcXExLC0t8ejRIwQEBGD+/Pk4cuQImjdvjlmzZiE8PByZmZllGoz5\nzGLjxGcO1/f1aycGnvvVV9m10yCJID4+HoMHDwYAFBUVYdSoUZgzZw7u37+PoKAg3L59u8Luo0wE\nxomJoL6vXzsx8NyvvjqXCGqCicA4MRHU9/VrJwae+9VX2bWTdxYTERk5JgIiIiPHREBEZOSYCIiI\njBwTARGRkTPYDWVkXKysbP4cJoKI6hp2HyW9YPdPY1+/dmLguV997D5KRM8BDlqnK6waIqJ6goPW\n6QoTAWnF+n2i5xvbCEirmtfvA4avo+b6z0MbQc3WN0NpqaJ6LC2bISvrfg32b1iVXTtZIiAiI8Gq\npYqwsZiIyMgxERARGTkmgnqgpg9uZ7c5IqoMG4vrgdq4Gasmx4yNxVz/+WgsNu4b2nhDGRERVYiJ\ngIioSp7fO5vZfVQPDH9Dlumf1TtEVH3Pb/fTOlciOHjwIFxdXdGuXTt8/vnnhg6nVpQmAVGDv5p6\n8gU21P6JqKZq2mmkMnUqERQXF2Pq1Kk4ePAgYmNjsWXLFly7ds3QYen0P+D5oDZ0AHWI2tAB1BFq\nQwdQh6j//LdmVUu6/EFZpxJBTEwMXFxcoFAoYGZmhhEjRiAqKqrMcvqumyv/P2B+OdOM9Re12tAB\n1CFqQwdQR6gNHUAdov7z37pbMq9TiSA5ORnOzs7Se7lcjuTk5HKW/HsHMDs7m7/oiYgqUKcai3V3\n0a1ZI09p/2MioudTnUoETk5OSExMlN4nJiZCLpdrLOPu7o5Ll6pzYa7pxby89cMMvP+6tH5VjkVd\n/wy1tX5Fx6K+xF9b65d3HOrbZ6it9Z8cC8PF7+7uXvFW69KdxUVFRejQoQOOHj2KFi1awMfHB1u2\nbEHHjh0NHRoR0XOrTpUITE1NsXLlSvTt2xfFxcWYMGECkwARkY7VqRIBERHpX53qNVQfLVmyBCYm\nJrh/v/4+uaimZs6ciY4dO8Ld3R1DhgzBw4cPDR2S3j2PN0JWR2JiInx9fdGpUyd07twZK1asMHRI\nBldcXAxPT08EBgYaOpQKMRHUQGJiIg4fPoxWrVoZOhSDCggIwK+//opLly6hffv2WLhwoaFD0qu6\neiOkIZiZmWHZsmX49ddf8csvv2DVqlVGeyyeWL58OZRKZZ3uis5EUAPTp0/H//7v/xo6DIPz9/eH\niUnpV6lbt25ISkoycET6VdUbIY2Bg4MDPDw8AAAWFhbo2LEjUlJSDByV4SQlJWH//v2YOHFinR7C\nmomgmqKioiCXy9GlSxdDh1KnrF+/Hq+++qqhw9Crqt8IaVwSEhJw4cIFdOvWzdChGMx7772HRYsW\nST+U6qo61WuorvH390daWlqZ6Z999hkWLlyIQ4cOSdPqcravDRUdiwULFkh1n5999hkaNmyIkSNH\n6js8g6rLRX5DycnJwbBhw7B8+XJYWFgYOhyD2Lt3L+zs7ODp6Qm1Wm3ocCrFRFCJw4cPlzv96tWr\niI+Pl27QSEpKgre3N2JiYmBnZ6fPEPWmomPxxNdff439+/fj6NGjeoqo7qjKjZDGpLCwEEOHDsXo\n0aMxaNAgQ4djMD///DOio6Oxf/9+5OfnIysrC2PHjsXGjRsNHVoZ7D5aC1q3bo1z587BxqbuPnhC\nlw4ePIgZM2bgxIkTsLW1NXQ4escbIf8ihEBISAiaN2+OZcuWGTqcOuPEiRNYvHgx9uzZY+hQylW3\nK67qCWOvGpg2bRpycnLg7+8PT09PTJkyxdAh6dXTN0IqlUoMHz7cKJMAAJw6dQqbNm3C8ePH4enp\nCU9PTxw8eNDQYdUJdfk6wRIBEZGRY4mAiMjIMREQERk5JgIiIiPHREBEZOSYCIiIjBwTARGRkWMi\nIINo0KABPD090blzZ3h4eGDp0qVah+lISEiAm5sbAODSpUs4cOBAlfd35coVqV978+bN0aZNG3h6\neiIgIKBGn0NXoqKidDpq5507dzBgwADpfUxMDHr16gVXV1d4eXnhrbfeQl5eHqKjo/HJJ5/oLA6q\nG5gIyCDMzc1x4cIFXL16FYcPH8aBAwcQFlb1Z0BfuHAB+/fvr/Lybm5uuHDhAi5cuICBAwdi8eLF\nuHDhgsZ4UfpWUlJS4bxdu3YhNjb2b22vqKioysuuXLkS48aNAwCkp6cjKCgIixYtwm+//Ybz58+j\nX79+yM7ORmBgIL777jsUFhb+rViofmEiIIN78cUXsW7dOqxcuRJA6fj+M2fOhI+PD9zd3bFu3TqN\n5QsLC/HRRx9h27Zt8PT0xPbt23HmzBn06NEDXl5eePnllxEXF1fpPp+UPg4dOoQePXrA29sbQUFB\nePToEQBAoVBg7ty58PT0RNeuXXH+/HkEBATAxcUFa9euBQCo1Wr06tULr732GlxdXTF58uQqbXf2\n7Nnw9vbGjh078NVXX8HHxwceHh4YNmwY8vLy8PPPP2PPnj2YOXMmvLy8cPPmTahUKpw7dw4AcPfu\nXbRu3RpA6RhPAwcOhJ+fH/z9/ZGbm4vx48ejW7du8PLyQnR0dLmff+fOnVKJYNWqVRg3bpzGKKFD\nhw6FnZ0dZDIZXnrpJYMmTNIDQWQAFhYWZaZZW1uL9PR0sXbtWvHpp58KIYTIz88XXbt2FfHx8SI+\nPl507txZCCHE119/LaZNmyatm5WVJYqKioQQQhw+fFgMHTq0wn2PGzdOfPfddyIjI0P06tVL5Obm\nCiGECA8PFx9//LEQQgiFQiHWrFkjhBDivffeE25ubiInJ0dkZGQIe3t7IYQQx48fF40bNxbx8fGi\nuLhY+Pv7i507d2rd7qJFi6RY7t27J73+8MMPxRdffKER4xMqlUqcO3dOCCFERkaGUCgUQgghNmzY\nIORyuXjw4IEQQog5c+aITZs2CSGEePDggWjfvr149OiRxudPTU2VjqMQQgwZMkRER0dXeLzWr18v\nPvjggwrnU/3H0Uepzjl06BCuXLmCnTt3AgCysrJw48YNuLi4SMsIITTaFDIzMzF27FjcuHEDMplM\na1WGEAK//PILYmNj0aNHDwBAQUGB9BoABg4cCKC0WunRo0do0qQJmjRpgkaNGiErKwsA4OPjA4VC\nAQAIDg7GyZMn0bhx40q3O3z4cOn1lStX8OGHH+Lhw4fIyclBv379NGKsCn9/f1hbW0vHbs+ePVi8\neDEA4PHjx0hMTESHDh2k5W/dugVHR8cyx6MiLVq04HhBzzkmAqoTbt68iQYNGkjDeK9cuRL+/v4a\nyyQkJFS4/rx58+Dn54ddu3bh1q1bUKlUVdqvv78/Nm/eXO68Ro0aAQBMTEzQsGFDabqJiYlUH//0\nQGJCCMhkMgghKt1ukyZNpNfjxo1DdHQ03NzcEBkZqTFu/dPbNjU1ldoU8vPzK9weAHz//fdo165d\nhZ/5SaxPdOrUCefOnZMS37NKSkrq9IBpVHNsIyCDy8jIwNtvv41p06YBAPr27YvVq1dLF9u4uDjk\n5uZqrGNlZYXs7GzpfVZWFlq0aAEA2LBhg9Z9ymQydO/eHadOncIff/wBAHj06BGuX79eZtnKfi3H\nxMQgISEBJSUl2L59O1555ZUqbxcofYCLg4MDCgsLsWnTJumCa2lpKZU6gNK2hbNnzwKAVFIqT9++\nfTUeGH/hwoUyy7Rq1UrjIUNTp05FZGQkYmJipGnff/897ty5AwBITU01+udyP++YCMgg8vLypO6j\n/v7+6NevHz766CMAwMSJE6FUKuHl5QU3NzdMnjwZxcXFAP76lezr64vY2FipsfiDDz7AnDlz4OXl\nheLi4ir9grW1tcXXX3+N4OBguLu7o0ePHvj999/LLCeTyTS29/Trf/zjH5g6dSqUSiXatGmDwYMH\nV3m7APDJJ5+gW7du6Nmzp8bQ1SNGjMCiRYvg7e2N+Ph4vP/++/jvf/8LLy8v3Lt3T4rh2djmzZuH\nwsJCdOnSBZ07d8b8+fPL7NPBwQFFRUVSA7adnR22bt2K999/H66urlAqlTh8+DCsrKwA/NW1lJ5f\nHIaaqJrUajWWLFlSZx82UpnQ0FB07NhRo72iPCUlJfDy8sLZs2dhasqa5OcVSwRE1fTsr/H65J13\n3kFkZKTW5fbu3Ythw4YxCTznWCIgIjJyLBEQERk5JgIiIiPHREBEZOSYCIiIjBwTARGRkWMiICIy\ncv8fODqYcsc6ooIAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mapping the Data\n", "\n", "The final step is mapping the data. There are three major steps to this process: (1) Setting the color of each marker, (2) creating the temperature over time graph for each marker, and (3) plotting the markers on the map.\n", "\n", "#### Coloring the Markers\n", "\n", "We want to assign a color, blue or red, to each marker depending on its change over the period. The first step is to discretize the deltas into 7 equally-sized categories. Because most stations had a change of less than +/-5° C, we will use that as the limits of our scale. The first two lines find any changes of greater than 5° and cap them. After that, the [`cut`](http://pandas.pydata.org/pandas-docs/version/0.13.1/generated/pandas.cut.html) method is used to discretize the `delta` column into 7 categories. I chose a 7-tone blue-to-red scale from [ColorBrewer](http://colorbrewer2.org/) to represent the colors of the markers, which will be assigned based on the discretized delta." ] }, { "cell_type": "code", "collapsed": false, "input": [ "change.delta[change.delta > 5] = 5\n", "change.delta[change.delta < -5] = -5\n", "\n", "change['cat'] = pd.cut(change['delta'], 7)\n", "cats = change.sort('delta')['cat'].unique()\n", "colors = [\"#2166ac\",\"#67a9cf\",\"#d1e5f0\",\"#f7f7f7\",\"#fddbc7\",\"#ef8a62\",\"#b2182b\"]\n", "cat_colors = {cat:color for cat, color in zip(cats,colors)}\n", "list(cats)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 11, "text": [ "['(-5.01, -3.571]',\n", " '(-3.571, -2.143]',\n", " '(-2.143, -0.714]',\n", " '(-0.714, 0.714]',\n", " '(0.714, 2.143]',\n", " '(2.143, 3.571]',\n", " '(3.571, 5]']" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Creating the Charts\n", "Visualization is done with [folium](https://folium.readthedocs.org/en/latest/) for the map and [vincent](https://vincent.readthedocs.org/en/latest/) for charts. These both have simple and well-documented APIs. \n", "\n", "To create the graph for each station, we first select the readings that pertain to the station of interest, using the familiar bracket syntax for selection. The next line does several things.\n", "```python\n", "weather_frame = weather_frame.rename(columns={'tmax':'Max','tmin':'Min','tavg':'Average'}).set_index('date')\n", "```\n", "\n", "We first have to rename the columns to human-friendly names like *Max*, *Min* and *Average*. Vincent will use these column names in its legend. Next, we set the dataframe index to `date`. This lets vincent know that we want a time-series graph with the date on the X axis.\n", "\n", "Vincent doesn't need any more configuration than that, so we can simply create a line graph with `chart = vincent.Line(weather_frame)`. The rest of the code performs fairly standard formatting of the chart, and is all found in the vincent API linked above." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import folium, vincent\n", "vincent.core.initialize_notebook()\n", "\n", "def popup_chart(station_id):\n", " weather_frame = data[data.id == station_id][['date','tmax','tavg','tmin']]\n", " weather_frame = weather_frame.rename(columns={'tmax':'Max','tmin':'Min','tavg':'Average'}).set_index('date')\n", " station = index[index.id == station_id].iloc[0] # Find the station in the index\n", "\n", " chart = vincent.Line(weather_frame)\n", " chart.axis_titles(x=station['name']+\", \"+station[\"country\"],y=\"Temp (C)\")\n", " chart.legend(title=\"Key\")\n", " chart.width=400\n", " chart.height=200\n", " chart.title=station['name']+\", \"+station[\"country\"]\n", " chart.padding={'top':20,'right':80,'bottom':40,'left':40}\n", "\n", " chart.to_json(\"viz/historical_map/data_{0}.json\".format(station['id']))\n", " return chart \n", "\n", "chart = popup_chart(42572734000)\n", "chart.display()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", " " ], "metadata": {}, "output_type": "display_data", "text": [ "" ] }, { "html": [ "
\n", "\n", "\n", " " ], "metadata": {}, "output_type": "display_data", "text": [ "" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Putting it All Together: the Map\n", "\n", "With folium, creating the map itself is fairly simple! The configuration options, such as the default zoom level and tile set to use, are found in the [folium docs](http://folium.readthedocs.org/en/latest/).\n", "\n", "We then iterate through all our stations, and add a marker for each one. The popup is set by generating a vincent chart in the above manner for each data point. The fill color is set with the previously generated `cat_colors` dictionary. We save the final map to disk, and we are good to go!\n", "\n", "(Note: the final map won't show up in the online ipython notebook viewer because it relies on Javascript. However, you can see it [here](http://corbt.s3-website-us-east-1.amazonaws.com/weather/historical_map/map.html).)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "folium.initialize_notebook()\n", "f_map = folium.Map(\n", " location=[0,0],\n", " zoom_start=2,\n", " tiles=\"Mapbox Bright\")\n", "\n", "for station_id,station in change.iterrows():\n", " f_map.polygon_marker([station['lat'],station['long']],\n", " radius=5,\n", " popup=(popup_chart(station_id),\"data_{0}.json\".format(station_id)),\n", " line_color='grey',\n", " line_weight=1,\n", " fill_color=cat_colors[station['cat']],\n", " fill_opacity=1\n", " )\n", "\n", "f_map.create_map(path='viz/historical_map/map.html')\n", "f_map" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] }, { "html": [ "" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] }, { "html": [ "\n", "
\n", "\n", "" ], "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "" ] } ], "prompt_number": 13 } ], "metadata": {} } ] }