api-examples/ERA5_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ERA5 tutorial\n",
    "\n",
    "ERA5 contains historical weather data, which can be used to analyse very wide range of problems. In this tutorial we briefly demonstrate how to:\n",
    "1. Get timerange data for single location and single variable\n",
    "2. Analyse possible wind farm production potential for three different locations\n",
    "\n",
    "To read more about the ERA5 dataset, please follow these links:\n",
    "1. https://data.planetos.com/datasets/ecmwf_era5\n",
    "2. https://software.ecmwf.int/wiki/display/CKB\n",
    "\n",
    "To learn more about Planet OS/Intertrust datahub, please refer to documentation in https://data.planetos.com/datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialize notebook environment.\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Import Planet OS API\n",
    "from API_client.python.datahub import datahub_main\n",
    "from API_client.python.lib.dataset import dataset\n",
    "from API_client.python import package_api"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize dataset, print available variable names. In order to be able to access the data, save your API key in the file `APIKEY` and put it into the same folder where you run this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "apikey = open('APIKEY').readlines()[0].strip()\n",
    "dh = datahub_main(apikey)\n",
    "ds = dataset('ecmwf_era5',dh,debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['100_metre_U_wind_component_surface',\n",
       " '100_metre_V_wind_component_surface',\n",
       " '10_metre_U_wind_component_surface',\n",
       " '10_metre_V_wind_component_surface',\n",
       " '2_metre_dewpoint_temperature_surface',\n",
       " '2_metre_temperature_surface',\n",
       " 'Maximum_temperature_at_2_metres_surface_1hour_Maximum',\n",
       " 'Mean_sea_level_pressure_surface',\n",
       " 'Minimum_temperature_at_2_metres_surface_1hour_Minimum',\n",
       " 'Snow_density_surface',\n",
       " 'Snow_depth_surface',\n",
       " 'Surface_pressure_surface',\n",
       " 'Surface_solar_radiation_downwards_surface_1hour_Accumulation',\n",
       " 'Total_precipitation_surface_1hour_Accumulation',\n",
       " 'lat',\n",
       " 'lon',\n",
       " 'time0',\n",
       " 'time1']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.variable_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Choose location, variable\n",
    "lon = 26\n",
    "lat = 58\n",
    "variable = '2_metre_temperature_surface'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "d1 = ds.get_json_data_in_pandas(count=100000,**{'vars':'2_metre_temperature_surface','lon':lon,'lat':lat})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAADyCAYAAABakyJyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXeYVNX5xz/v7tKWXpeqSJdmYQUsyKKiIhpb9KcplsQQjTUxsUZj7NHEGFsSYzcqGo1iRLEgK70KgggoSpXeWcrC7p7fH3dmd3b2zsy9M7ec2T2f59mHmTt35n655bznvOc97ytKKQwGg8FgcEpO2AIMBoPBkF0Yw2EwGAwGVxjDYTAYDAZXGMNhMBgMBlcYw2EwGAwGVxjDYTAYDAZXGMNhMBgMBlcYw2EwGAwGVxjDYTAYDAZXGMNhMBgMBlfkhS3AD9q0aaO6du2a1nf37NlD48aNvRXkAUaXe3TVZnS5w+hyRya65s2bt0Up1TbljkqpWvc3aNAglS6TJk1K+7t+YnS5R1dtRpc7jC53ZKILmKsctLHGVWUwGAwGVxjDYTAYDAZXGMNhMBgMBlcYw2EwGAwGVxjDYTAYDAZXGMMREPsOlLNtz4GwZRgMBkPGGMMREOc+NY2j7/k4bBkGg8GQMcZwBMTSDbvDlmAwGAyeoI3hEJGGIjJbRL4QkcUi8sfI9sNEZJaIfCMir4tI/bC1ZkJZeUXYEgwGgyEjtDEcQClwklLqCOBI4HQRGQr8CfirUqonsB34eYgaM+aY+z4JW4LBYDBkhDaGI7LivSTytl7kTwEnAW9Gtr8InBOCPM/Yvvdg2BIMBoMhI7RKcigiucA8oAfwJPAtsEMpVRbZZS3QKSR5rlFK0e8PHzLo0JZhSzEYDAbPECuvlV6ISAvgbeBO4HmlVI/I9i7A+0qpATbfGQOMASgoKBg0duzYtI5dUlJCkyZN0pVejfIKxc8/2ltj+wunu89c6aUuL/FD15wNZTy5oJSnTs4nv56k/Tt+aFu2rZwHZu/nr0WNaNkwvQF7XbqWXmB0uSMTXSNGjJinlCpMtZ9WI44oSqkdIlIMDAVaiEheZNTRGViX4DtPA08DFBYWqqKiorSOXVxcTLrfjeehCUuxBk3VSef3vdQFsHrrXu57/yv+dtFRNKyXm/bveK0L4LJbxgNQr1NfivoWpP07fmh789XPgfXktO9N0ZHpDX790PWPz75l34Fyfj2yV9q/4YcuLzC63BGELm3mOESkbWSkgYg0Ak4BlgCTgB9GdrsUGBeOQvc8VVzTaOjC7e8s4sPFG5n89eawpSRkw859YUtIiEj6IyGvmbR0Ew9+sJS/TfwmbCnacqCsgv0Hy8OWUWvQxnAAHYBJIrIQmAN8rJR6D7gZ+I2ILAdaA8+GqLHWMOWbLQC8NGNVyEoSk5OjT+McJSzH7pmPT+Gsx6fafnb92PkBq8k+ev3+A/rcMSFsGbUGbVxVSqmFwFE2278DBgevqHby4eIN/PLleZXvd+7TN8orV6NefTxBK/vy+10JP9u1vyzhZ2Gz90AZk5ZuZvTADmFLyRo++WojB8orOGOAvudMG8Nh8Jdzn5rGOUd2YvaKbdW2L/p+Z0iKUpOjgeGYvnwLzRrVo3+n5tW26xdSoid97/wQgLXb+/DL4d1DVpMdXPHSXABWPjg6ZCWJ0clVZfCR+at38Id3F4ctwxUa2A1+9Mwsznx8Kht27q+2PSxpqTIPrNlWM4pPBx74YGnYEgweYgxHHUNlUV85V6M5jqEPTATg++3WhH1YZ/HhD5cl/fzKf89L+rnBWl+1eJ2+I+0ob81bG7aEhBjDUcd4f9GGsCUk5M5xXzIoJoOwDq6qeBas2QGE17NfvC7xXIeTzw3w7NQVjH5sag23rW7c/NbCsCUkxBgOgza8NGMVW2NqlmhoNyoJa+HstG+3VL7+y0fLuPS52aHoSIePFuvRaYkaV13delF0vv+N4agD/GfumrAlpEUYayUOlldQXqG45LnZfPLVxmqfxc4vhJVwQSmYu9LqKT/+6XI+03gdTjxjXg7fjbbvQPas5ZDQZtJSYwxHHeB3byYf8s5bpeeQPYzHpuftHzDo3o+Z/PXmyuiWKH+aoMcE763/XRS2hKRc9e95jFvwfdgy+O/na+l6y3jWxywkvfrVz0NU5BJ97YYxHAY4/+8zwpZgS1hzHDsSZDBOtpYiSHQOb6ioUHzw5QauH7sgbCm8GZlcHju7asQ9bfkW3p4fvlFzgsZ2wxgOv9hu6otnzLKNwVZNLC1L7saItWM6N95hslWj+376t1sBTCoWHzCGwyfOTJAewuCcxwJ+4Hv/PnlKimhDBLDnQHirtXXMaJ0tZNOZ0zGqMIoxHB5wsLyCjbuqLxD7foe+Cfrs6BrJSGtwxr8mfxfasY3d8IaNu/czf/V2vgl4ZOsUje2GMRxe8Pu3v2TI/RPZG2IvNNtYs20v05ZXhZZWVGRXa5hlckNB91DhXfvKOPep6Yz86+SwpVTyzJSqDonGdsMYDi94f9F6AP6hcRp1J1zwj+mBuUGGPTSJHz8zq/L9c9NWBHJcLwnL2H23ZU8ox3WL7qHCOvbo7x2/JGwJjjCGwwN2l1ojjcc+XR6yksyYs3J7aMdetkFPd0Ey3ot0GAzWvMuL01c6yrZcVl4R2jzNgbJw1+Js3LW/MvtAKvZovObEGA5DNcKKislGz0+JxunMg+bz1dv5w7uLOeWRz5Lut37nPnrc/gGvzdZrUeqXAWWJLnq4mHOenBbIsfzEGA5DNdbv2J96JwMAq7Zmh8soCLaUOOtwRCtOvvtF+Gsp/vFZlWv589XBjLb3RaoQvjJL3wJqTjCGwweyuUGpMCE7jpn8zZbUO9URnIaOvhoZaeiWTiNoNbe//SUHyyvYU5qdo1ZTyMlj/vnZt4zo0y5sGQB8/NVGVrtM5FZuDIdj9Gr6wsXpuYg2lH5PTK/e6u6+X+lyfy/oefsH1rE1LtiUCGM4POaBD5ZqU7TmF3G5lpwQVqRQNtorHaNydKdP+6Ys31Ti67nbd6CcbzeXuPrOs1NXcMeZfX1SlJpXZ62ulvlYd7RxVYlIFxGZJCJLRGSxiFwf2d5KRD4WkW8i/7YMW2umbC0pDVtCQoK2G9HommxcDa1z7YvyChXoOXVqCIKQdMlzs7j8hTmuv3fZ8+GsO5m7chu3vb2I8QtrRund+l89a3JoYziAMuBGpdThwFDgahHpC9wCTFRK9QQmRt5nNYPvnxi2hIQEMccRGxJZUlqGUoqP41KYR+l6y3jWbNvLii17sjJk1y/Kk1j4HXsP0P2293lmSnBrY5LpiSV6f01bvtW3fG7phpUXLwtn3UmyLBOvzV7Drv2pQ5yDRhvDoZRar5T6PPJ6N7AE6AScDbwY2e1F4JxwFNqTTmoRpw9ZGGzctd9RLH4mHPtAleFUwGG3vl+5FsaOZ6euYMSfizntUX1W+EYJy5jt2Ju40d0QSX/z4oyVwYgBxs5xFl4b2zFxO/9WW5m+fGvSzwfe9VFASpyjjeGIRUS6AkcBs4ACpdR6sIwLoMfMc4RsjYpIxPVjFzDU5xFR7FqR7zanjkDTeS7htEcna1F7wo6124PJl6aU4tOlmxztu+9g1WhT5+saJK9nYaE17SbHRaQJ8BZwg1Jql9MqcCIyBhgDUFBQQHFxcVrHLykpcfXdNbsrUu9kg1t9bnVlwr6D5Y6Plamu1yem9kWvXbu28rWbYwV1zq4fu4DmO5xn8vVC17Tp0xN+9uKEqlQuQZyvT1c7H6F+vbZqAnjevHlsW57rm6500PH+Av10aWU4RKQeltF4RSn138jmjSLSQSm1XkQ6ALZdG6XU08DTAIWFhaqoqCgtDcXFxbj57lfrdsG0Ka6P41afW10ATEg/4+1lE/bwy+HduHXU4b7qem1paj9363YdYJXVK3NzLNfaMjhffQcNpV3Tho729eJaHn/ccfDpJ7a7xp5TX89XhGefnQU4iwjaXZYDWIvgBg0axMDOLbzXlcF19P18paktiOvoBm1cVWINLZ4FliilHon56F3g0sjrS4FxQWtLhsrKZBnO+Odn4aUOjyXWf971lvG8PX9tkr3D4ewnvE8j0fWW8VmR7j43x7nPye/8SyZDdTBoYziA44GfAieJyILI3xnAg8BIEfkGGBl5rw2ZBCGt37mPD3xIlFdaVs4f/7fY89/VhXcXrAtbQg3W76y7qVo+XxVecsx4Nu/WN9S9NqGNq0opNZXEC1BPDlKLG9ZkEBnyw7/P4Psd+1jxwBk4nctxwhtz1/L8tJWe/Z5u7DbJBbVKQ78rzeuRhUt3QuP5aSu4/PjDwpZRiU4jjqykLIPQ2mgor9fRudleFyQVc33o4eocIv0zm8VsT05ydo1fnrHSWzGas6GWjvz++L+v2KZRPXdjODIk3ebmPz6G4GVb2VodiGYt1RGnoa523DGu9ros7ZiXYZZbHRfbRdEpu4IxHBmS7sX83ZtVqQR0uiHqKiYrcHos36TXav5Ms+7e+t9FHinxnky8G15jDEeG+L3K2hAM2VbzXAfuencxpzyi12r+jbsyc1XpPLn+0oyVYUuoxBiODLnTA1dAXWmynpu6QtvwUmM33PPC9JVhS6iBjpqivJGhe/qVWas9UpI5xnBkgOmluuPu974KW0JCdJp4NLhj/8FyBt/3CcXL0p8LCoKb3sws0+2Ovfp4N7QJx81G3p7vTY4i414Pl8J7P2GLxqnudaO8QtH9tvcz/h2vbvvvNu9h0+5SHtSkDk5dwIw4MsCrXmptXn2eDRij4Y6D5enlZ4vHq4CE6PPjxVqo2Su2ccWL7mt51DWM4TAYDK5wWl88FTO+TZ5O3ClR++PVEtpPlnjn8jpQVlEr06AYw6EBOruqZn7nzcPtNf/7Yp0noaB1IRTabRnVVGzd480I7eEPl3nyO9GRS46GrdlZj0+l750f1rr7TMNTrT8/eWYW05dnT33gTNC12M61r83XLhRUV07+y2ee/dbyTbt5bbZ3i1e9cBNGY1S8Ggl5ybKNVuemtgVfODIcIpIvIo39FpMNbN9zgKnLt/CjZ2al3jnCzzTKMeMWrx7GVVtTF2wypE+3NskfT6/CoE95ZDKPTXReeyQVC9fuyPg3or35hWt3ZvxbflFe10YcInIdVrrzf4nIr/2XpDfpTOiN6NPWByXB4CJjdlJG1oHRQZhp0CfccGIox9UBXaPi19Xi1D9ORhyXAj8CfoKV9rxOc7Dc/V06rGdb3rv2hISfl9Sy8rN2HPAoEsdr/NA1P8N8SelQP6/uep11nT+48J8zPP/NX7++gDcc1nf3Eyd328PAf7Eq8z3qrxz9iQ1FfN5Fauv+nZpzwaDOtp/993PvChN5FakSpbbVVI9n/urMXSXxnPvUdG2DCmojuo44YtOf7C31Jonm2/O/56a3MltI6AUpDYdSaqxS6tzI30tBiNKZ2M7NOgcpnOvnVp3iRNMFXubHmeGiwerUolHKfdKttVDXWbDGe4OUiC6trOv4r0sKAzumTvgx4vC6A6Zz1oR0cDLHcYeI3BiEmGzA7RzHqAHtK18nmmieF0IFtZF9C7jsuK6BH1c33E7hnDmwg6P9glzFPOWmkwBo3qheYMfUCT9GHDv2Zh4FFdtUbK2DUVU/Bf4ev1FErhCRW72XpDduDcfdP+hf+TrRytZ05k0S4bQhHNW/PfVyU+89do4+idX84OnJ7uqqP3LhkT4pseh+2/vcnGZOowYp5jk2ZZg5Vlc2Z8PKf03nYdLFieHYp5SyC+Z/GWvCvE7htnfTPL+qF5jIVXWgzLsJ2uWbnC32Ou/ozuTlpr78a7bV3sgQcF+Gtn5eDkcd0sInNVYeqNfnrmHZBveLGwd2bp7088H3T+Sud9PL5rx8U4kvLqFvNpawp7SM3/7nC3ammcRvu8ve/N1n90u5jxdR6LXLVFTHkeEQkRrjc6VUKVDnHOCZPDyJQluXe7iy182CvVH926feqZZT0Lyho/2OOqQFFw/uAsDbvzreT0kAaWV6dZKrKZ204y/PXMUpj3zmi/vtgy838PLMVbw5by1PFS93/f3SsnL+4NIYXnJsV6bfcpLrY2WEhosTM8GJ4fgLME5EDo3dKCLtAE9jGUXkORHZJCJfxmxrJSIfi8g3kX9benlMt2TiTz36EHvp8fWug8pt07pJA0f7jfrbFMYt8CYTsG446QjMvPVk3v7V8Txw3sAAFFno1Ft94P0lAPzTpVsvCNykGr/suK4suutUADo6CAzJFF3DhL3ASVTVf4AngXki8p6I3Csi9wPTgD97rOcF4PS4bbcAE5VSPYGJkfehkUlGz/OOtg/HjeWrdbvoe+eHvLdwXVrHiDdCqXjrquNS7rNk/S6uH7sgLT0VFSqtB+i2M/qkdTy3OFHW3uGoxEsOeui+zBSd67G7ubXu+kE/mjYMLoAg3ZZi7JihnurwA0erhpRSLwKHAW8A9YD9wMVKqVe8FKOUmgxsi9t8NvBi5PWLwDleHtMtftem/nKdlTaheNnmtL7/1fpdKffJjfGZDejUnFMOb5d0gWImdLvtfc5+cprr7108+BAf1NigaafwLx9/HbYEAOat2ub7vO7clVZUYTqT3DqXJIg9b1+4CM8e2q21D2q8xXEhJ6XUbiCMdRwFSqn1EQ3rIy6y0Khto8/6eTk8c+kxvh7DTQ6hwkNb8sylhYH1DHVqeNJxB54xwP08lVLKce2KOSv9DRUvKS2rXPMyN41j6br4r7ZTayoAisgYYAxAQUEBxcXFaf1OSUlJ0u+u2Olu2O5Ux0cTJ1E/V1i21vLZbtiwgeLiqgcplS43VFSotH7L7jte6gJg/y4WzJ6ekaYoTrRt2mQfonp8xzymrStLeIy2jYTN+5K3WomObadrz0HF9RPdZyIu27Wl2m+d06Me7yxP7vc//68fct3RNd1vdrq+/S6z9QePjcjnw5UHGb/CXlNsFODefftc32Ob9zp36SX6jdYNha37q1/LRV8upuGW5GnfPb/3IxQXF1NYkMvcjYnbmkzv+0zJBsOxUUQ6REYbHQDbcBOl1NPA0wCFhYWqqKgorYMVFxeT7LvNV2+HGc4btvjfeqbdRq54aW6N/cZ8vJeVD45m09w18OVC2rdvT1HREY51VTIhcZK9hXedysC7PmJIt1YUFR3r6rtQ8//iSFeK34yndZs2FBUVOv5usmM7OWevr50HGzfU2P7SNacy+L5P+O1pvSmycZvNHq447Nbk5VO3NetB347N6NO+WUpdm3eXwsRPkv6eHZecMogTerapfL+v9XreWf550u98vqmcHkcMpnPL/JS6vmI5fO2+bsZ7157A5pJSRvRux745qxm/YhHnHNmRdxYknrvbsk9RVFTE6q17+Wr9Ln77ny+YeONwlnw+M+F1XLNtL0ye5EhTjd+I3F+nDuxcI1X8atWa3xUdlfT3vL73Abq2zqeoqIjJu79i7sbEKY0yve8zxXFmNLH4iYjcGXl/iIgM9k9aJe9iJVok8u+4AI6ZkEyHxkcmWQOweXepr2VMmzWsx4QbhiV0Tf1+9OFJv79yyx6mfLO5MkX6pKWb+Hq7txOnQbseErkec3OEeXeMTDjX4sTV85s3vuD0R6dkIi8lsUYDoF0zZ5FyThNrpuua7d+pOSN6W15liSxLredg3RDAiQ9P4sp/z6OktIyXZ6zyRV8q/vdFesEpXnHmEdYKiMNSpMsPCzcpNZ8CjgUujrzfjRVt5Rki8howA+gtImtF5OfAg8BIEfkGGBl5HxpuIoTs1knkJclTfsx9n/DQhPSrou3anzo0sU/7ZjRpYD/QbJMiPPdgeQU/fXY2wx8uBuDyF+Zw/yxvVyP/9tTela9/d1rvJHt6gy5zHF7pcJp2xG30XVg8MSn52o5MglWO72FNQh/Zxb8FnW65NJIGKNpK5CZoL3rc9j7XvTY/GFE2uDEcQ5RSV2NFVKGU2g7U91KMUupipVQHpVQ9pVRnpdSzSqmtSqmTlVI9I//GR10Fipvn7e8/GVRjW4t8T09ZNc74W2a92zMGJM/DtPeA/2GZvds3rXx9ybGHJtnTQinFQxOWsnZ7epUKtQl28EhHg7xcR/uVOUxzU5GGgZly0wjb7X6c6kwMxytXDGXlg6MdBwr4zcoHR3N5XNG3RMrKKhTvhjgqcmM4DopILpHrLyJt8XgBYDbgdzhuJqzdnll6kPp5Ofz6lF4JP08nrNYpw3q24R4HqSDieXH6Sp4q/pYT/uTMzx1PJlezd0HT1DsFoCOWzi2dLWwr83HE0aVV9bkTN5kk3aYPcZMNOhFHJRhxLN9UwuJ16VUVXOIgLB7gzSuP5Z2rj+fJHx1dbXvPgqY0yMvhhiTPY5i4MRyPAW8D7UTkPmAqcL8vqjRGZ8ORzRzTtRU/PbZrtW1O2ra7/pdZuupMVvf26eDMcKTqtc9duY073vky6T5Ocdp7duqq8uJuPzRiSPp1bMbc35/Co/+XOFHkWy5r08xZkbkDomdBU35V1L3G9lMe+YzRj00FYOmGXY7ulXELvmfphl2OR+eFXVtxZJcWjI7LutykQR7L7h1FUW89q4c6iqoS626cDMwDTsbqQ5yjlFriozYtKT1Y5wZZgTDmxG41N4Zko53UKQE4oUcbxiWJEoqyZMMu+nVMnIDwh/9wXiluwg3DUk64v/izwVz+/Oykhreswtl97EU/aUi31nxw/TD6tG+KiNAiP/E8zL3j3TUpQUzVvDxjJXeMW8zfLjqSs4/sxK79B9m1z34+Md0MC4lwGlAQNE5XjivgHaXUUqXUk0qpJ+qi0QC48T9fONqvT/vUvdGHf5g495HXA5voRGAqDmltNZo3nNLTWwEpaFivpm8+OrpL1tD4Qc+CJo72++Ggzvz750NS7rd2+z7mrNzG1G+2ZCqtRmivHcN7teWSuNFbPM5HHN7ciId3aFY5GrK71unilQcg2UDtjnFWEsUnPrUm6kc9OiVt16hb6uflMO/3pwRyLDe4MWczRcTfJcZZwDYHPtji3xYx/rphKfdLFDEB7ofsqRCHjuZzjuzEK1cM4cLCLp4ePx0a1bcamHOP6lS57aEkxjYd7NodJ40zWG6h+HBYO16dtZoL/jGDnzw7i4279jNxdXrpw+NJVk8lUdGwKH7Ncbx5pc36oDgGd22V9u8vXLuDTburIvm86mClMrQA2yPFnb7fEWypgdZNGiQ0bOkEL3iBG8MxApghIt+KyEIRWSQi4Re/1Yy///hourZpnNQoREm2T35973plbhARju/RJmXDs3Ct/6VRG9bL5cs/nsYdo/sCMKJ3W9o2dbZOQSe+iDlXQ+6fyMtfHeDL79ObdI3lHzZRe1FSeTi2ljibhJ6z0vkcwsoHR1PowCjkOHg2EvGDJ6Yx/KHiyvdejTgKmjWkaYIw9SiHxE36z9sYXFWJRy48wnZ7t9uqFqFu3l2aVh2XdHBjOEYB3YGTgLOAMyP/GmJo5KLBT9Y4ex36+uMh7pIGFqRYSLYzgY83Fi/SwzdpkEdOjrDgzpH886eFnNAjdQ9fNwoPrdmYelG17uTDCxJ+1r1tcnfbbyMu1yc+/Yaut4znizU72FdWsxGettzb2tuZEI1Uis3Wm8pwOBkBRRmeYiI6vlF++avgysG2apz4eXxuqrXC/KQ/F3Pao5MD0eMmyWHyJZwGwN18rpNRiRc0aZDHqBRrNOJJFZ3jxPV1d4YRT7E4Wf/iJnlf0t/xeFb+kyUba2yL7TSIeD+nFd87TsRz01YCVqh1jxY5jNLPnQ5Y7qFRNuuUvtloXwTtsYuP4gdHdASgW5vGXDQ4teu1b8dmvLdwfcLP98R15mLniWav2Mb1Y/1bkJesqbj7va/o1rYxux1mA/ACx4YjmmokHqXU3d7JqVsEZDf44PrU8y1ucdI+b94dbC3o8/4+3XV1vinLa05Yj+rvzsimw/PTVjC8l9XD9SPC+ziHI7PYy7h8R/VIq7Jy5xGEfo8E732vZiekrLyC77bssd2/PCZq7NPfFjk6RkuXi3Njz85DE5ayfqd/Nd1zUzxwlz0/x7dj2+EmyWHsFWqI5aqqk5FVmTL+uhNYs21fIIajU4tGNRdkecBuB+lNgmb+avfzLnb13oNIQbEhg0bm9jMOd7RSvkmDPMc5qewY9pDzyKHBh6U/4e0Eu0p/nyxJXF739H7ujf+FhV0oK6+ojKJyg99T1JnMC/mB4zkOpdRfYv7uA4qATim+ZrChX8fmnN6/fSDlK8872p9LdOW/qzKwLlm/y7YBduo1mnnryV7JyhpSBR8k4xcnduOPZ/dPuZ+TI8TLmB2zoM7PHvTQbu4MjV34cLL5DTdzjVFyc6TGItR4Ln9+duXr2FOXziR9RxeVJTO5X/wgk9Ul+YDNqq26zWGtnWez7N8p8aIwtySaiA5iofuov03hHhtXwpptzsIWszFSKlMCaQeSHCNRHrAL/+l8MWImDO/lrh6b3YLFMJI4TEpQmTMdLQM7Ox/ZBjUf6hQ3cxyLqBqR5QJtgXv8EKUryVIO/PEH/Tjv6E6B1jSO8sqsVdz+tn3KiqCyv85fU7N627KNwYQGxvPl9zv59esL+EUffVf5R28lP7PUJmtqmjbMQynFFodhuV6TLEu0HTpm8y2J8Z65TV3z6hVDGOjCJerGcGxyUdwqXdzMcZwZ87oM2KiUCm4aXwMmJvGpVigVitEAEhqNIEm39/eXC47wvDd15uNWfqG383K5MMl+Wz2qfTKwc3NX5XGhqgf23kLnGU6P6drS1TGSRZhVKJi3yr5U65OTlvvew22cYs1EPF+4PL9B4/b2dxq8EMWNofW6Ro4dblxVv1JKrYr8fa+UKhORP/mmTEOSrRj1e9g8e0MZXW8Zzx6Xk51BDefTdb2cP6izt0JiSBUy/LMXa1ZiTIfWjd2nyo/2UN3kNkqW78qOZNekQinWJZjDePjDZTz4wVJXxzrCZUBBxxbO/fuJCLuWSkG+8O4X6/ho8QbfnzOndVYASv23G64Mx0ibbaO8EqIzi9bu5Iy/TalWHzkeP+6bgzHhkE8tsHrHa7bvRSkVeKhrNpLKln2xpmYUVjopTQ51Ma8VpbxCuU5dkSzFiFt4oDqkAAAgAElEQVSmLd/iaSGgaGixU/JyMkve9/XG3VzzaniFjKJc99p8xrw8r1oaFD/o0iqfy4/v6mhfmzgVz0l59UTkqsj8Ru9IqpHo3wpgkf8Sw+fS52fz1fpdvDwz2DWQj39as/rZ6Y9O4blpKznmvk/4dnMJG3clv2HD6JOd9fhUrnXQKN2dRv0Nvxl0qDt3ULqUVShbw5WM613WZkhmZr783lm9CL9wajf6d7LPGzbNZv1N0MQ+Wxt3+d+R+8NZzp4XXeY4XgU+AB4AbonZvjvsanxB4SSxYSZ1HRKxeqv94qbJX2+OfL6XJg3d+Yr9Iva/v+j7nSxykIsp3fUleTniKElfReTRHrfge4Yc1pr2zRuyp7SM1dv2cngH+wbJaUr1WNIJlWyRX89NfSOAhCV/E6FLZTs7nJ6zRI+VDmVxNu3VQIQNE1f7P/Wc8k5USu0EdgIXi0hLoCfWAkBEBKVUMMlRNCfdG/knQw/h3zNXu/rOZxHDoVCVOYe81gVQPzeHAw5XDy9e57wHe+uoPny5bhfHdXeW6j2e8dcNc5STZ8a6cvYdKK+cR+jaOp+VW62Fc1/f652XNZ32+dJju/oekhuU2WiURpp0p5O9CQ2H6yM6p0PzhtXWsLRqXN9R57Eu4djRKCJXYBVz+hD4Y+Tfu/yRZXv800VkmYgsF5FbUn8jWNJdi3DZcYel3ikBSsGqrclXEGcygdjBgwlMO7q0yufxi49yXB87nt4Oap1EiY3/Xxlzrrys5GjXBtZPkZ5WBN6Z72/N6JtO7+3r7wMU9W7Le9ed4Pp7TldCR69TfADCnya4m7x3Q7+O1UejxmjUxM0M1fXAMcAqpdQI4CjAfjWMx0RqnT+JNRnfF2v00zeIYzulQV56k33d2yaeWE2VqM5J25cqx00y+ruM4nFKCxcRIm6wcz8lKnS3IkGOo3Swc7uk8uHvKS2vHDn6xf8dcwj/u+aEhIv9vKCoV9uUmXjtcOveiw/ftctU4BV+rxlJNG+TTbhp7fYrpfYDiEgDpdRSwP8ujcVgYLlS6jul1AFgLHB2QMcOjaHdkrtytjhYhxBfy9gND18wkLeuOpaVD45O+zfscBvDbsdZkcynsRzRuaahO5jAcpSVe9c42K0A/tclhUm/c9vbizj5cOerp9OtzzKgc3OuP9m/ao4n9Umc2j0ZTpclRDtHbkbObkOD40m0OtwrdJifyRQ3hmOtiLQA3gE+FpFxgL9j7So6AWtitVBL8mQlncBM8XBF6xMko3WSPP6pyK+fxyCbWhLJCKo62i+G1XTxjZ2zpkaY8kSblOZg1QH3CjvjfFz31Maxb0fnPc9MKvZlGvqaiMcvPopDWqcX4OB0xBF1VblpbE/rl54xi3JsXIetscdF1TRcBO8aR2EaYrVu1ymldgB3icgkoDkwwU9xsRJstlU7/SIyBhgDUFBQQHFxcVoHKikpSeu7Xy5eTKOty9I6ZiIWLPiCA2sS37Tzvk5dXnbp/Jn45w2uTnFxMd/tTL36KN1rE8uCTfaRI8fc90m195/OtU/gfNOb9sUrJ0+eTP001kv0bJFDuYLvdlojnMmfFaf8zqJl3zn+/RxVkfZ5O+Dh6AqgWX24sHd9mmxbRnHx12n9xupdzlap7dmzh5KSCpqK88b7u2+/o5j0Sy/vjgtxj6/DkSkd8vZ68gwkon6O8vX3waHhUEopEXkHGBR5/5mvqmqyFoitxNKZuNGOUupp4GmAwsJCVVRUlNaBiouLqfHdCeNTfq9fv34UuSyWlOr3H5qznx8O6syfLzjCdp+l21P7edM9DzVwcA669CukVWkZzJjmu6ZNc9fA56krF3fs3BlWrXT8uyeeeCIN04gSiv6Xut5inacRI0bAh4nPWW6O8MEK56npmzRqkPZ5U0rBx++n3tEh/7psKENSuFFTUV6huHN6ak2N8vNp0gRGDerEko+dGanDunWjqKhH2tq6D9zrKqW8Wy44caDnbUUszRrkePfcJ8DNGHamiBzjm5LkzAF6ishhIlIfuAh4NyQtgfLmvMQ9Jz/qbGTC+h37HVUG9IJUUUtR3PqT0zEa6fCrou6O9jtjQHsgs2y6Xq7neOWKIRkbDXCetC96+b52kTAz0zVVfj9XfruqPEwwkBA3hmMElvH4NrJyfJGIpO7yeUAkmeI1WCHAS4A3lFLuq614yMe/PjHMwwP6ZQx1EuIabQgzpb7DKDY3YbdOjZEXONUVrVfuVVvwk6Huas/H07BecOcIqgx/spKuib5TV+nUxP9r5GYpaqh5qZRS7wPejbczpFtcCGIYN+v2vcnjy6fcNCIgJRZOTsHfLjrKk2PVc9jIuzGuQax7iOL0fokmt+uTYKW7U07o0Yapy7dwWr/2rhecVse/7myjernsO1h9PkEphZPM7w3r5bD/oOW67drGfe6wIPE7OeMvBvpf38aNaVoNDAMuVUqtwmonMgtf0BilFOMWfJ8wXjz+8clvkLmLo6/LxqEiRaMYtCtr/MJ1KV0qThv8VDhN+OdmxHHFsODqkjmxZ8d1b825R3Xi9TFDefJHmRlcrxorr1e7x5actQs5VsDry1JbjvvPHQDAb0/txZkZhKBHGXe1u9r1bvC7k9koz39flZun+CngWODiyPvdWIvyahWXTdjDPe99xYeLN3L92AU8/uk3tvvFPkAPnDeAIpfZQe1oke9uYVx5kjvwd6cF13uO8sbc9CNZ3OLUAC1YE14dh09+U9OdeVVkbsOJH/7K4d3JyRGGdGvtWa2X+DkotyVc3VS4TMWUm0bwwuVV06Z2i2hXbd3LQQdW9rA2jVn54GiuOamnJ3M6R3Rp4XitiVs6t/SuZPTFgzNzPaaLG8MxRCl1NbAfQCm1HXBfiCALeHbqiko3kJNFdhcPPiSjm/WaEVYESKHLzKzJ2h6nk6/pUNQ7cyOZKU4Nh5O1LgBL7j49Ezm2FDSrmbLl5tP7UD+3ulvvsYuPct2AuyXqxonH7QiwZRq1RxLRpVU++fWrvOW5CUaRO0tTGw4/anK7WWfjhm5t3K+0t+PJHx3Nfeekrj3vB27umoOR1B8KQETaAvrW5syQqG985Rb7XFBeRqo0a2Q9PA1cRvQEOa8SmyvohcsHc1jIfuSmHmcFbuTxIq9k5FDdzXh4+6aMHXNsjf0ae+D+jBKt9rcyLuPykRmusvaSwV3to7WcZKH1I2GkH8bIS3JzxHHOL69xYzgeA94GCkTkPmAqcL8vqjQg6huf8d1W34916XFd+c3IXvz8BGs19JDDqvc+v0kQipjovh7arZXnKbWLf1dU7X3vAvtEg9O/DaZOQjr5kYIm0TUQqT7HkehS+ZEWfUdcQIXbqoJ+culx9jm1tu5PbTjSTZiZDL86Zg08ikyLX+EeJI7/B0qpV4CbsIzFOuAcpdR//BIWNk6ica49qQf/u8Z9ZtB4GuTlct3JPWlYL5eVD46uUbDl759960qjXe81U+J97OcebZ/x5f73g1mnHlJHyxWJJArVJ+2DrJsRn7pEp/M4sHMLLjuuKw+eN8D1dw9NM/VJMrzIoHz1iCqX8XvXnsCto/p4tlaouc2c6NJ7vHe52uEmrXpD4AzgFOAk4PTItlqJk6iXG0/tzQCbxHqZEt+OJFpUl0n+okw5rV97fjwknIk50LdI0YQbhvFwpPxsoozJItUnxxP9T3q2835UVVGheOj8qvK4urlj7vpBPy5KMuHbq8D+nPjx3/jZ8emXPIjyw0FVCS/6d2rOL4dnNvd43UnWfOgbv7TvHAa1gNXNmOkloB+Wy+oJ4HDgZT9E6UCYj1P8w7zZwQR9EMQ/tGGus4rtKX94Q/iLMaP0ad+MCwqtxiIvwcSzNeKoem/XeM++/WTPIqliOVCuuPCYqsbMyQpur7Mjx/PcZYW8ddVxCT/v2aLqPJ59ZHC5TYdkELBwVVF3Zt12Mh2aW33rh9OoZW/Hr0f2Yvl9o6qFMYeBmxnG3kqpI2LeTxKR5OXnsozYXqDTdQJ+oFknsJJxV59QY4FWWMSOONwUdgqa5y8/hsufnwPAOUdaqeBFqq+rsLve7Zr6M5ivsS5Jg3vNLjX7+OtOYPRjU4HqnYT+nexH+H6kuknXU3XjyF5cG5PK3kvDKyLkhdg2RXEz4pgvIkOjb0RkCJA8m12WEb1RRCDXp1TUTojvBE72ueCPUxrVz6WVh+GYdYERvatqbtw8qg9gNXJ2I47oAtB0q0k64WBcKeC125JXkGzno5ZkxE7ahzUPk67huNbH+ifJuDVyfwWBq3UcwHQRWSkiK4EZwPAgc1b5TXRBXY6I45rI/hB+j8IJ8SU2DcmJGohdBxSfLtlU4/Po+phLfazY16ZJdUNwMEXKdafJCP0kLAlelhcOgkznT9zgxlUVzHR9iERvlPIKxc59zlNee40Gz6ojfjT4EP7ywWK2OQiXBCsXUV0m1k20IabmQ9RVFe2s+BH0MLRbK2Z+t43CrtUXmU5cal/oKkrU2HVs3rBGfragWLy16rwF2Zi3zHc/up5z+yk+KNEPx4Yjkp+qVhN7T973vn0BoCDQNWIoHhGhS9Mctu13Nu/xv2v9y/+TzUQb5wuP6cIbc9dy/tGdPT9G/cg6hwNxrqpU7XDUYzv91pM915QODRIEHPhhUBKtZAfLhbdpd82gFT/djDrhJhy3UETeFpHPg06rHhS6pCnP5CEIOu21G3q003MSe0RAKVTiG+0o0X5C55b5zLztZF+SU0arGsbXWk91z+sWrnts99bcY5NmI1HocyaE667WGzdn+xXgeeB84CzgzMi/tYZkSQODZMfe9N1kfoRw1nZaNwmml5jI/RlE8atoDfT4BHspRxyaGQ4R4adDa84B+TFKT7Qm4q6z+oYaiq4DbgzHZqXUu0qpFUqpVdE/35SFgNIk81Ymz0CQxYgg3LUctYUgAvguP74r0245icPjUven6iyFaTf+r7BL6p185saRvYDqJQ8uPa5rSGr0wc0t+wcReUZELhaR86J/vikLAV1GHHr18WoHRyRJ5hfUZU+cgsT/Ky4idGpRM513KrdomCMOrxNZpsOIPlY49fmDquadRKTOVxl0YzguB47Eiq46iyp3Va1Bl/C7TIbdmdZbdkvDkAOlEkVqDY+rj6LBmqmE1zXMxaapQqrDdPPrEArcv1NzVj44ujIBaRVVz1mzkA3c6AGZF65yi5v/8RFKKffZx7IIp4bjiQyrsaUi/MfFOZf0bcDy3WVs22Nfpe2oQ1owf/UOX479/GXH0CNBPqf4NkcHX30iBWE2kH3apzIc4WnTwXAkIrapOOuIjrwyK5NyvJnx5I+PDryinpsRx0wR6euHCBG5QEQWi0iFiBTGfXariCwXkWUicpofx49S4XCO48yBHf2UoUUj55Qm9YWrI4Wo7KioUPznymMZc6L3ZVlH9GmXMAIp/hzqcEoTpdP2qpxuOgxIkMJDB7LlORjksgBbbcDNiOME4FIRWQGUYnWglFLKi+xdXwLnAf+M3RgxVBdhJVfsCHwiIr2UUr4kTNLFVRVitpO0SOYeq1BwTNdWHNM12KRs8QVukl1ar+pxJ6Jt0wZs3l2a0K0WRthn/dwcDpRX0KdD8hDpMBtvjQcc1e4YP8KndUeLleNKqSVg6wM+GxirlCoFVojIcmAwVroTz9FlHUcQk6VekszgBmWM6+flVFuZHX8Gw7yyqdbW6OaSef6yY7j8BSsxY6hPhAOj1bllI9Zu3xeAmOrEdpY0u3yB4KZvuxoYBlwaCcNVQM20lt7SCVgT835tZJsvhFnfIpZsG3EkO21BDeLiqybGN8bJRkXXJHG1eUFl8syIOTu3R9VamwsGddYqU0DhoS0rSxkDnJ+gYFcQOGmQJ/22iK/vHeW/mDiid9O5R3Uiu2YlvcHNiOMprBrjJwF3A7uBt4BjnHxZRD4B2tt8dLtSalyir9lss20BRGQMMAagoKCA4uJiJ7Kq8X2Js0mOdH7blY7d6S8o2V9a6ru+WEpKSlClKxN+vrukJBA9Fx2imPJN1fstW6pnFN6xc1fC765ePBc/pzb377fyUs2aNZPv8nPIKT9A9NbesWUDxcXbfTy6PUe2FWZvgKmTJ9OrZQ5fb7fuuct7lDJ//vzK/UrWf0dxcTgTv8321PRIx99L06ZMDkgN3HN8IxrXszQcPGgt5uzfYCsL5leVlw7y2UtESQDPnBvDMUQpdbSIzAdQSm0XEcdZwJRS6WT/WgvErgLqjFW21u73nwaeBigsLFRFRUWuD7Zk/S6YOiXlfun8thuWb9oN09J7IBo0aOC7vliKi4u56swi/jRnvO3njfLzKSoaHoiWGyd/wP6DVgO4+WADoCpleNNmzWCnfXSX3+dr9M4veXHGKkYWnUCL/Pq8tuRDoAyALp27UFTkS8xJUo49oZytJQfo2KIRf182HbZbxuvUk4azeN0umGlVTBgwYABFh/vtWLBn0P6D3D3jo2rbKq/VhPHV3wdM3mcfwcGDFA4aZG2InK+w9MRSXFzsuw43TpGDIpJLpMcvIm2xRiB+8i5wkYg0EJHDgJ7AbL8OpsscR4fmNRdqZSu6uP/CXLF1x5l9mXXbybSIZFv9cFVZ5Wd+T8wnokFeLh1tFgTmiGjjs9c5qkrFlGD4ZlNJyGqCx43heAx4G2gnIvcBU4H7vRAhIueKyFrgWGC8iHwIoJRaDLwBfAVMAK72K6IKaha5CYvGDdwvKBoccNSSU+KT6vlJbFBB/LxBmFkB8nJzKGhmX9FPh0CI2DOTI/o02LrosCN6zkRgvyZVMYMkpeEQkTwApdQrwE3AA8B64Byl1H+8EKGUelsp1Vkp1UApVaCUOi3ms/uUUt2VUr2VUh94cbxE6DLiSIdHLzoybAm2BHlOY9uZFVv2VPvs1L5202twaOu6F0oZT2zggIhUO49hRqhrbDfYvb9q1KizgfMLJyOOSteQUmqpUupJpdQT0RDa2kQ2GY74PD7tmzXkx0MO4bnLHMUqBMaVw71f+JeInyapnHd8j9a22/uEXK9chzYn/q7XpSHURUcyNu8u1ca1FyRODEedOS3a+ONtOLNb9XTp8RclJ0e479wBNbKfhk2j+sHl8bnl9OBqLmfCwLZVCwF1eLjiRxWxDXaYT0Q2NMi5OZIVBs5rnDzVbUXkN4k+VEo94qGeUNFljsMJ8SujdaXcaR4XD0i2HiKRyyWMdB8LN1f5xFdv25tkz3Co7qoKz3RkQ4Oclyt6WP+AcfLU5AJNgKYJ/moNuruqPvnNcC46pgvDerap9lBdWOh9qVGv0GUUdzDBJH37BJPWftK9edVjN+PbrUn2DIb4M3NYm8aVr5s3Cq8wWBbYDfKybbWuRzgZcaxXSt3tuxINGLfAdomINvRo14QHz7dSgxXe+3Hl9l4F+trvCk0Mx7Hd7ec4wqBdvvDtTuv17tKy5DsHQPyool5uDsvuPZ3iZZsZ0i2886bTivpE5OmQrz8EzBxHDNOWb7Hd/uylhbbbwyQbHirQY8Rx2xmJ5z6STajXFU6xWeDXIC+X0/rZR6IZqsgVHQKqg8fJiONk31VoQqJGLszhepQak+GahEymQnf336GtG6feyWs0a2muGdGDCwu7MHeWL7lDazW6JagMipSGQym1LQghOlCWYHJcx6Zv467SytejBujbM9TBcOhsWHUgJ0do37whTepnTyM45sRuWtQS2by7tE7eX+EX9dWIgxo0cunQuaW+i9jy01gFbzCk4rYzDg/1+M0a5rFrfxn59XNp26wBAC3yw/dMBEXdDAlIgJPecVgrjbPNpA2OpDlvle84D6ZvaHfutBNkcEufyHopRVXYsA4joKAwhiOGRBc+dij62e9GBKQmsYZsoGGk2l1YSfwuHtwl9U4hkWWX0mBD1KmnlHZTVoFg/AgxnNynHQvW1Ey93SK/Hu9fN4xd+w+GoMoi2xqbaJRi0N6/ow9pweerd9AxJsOwbkY3Sz2iBhtiO0a63Wd+YgxHDE0a1jwd//75EK3XSehKdPge9Mrjod1a8/nqHVovHmuaRZPQBnsq7y8F/SOeiiuGHRaeoIAxrqoY7Nq4od3CTVf+4s8GAzCsU3bZ+GhKlKB7YVXprms2zl/fO4qXfz44WEE2nN617kyi1laiqzcU0KpxfVY+OJqi3u3CFRUgxnDEUGHTyoWdL2d4r7asfHA0HZpk16X63Wm9Oa57a07pG2z1ODtDFXUn1M/LoVG93Jo7BEzb/Oy6loaaRDON1CX3VCzZ1Y31GTvDobPLQ2cOaZXPq78YGvhxo0YiUU2J2O3180wDbkiPqhFH3bQc5smJYfTAjjW2ZUtqD90I67RdeWJ3Rg/swE+G2qcSadqwyk208A+nBiXLkCZFvduGLcGW6P1dV0ccxnDE0KlFIy7uE/66g9pAWBl8Wjauz5M/OppmDe3nEQ5pVbUOp6EGbitDdjI0kvyxY4vgsyvrgDEcBl/QYaA2vLO+ntgjOluROMeGmH02G9C1R3/V8O48dGIjerSrmxGXWhgOEXlYRJaKyEIReVtEWsR8dquILBeRZSJyWrLf8USL3weoI+hwHqNhr2EWI0rEb07tDdTdtNzZTk6O0K4OBzno0iX7GLhVKVUmIn8CbgVuFpG+wEVAP6Aj8ImI9FJKlSf5rVrBA+cNIL9+9rpSwo5GMxgM/qGF4VBKfRTzdibww8jrs4GxSqlSYIWILAcGA7U+//PFgw8JW0JG6Go3ommww3YRaXp6tCbbn4nahBaGI46fAa9HXnfCMiRR1ka2+YZ5oL1B12i0erk5jL/uhHDqcNigoRdNS56/7BhG9Kk7C+x0JzDDISKfAHaFI25XSo2L7HM7UAa8Ev2azf62j5qIjAHGABQUFFBcXJyWztIDpdUOm+7veE1JSUlCLWFqTKRLh/N28OABQPhuxQqKi7+v9tnmcCQB1jlbufALALZv36bFuYLk91hYbNu2H4BFixYiG/Tq5+p4viAYXYFdCaXUKck+F5FLgTOBk1XVbOZaIDbNaWfAtjC4Uupp4GmAwsJCVVRUlJbOj1d9DByofJ/u73hNcXFxdS0Txle+DFNjpa4YPaDHeXvrm4+Ag3TtehhFRT3DllNJcXExPzv1BD7eMJP7zxtAv456pOOucY9pwPPfzYYtmxkwcKB2KT10PF8QjC4twgJE5HTgZuAHSqm9MR+9C1wkIg1E5DCgJzDbVy1+/ngtpnEWT+SHQX79PMZdc4I2RsNgcIMuY78ngAbAxxHf+Eyl1JVKqcUi8gbwFZYL6+q6EFGVjdTPy2HPAXNpDN5ipoD0RAvDoZTqkeSz+4D7gtKi6Zyu9tTLrRq8njmwQ4hKamImoLOXHXstt7GOa3HqMlq4qnSibSNjOdLh3KOqgt3uObt/iEqqqCqZYBqdbGXh2p0ATFu+NWQlhliM4YhjYNs8/vp/R4QtI+s49+gqw6HLqE0TGQYPOFheEbYEQwzGcNhwTNdwize5YbAmWls11jc5pPFyZD92JQ8M4WEMhw3R1cXZQNumDcKWAECrfH0NhyH7MXZDL4zhsCFXF1+LA07s1SZsCUD1leJhpVRPhGlzsh9zDfXCGA4bGmbRmgRdUnvooaI6mpwagweYEYdeGMNhQ6IiQDrRuWUjQJ8GW+tG2rQ6WU+rxvo/k3UJYziylOikuC7py3UZ+cSinyKDW/p1bAbAKYcXhKzEEIsxHFnKlUXd6dKqkZ4ZQ02LbfCIvMjCUjNm1AtjOBJweIdmYUtISq+Cpky56SStw2ANhkypXMRpLIdWaJFyREfeufo4ysrN3ZoOuqWH0EuNwZD9GMORgAZ5uTQwZyctyir0aKqj0y6a2TGDIesxriqD55RrYjgMBoM/GMNh8Jz8LFoHYzAY3GMMh8FzmmqyDiYv4qrKphQyBkM2YLz4hlrLKYfWo3lBZ345vFvYUgxp0qNdExas2UGzhqap0glzNQy1lvq5wu2j+4Ytw5AB957Tn66ymZ4FTcOWYojBuKoMBoO2NKyXS/82pn+rG8ZwGAwGg8EVWhgOEblHRBaKyAIR+UhEOka2i4g8JiLLI58fHbZWg8FgqOtoYTiAh5VSA5VSRwLvAXdGto8Cekb+xgB/D0mfwWAwGCJoYTiUUrti3jamKkvE2cBLymIm0EJEOgQu0GAwGAyVaDPrJCL3AZcAO4ERkc2dgDUxu62NbFsfrDqDwWAwRJGgEtKJyCdAe5uPbldKjYvZ71agoVLqDyIyHnhAKTU18tlE4Cal1Dyb3x+D5c6ioKBg0NixY9PSWVJSQpMmTdL6rp9kg67LJuwB4IXTG4cpqZJsOGc6YXS5ozbqGjFixDylVGHKHZVSWv0BhwJfRl7/E7g45rNlQIdUvzFo0CCVLpMmTUr7u36SDboOvfk9dejN74UnJo5sOGc6YXS5ozbqAuYqB+20FnMcItIz5u0PgKWR1+8Cl0Siq4YCO5VSxk1lMBgMIaLLHMeDItIbqABWAVdGtr8PnAEsB/YCl4cjz2AwGAxRtDAcSqnzE2xXwNUByzGkyfw7RnKgvCJsGQaDwWe0MByG2kFLU8bWYKgTaDHHYTAYDIbswRgOg8FgMLjCGA6DwWAwuMIYDoPBYDC4whgOg8FgMLjCGA6DwWAwuCKwXFVBIiKbsRYSpkMbYIuHcrzC6HKPrtqMLncYXe7IRNehSqm2qXaqlYYjE0RkrnKS5CtgjC736KrN6HKH0eWOIHQZV5XBYDAYXGEMh8FgMBhcYQxHTZ4OW0ACjC736KrN6HKH0eUO33WZOQ6DwWAwuMKMOAwGg8HgCmM4DAaDIcsQEQnz+HXScIhIbuTfUE++HSKi3TXR8TxB1XXUDRFpHvlXq2spIu0j/2p1PUWkn4g0DFuHHSJyvIh0D1uHDY3CPLhWN7bfRG6CF4Hfi0grpckEj4gMFpHrAJRS2lRCEpEhIvIv4GYRSbkoKChEpFBEXgbu1OWhFpEcEWkmImFWRj4AAAdHSURBVO8Bj4E+11JEjhKRicA9UFkgLXREZKCITAXuBVqHrScWETlaRD4CPgWah60niogMFZG3gCdF5NSwOk91xnCISDfgKWAScChwj4iMDlcViMgNwNtYxmxUZFuoPWkRyRWRB7CiM6YBRwN/EJGCkHXliMgTwD+BiUAH4C4RyQ9TF1Qaid1APaCTiPwfhDvqEIu/Ai8BLyqlfhGWlgT8HnhTKXWuUup7CH80JCL1ROSfWPf+Y8CHQFHks1DbSxEpwmrD/gssA34CtAxDS50xHMAgYIlS6gXgRmABcKaIdAlVlVVP/UzgKuBWAKVUecgPUA6wGrggcr5uAIYS8vA40jh/Cpwc0fUQoICyMHXF0Acr1cOjwI9FpKlSqiKsaxkZWTQB5iulXgIQke4aNIA5kZFiiVLq0ci2kSLSAgjbjdwA+AwYppR6D3gLOFxE8jQYQQ4A5iilXgFexuqklIQhpNYajsiQrlfMpjlAZxHpopTajtWT3gGcG7Ku8cDCyL8lUZcVkQcoJF0VwGtKqa9FpIFSah2wFisHTqDEny+l1H+VUjtEZCQwF2vUcb+IHB6WrphGbjlwAFgR+btURA4J0jVkc3/dCAwRkTtEZBrwMPCCiAwKSlO8rkgDvAkYJiKjReQd4LdYPfzfRfYJ65ztUUq9qpTaF3mfB5QrpcqCNrg213IKcIGI3Al8jnXvPyUiFwSpC2qh4RCRFiIyHvgYuFBEmkQ+2g9MBS6MvF8GfAW0DmJizkZX4+hHSqlypdR+4C/Az0WkjVIqkF603fmK6NkBoJQqFZGmwGHAuiA0JdDVOLI92khvB36klBoJ7MVqpH13pdnpimnkCoFdSqnFwGLgD8DfI+4PX5+1ROdLKbULeBI4H2tEezGwHjg/iHmrJLp2A89jzbs8p5Q6DXgGGCoiQ/3WlUibUkpFXHzR6/UZcK6ItAxqxJGoDVNKLQBOB7oCv1JKFWF1gE8PuuNU6wwH0BjLL3lt5PWJke2bgZnAABEZrJQqB74Hjo802qHoirsZiyMarwVr0jwEXcNs9hkCLFZKrRORJiLSMwRd0fOlIv/OVUq9H9n3feAoLAMSiq4Iq4GmIvI6cBMwD/haKXUwgEYnoS6l1GPACKXUZKVUKfAOlpEL+3y9h9UIRv30c4GNQGkAuhJqUxYVEeOxMrLP8IA02emqfCaVUrOBthFdYLlumwJ7AtRXOwyHiFwiIsNFpFlkku1p4A2sUcZgEekUMRQzgfnAXyNWvB+w2q/J1RS6hohIx8h+AtbcBlaEyc0ishM42g9frwtdeZGvtADWiMjlWC6/I73W5EaXDYOwetG+jNJc6GqJ9VBvwDJkVwG9/eoNujlfEfdslEFYrsfykHR1imhaiOWaukZE2mBN9vYHtvqhy6G2ymcyYuyj3oj90e0h62oATAeujnz1ZKyItCA6v1V6NYnMc03kArYHXsXyyX+LZZ2vV0ptiexzPJZraq5S6uWY7z4CdMaKrrpEKbUsJF1zlFL/jmzLAbphDd8PADcopRaFrSuy/WXgx8CLwF8jD3youkSkGdZI6H6shvpGpdTXIemqvL8ibsbo502A+kqpbSHpij1fDYBjgT9jGVktzldk+2+w7v2ewK+VUl95pSsNbbHnLDcSqPIy8K1S6q4QdcXeY/2w3KDtgYPANUqpJV5qS4lSKuv+gNzIv72Af0de5wGPA/+N2/fXWL345kDT6PejrzXRlR/Z1g7LpaCDrmZAk8i2i4AfaqKrOdAwsu0U4GyNdDWOub9yNNLVKLLtOOAcjXQ1jdlez2tdGWrL91NbmrpaxFzLRkA3P86ZI/1hHTjNk52H1cP8E5bP8Sys+PTo54LVmxoes60JVnjkbCz/aUdNdXXWTNeciK4OmurS9TrqpqvOnS+dtXl0LTv5cc7c/GXNHIeIDMeabGyJFfZ4D9YwbUR0EllZZ/lu4K6Yr44GfgV8AQxQVmipjrrWaqZrQUTXek116XodddNVp86Xzto8vJbfe6krLcK2XC4s9TDgpzHvn8KafLwMmBfZloPl93sD6BrZdjZwotFldBldtV+Xztp01ZXW/yVsAS5Oej7Wqs6ob/DHwAOR1wuAayOvC7EWrxldRpfRVcd06axNV13p/GWNq0optVcpVaqskFWAkVhrMwAux0oL8B7wGtaqykDSFhhdRpfRpY8unbXpqisd8lLvohdiJQBUQAHwbmTzbuA2rBjwFSriA1QR8210GV1GV93SpbM2XXW5IWtGHDFUYCX32gIMjFjoO4AKpdRUFd7EkdFldBld+ujSWZuuupwTtq8snT+sTK0VWLmnfh62HqPL6DK69NOlszZddTn9y8qV4yLSGfgp8Iiycu9ogdHlDqPLHUaXe3TVpqsup2Sl4TAYDAZDeGTjHIfBYDAYQsQYDoPBYDC4whgOg8FgMLjCGA6DwWAwuMIYDoPBYDC4whgOg8FgMLjCGA6DwWAwuMIYDoPBYDC44v8BEIxjmogk60EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13167bbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure()\n",
    "plt.plot(d1['time'],d1['2_metre_temperature_surface']-273.15)\n",
    "plt.ylabel('Temperature $C \\degree$')\n",
    "plt.grid()\n",
    "fig.autofmt_xdate()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wind energy example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a more advanced use case, let's try to analyse potential wind energy production in three different locations. \n",
    "For this, we take wind speed at 100 m height, apply a simplified energy production curve to hourly data, and get a statistics for hourly and weekly average (wa) data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4HOW59/Hvb9cV3G3ZcbcxNphqwLHpvRgCB0JvoYQEkguScJJwQnJOXkjhJCGUc9LIC8GhF7+UYEI1hN7cuzGWjbFluci94Sbd7x/zyF6EtDuStRqtdH+uS5d2n31m5p5de2/N00ZmhnPOORdXKukAnHPOFRZPHM4552rFE4dzzrla8cThnHOuVjxxOOecqxVPHM4552rFE0czJGmAJJPUoh72tVDSyeHxzyT9bfcj/ML+j5dUUp/7rOYYD0j6dZbXN0raK58x5EOu89qN/V4m6dX63q8rHJ44mqjML/SGYmb/bWbfashjhgS4dz6PYWbtzGxBPo/RWFX3R4aZPWpmpyYZl0uWJw7XrElKJx1DHPVxddhc+HuVf544mgFJaUl3SFopaQHwtSqvd5R0v6SlkpZI+nXlF6qkQZL+JWlV2P5RSZ1qOM6tkh4Jj/8Umngqf3ZIujW81kvS05LKJH0q6fsZ+2gbmljWSJoNfDXLeb0dHk4Lx7hI0lWS3q1Sb+dVSdj3PZJelLQJOCFU6yZpnKQNkt6S1D/L9n+W9EKo+5GkQRl19w37WS1prqQLs8T/pqTfSBovaZ2k5yR1Ca9V/qV/jaRFwL9C+b9JmiVpbdh+aMb+DpE0OcT1JNAm47Vc70tbSXdK+izE8q6ktkDle7w2vMdHVN2XpCMlTQjbTZB0ZJVz/JWk90Jcr0rqluU9OVvSVEnrJc2XNCqUf+EKusq/tS+9V5JelnRDlX1Pk3RubT8n92WeOJqHbwNnAocAw4Hzq7z+ILAD2DvUORWobHIS8BugFzAU6AvcmuuAZnZDaOJpBxwNrAGek5QCngemAb2Bk4AbJZ0WNr0FGBR+TgOuzHKMY8PDg8OxnswVV3ApcBvQHqj8ArwM+BXQDZgKPJpl+0uAXwCdgeKwLyTtCYwDHgO6h3p/kbR/ln1dAXyT6P3dAfyhyuvHEb3vp0kaAjwO3AgUAS8Cz0tqJakV8A/gYaAL8P+A87K9CVXcARwGHBm2/w+gAqh8jzuF9/iDzI1ConshxN0VuAt4QVLXjGqXAlcTvSetgB9XF4CkEcBDwE1Ap3DshbU4h53vFdFncEnGvvcD+ofY6vI5uQyeOJqHC4H/MbPFZraaKBEAIKkHcDpwo5ltMrMVwN3AxQBmVmxm48xsq5mVEX0xHBf3wJKKiL7QvmdmU4iuIIrM7Jdmti30HdxXebwQ621mttrMFvPlL9L68JyZvWdmFWa2JZS9YGZvm9lW4D+BIyT1rWH7Z8xsvJntIEoww0L5mcBCM/u7me0ws8nA03w5UWd62Mxmmtkm4OfAhfpi89mt4XP5HLgoxDnOzLYTfdm3JfqyPxxoSfQ5bzezp4AJcd6MkMy/CfzAzJaYWbmZvR/ei1y+Bswzs4fDOT8OfAyclVHn72b2STiHMex6v6q6Bhgdzq8ixPJxnHMIMt+rZ4FhGVeOlxF9blup2+fkMnhbYPPQC1ic8fyzjMf9ib5wlkqqLEtV1pfUnejL+xiiv9BTRFcPOUlqCTwFPGZmT2Qcr5ektRlV08A7MWKtL4uzlZnZRkmrq4ml0rKMx5uBduFxf2BklXNrQXQVECeWz4g+i241vN6LjPfDzCokLSa6cisHltgXVy2N+951I2rWmh+zfqYvxJRx3N4Zz2t6v6rqS3QVVVeZn+EGSS8Q/UHyu/D72vByXT4nl8ETR/OwlOg/ZaV+GY8XA1uBbuEv6Kp+AxhwkJmtknQO8KeYx/0jsAH4ryrH+9TMBueIdVY1scaxCdij8omkr1RTp7oloXe+P5LaETXXlNby2IuBt8zslFpsU/Vz2Q6szCjPjLUUODAjToV6S0K93pKUkTz6sSsZZHtfVgJbiJoHp1WJL9fy2aVEX8SZ+gEv59iuOotDDNX5QvxAnM/1ceAWRX1hbYE3Mo5T28/JZfCmquZhDPB9SX0kdQZurnzBzJYCrwJ3SuogKaWoQ7yyOao9sJGoc7Q3UftzTpKuI2rSutTMKjJeGg+sl/ST0CGblnSApMpO8DHATyV1ltQH+F6OQy0HMudYTAP2lzRMUhti9McEZ0g6OvQV/Ar4KDSV1cY/gSGSviGpZfj5amYHdjUul7SfpD2AXwJPmVl5DXXHAF+TdFK4mvsRUdJ/H/iAqI/k+5JahE7gERnb1vi+hM9nNHCXooEL6dAJ3hooI+rrqGkey4vhnC8Nx70I2C+8F7V1P3B1OL+UpN6S9g2vTQUuDu9pdf10NcXWn+h9fTLj32FdPieXwRNH83Af8ArRl8dk4Jkqr19B1Gk5m6gZ6imgZ3jtF8ChwDqiTtCq29bkEqIvm1LtGln1s/CleBZRO/enRH/t/g3omHG8z8Jrr5K7+eBW4EFFo4wuNLNPiL4oXgPmsavzO5fHiDrmVxN1El8Wc7udzGwD0cCCi4n+El9G1EzSOstmDwMPhLptgO/XVNHM5gKXE13JrSR6H88KfUXbgHOBq4g+w4vI+KxivC8/BmYQ9YusDnGnzGwzUef/e+E9PrxKTKuI+gx+BKwi6lQ/08xWZjnnms5vPFEn+t1E/97eYtfVzM+JrkbWEP0beSzG/rYSvQcnZ9av4+fkMshv5ORcMiS9CTxiZvU62965fPMrDuecc7XiicM551yteFOVc865WvErDuecc7WSt3kcYcjf20QjFVoQDTO8RdIDRMM014WqV5nZ1DAm/X+BM4gmCV0VZnQi6Up2zQX4tZk9mO3Y3bp1swEDBtTzGTnnXNM2adKklWZWlKtePicAbgVODLNwWwLvSnopvHZTWBIh0+nA4PAzEriHaHZnF6JhksOJJvhMkjTWzGqcvTxgwAAmTpxYz6fjnHNNm6RYqw3kranKIhvD05bhJ1uHytnAQ2G7D4FOknoSLVg2LqxdtIZocbJR+YrbOedcdnnt4wgzUKcCK4i+/D8KL90mabqku8PsVIjWtsmcqVsSymoqr3qsayVNlDSxrKys3s/FOedcJK+JI6yyOQzoA4yQdADwU2BfolVSuwA/CdVV3S6ylFc91r1mNtzMhhcV5Wyic845V0cNMqrKzNYCbwKjzGxpaI7aCvydXevplPDFBd/6EC0HUFO5c865BOQtcUgqUrhTnKI7iZ0MfBz6LSpX9jwHmBk2GQtcocjhwLqwAN8rwKlh0bvORGvMvJKvuJ1zzmWXz1FVPYkWn0sTJagxZvZPRbchLSJqgpoKfCfUf5FoKG4x0XDcqwHMbLWkX7HrpjS/DDcjcs45l4AmOXN8+PDh5sNxnXOudiRNMrPhuer5jZyca2DjZi9nRsna3BVdgztmSBFfHdAl6TAaPU8czjWgbTsq+METU9i8rRxVN17QJcYMPvp0NU9ed0TSoTR6njica0BTF69l87Zy/nr5YYw6oLq7n7qkXP63j/h8e003X3SZfJFD5xrQe8UrSQmO2Ktr0qG4KiQor2h6fb754InDuQb0/vyVHNC7Ix33aJl0KK6KdEo0xcFC+eCJw7kGsmnrDqYsWsuRg7olHYqrRkqi3BNHLJ44nGsg4xeuZkeFcdTe3kzVGKUkyiuSjqIweOJwroF8MH8VrdIphvf34Z6NUTqFN1XF5InDuQbyXvFKDu3fibat0kmH4qoRXXF44ojDE4dzDWDVxq3MKl3PUd6/0WilUt7HEZcnDucawLvFKwE4dogv+d9YpSU8b8TjicO5BvDOvJV0bNuSA3p3TDoUV4OUz+OIzROHc3lmZrw7byVH792NdMrXGWmsUinv44jLE4dzeVa8YiPL1m/hmMHev9GYRU1Vnjji8MThXJ69PS/q3zjaE0ej5hMA4/PE4VyevTOvjL2K9qRP5z2SDsVlETVVJR1FYfDE4Vwebd1RzocLVnHsYB9N1dj5BMD4PHE4l0eTFq5hy/YKjt7bm6kaO2+qis8Th3N59Pa8lbRIicMH+fpUjZ3PHI/PE4dzefRucRmH9u9Mu9Z+z7TGLuUTAGPLW+KQ1EbSeEnTJM2S9ItQPlDSR5LmSXpSUqtQ3jo8Lw6vD8jY109D+VxJp+UrZufq06qNW5m5ZD3H+miqgpBO+QTAuPJ5xbEVONHMDgaGAaMkHQ78DrjbzAYDa4BrQv1rgDVmtjdwd6iHpP2Ai4H9gVHAXyT5KnGu0atcZuQY7xgvCKmUqPBLjljyljgssjE8bRl+DDgReCqUPwicEx6fHZ4TXj9JkkL5E2a21cw+BYqBEfmK27n68s68lXTaw5cZKRQpeeKIK699HJLSkqYCK4BxwHxgrZntCFVKgN7hcW9gMUB4fR3QNbO8mm0yj3WtpImSJpaVleXjdJyLzcx4Z14ZRw3yZUYKRdo7x2PLa+Iws3IzGwb0IbpKGFpdtfC7uv9dlqW86rHuNbPhZja8qMibBlyy5i7fwPL1Wzl2iPdvFIqoqSrpKApDg4yqMrO1wJvA4UAnSZVDTPoApeFxCdAXILzeEVidWV7NNs41Sm98HF31Hr9P94QjcXFVXhhWePbIKZ+jqookdQqP2wInA3OAN4DzQ7UrgefC47HhOeH1f1k0jXMscHEYdTUQGAyMz1fcztWHN+auYL+eHejRoU3SobiY0ooyh08CzC2fg8t7Ag+GEVApYIyZ/VPSbOAJSb8GpgD3h/r3Aw9LKia60rgYwMxmSRoDzAZ2ANebWXke43Zut6z7fDuTPlvDd47bK+lQXC2kwiWHd5DnlrfEYWbTgUOqKV9ANaOizGwLcEEN+7oNuK2+Y3QuH96dt5LyCuMEb6YqKKlwxVHhCx3m5DPHnatnb85dQYc2LRjWt1PSobhaSIdvQ2+qys0Th3P1qKLCePOTMo4dUkSLtP/3KiQ7rzg8ceTk/7Kdq0ezl66nbMNWb6YqQLuaqjxx5OKJw7l69MbHKwA4bh+fS1RoKidq+iTA3DxxOFeP3vykjIP7dKRbu9ZJh+JqadeoqoQDKQCeOJyrJ2s2bWPKojUc581UBWnnBEDv48jJE4dz9eTteWVUGJzgzVQFKe2d47F54nCunrw5t4wue7bioD4+DLcQVXaOex9Hbp44nKsHFRXGW5+UcdyQIl8Nt0Dt7OPwCYA5eeJwrh5MWbyW1Zu2cbw3UxWsymk33lSVmycO5+rBa3OW0yIljh/iHeOFKuWLHMbmicO5evDa7OWMGNiFjnu0TDoUV0c+ATA+TxzO7aaFKzcxb8VGTh7aI+lQ3G5I+zyO2DxxOLebXpuzHIBT9vPEUcgqxzT4qKrcPHE4t5tem7Ocfb/Snr5d9kg6FLcbfJHD+DxxOLcb1m7exoSFa7yZqglI+42cYvPE4dxueHNuGeUVxsneTFXwfAJgfJ44nNsN42Yvp6h9aw7q3THpUNxu8lvHxueJw7k62rqjnLc+KePkod13fum4wrVrraqEAykAnjicq6OPFqxm49Yd3r/RRPioqvjyljgk9ZX0hqQ5kmZJ+kEov1XSEklTw88ZGdv8VFKxpLmSTssoHxXKiiXdnK+YnauN1+Ysp03LFEft3S3pUFw92LVWlSeOXFrkcd87gB+Z2WRJ7YFJksaF1+42szsyK0vaD7gY2B/oBbwmaUh4+c/AKUAJMEHSWDObncfYncvKzHht9nKOGVxEm5bppMNx9cAnAMaXtysOM1tqZpPD4w3AHKB3lk3OBp4ws61m9ilQDIwIP8VmtsDMtgFPhLrOJWZW6XpK123h5KG+NlVTsbOpyjvHc2qQPg5JA4BDgI9C0Q2SpksaLalzKOsNLM7YrCSU1VRe9RjXSpooaWJZWVk9n4FzX/TyzGWkhPdvNCE+ATC+vCcOSe2Ap4EbzWw9cA8wCBgGLAXurKxazeaWpfyLBWb3mtlwMxteVORLW7v8emnmUkYO7EpXv7d4k+GLHMaX18QhqSVR0njUzJ4BMLPlZlZuZhXAfURNURBdSfTN2LwPUJql3LlEzFu+gfllmzj9wK8kHYqrR5V9HD6qKrd8jqoScD8wx8zuyijvmVHt68DM8HgscLGk1pIGAoOB8cAEYLCkgZJaEXWgj81X3M7l8tLMZQCctr8njqYk5fM4YsvnqKqjgG8AMyRNDWU/Ay6RNIyouWkhcB2Amc2SNAaYTTQi63ozKweQdAPwCpAGRpvZrDzG7VxWL89cxmH9O9OjQ5ukQ3H1KOV3AIwtb4nDzN6l+v6JF7NscxtwWzXlL2bbzrmGsmjVZmYvXc9/njE06VBcPUv7WlWx+cxx52rhpZlLARh1gDdTNTW+VlV8njicq4WXZi7jgN4d/N4bTZAPx43PE4dzMS1d9zlTF6/l9AN65q7sCs6upqqEAykAnjici+mVMJrKm6maJu8cj88Th3MxvTRzGUN6tGNQUbukQ3F54BMA4/PE4VwMKzduZcLC1YzyZqoma+cEQL/iyMkTh3MxvDxzGRUGp3szVZPlEwDj88ThXAzPTytl7+7t2Pcr7ZMOxeVJ5eq43lSVmycO53JYtm4L4xeu5qyDeiH5LWKbqrTP44gta+KQlJJ0ZEMF41xj9MKMpZjBmQd7/0ZTJp85HlvWxBFWsL0zWx3nmrqx00rZv1cHH03VxPkVR3xxmqpelXSe/BrdNUOLVm1m2uK1/NvBvZIOxeVZ2jvHY4uzyOEPgT2BckmfEy1caGbWIa+ROdcIPD89uvXL1w7yZqqmrvJPY2+qyi1n4jAzH0bimq3np5VyWP/O9Onsa1M1dTubqjxx5JSzqUqRyyX9PDzvK2lEru2cK3Tzlm/g42UbOMuvNpoFb6qKL04fx1+AI4BLw/ONwJ/zFpFzjcTz05eSEpzhiaNZ2NlU5Z3jOcXp4xhpZodKmgJgZmvCLVyda7LMjH9OK+XwvbrSvb3f6a85kERK3lQVR5wrju2S0kS3ekVSEeALD7smbVbpehas3MRZPpqqWUmn5MNxY4iTOP4APAt0l3Qb8C7w33mNyrmE/WPKElqmxaj9fW2q5kSSN1XFEGdU1aOSJgEnEQ3FPcfM5uQ9MucSsqO8gn9MLeWEfbrTeU9vlW1O0pI3VcWQM3FI+iXwDvCAmW3Kf0jOJevd4pWs3LiVcw/tk3QoroFFTVVJR9H4xWmqWghcAkyUNF7SnZLOzrVRGLb7hqQ5kmZJ+kEo7yJpnKR54XfnUC5Jf5BULGm6pEMz9nVlqD9P0pV1O1Xn4nlm8hI67dGSE/YtSjoU18AknwAYR87EYWajzeybwAnAI8AF4XcuO4AfmdlQ4HDgekn7ATcDr5vZYOD18BzgdGBw+LkWuAeiRAPcAowERgC3VCYb5+rbhi3beWXWMs48qCetW6STDsc1MO8cjyfOBMC/SXqf6Iu8BXA+kPOL28yWmtnk8HgDMAfoDZwNPBiqPQicEx6fDTxkkQ+BTpJ6AqcB48xstZmtAcYBo2pxjs7F9tKMZWzdUeHNVM1UWp444ojTVNUVSANrgdXASjPbUZuDSBoAHAJ8BPQws6UQJRege6jWG1icsVlJKKupvOoxrpU0UdLEsrKy2oTn3E7PTClhYLc9OaRvp6RDcQmQRLlPNsgpTlPV181sJHA70Al4Q1JJ3ANIagc8DdxoZuuzVa3u8FnKq8Z5r5kNN7PhRUXeNu1qr2TNZj5csJpzD+ntN2xqptKpaPKnyy7OqKozgWOAY4maqP5FNMoqJ0ktiZLGo2b2TCheLqmnmS0NTVErQnkJ0Ddj8z5AaSg/vkr5m3GO71xt/GPKEgDOOeRLF7SumUhJ3jkeQ5ymqtOBycB5ZravmV1tZqNzbRTu33E/MMfM7sp4aSxQOTLqSuC5jPIrwuiqw4F1oSnrFeBUSZ1Dp/ipocy5emNmPDN5CSMGdqFvF18Jt7lK+QTAWOJMALxeUg/gq2GI7HgzW5FrO+Ao4BvADElTQ9nPgN8CYyRdAywiGqUF8CJwBlAMbAauDsdfLelXwIRQ75dmtjrW2TkX09TFa1mwchPXHbdX0qG4BKVTwvNGbnGaqi4A7iBqHhLwR0k3mdlT2bYzs3epvn8ColnoVesbcH0N+xoN5LzKca6unppUQusWKU4/0FfCbc5SPo8jljir4/4X8NXKq4ywyOFrQNbE4Vyh2LxtB2OnlvK1A3vSoU3LpMNxCUqlvKkqjjh9HKkqTVOrYm7nXEF4ccYyNmzdwUVf7Zu7smvS0pKPqoohzhXHy5JeAR4Pzy8i6o9wrkl4csIiBnbbkxEDuyQdikuYj6qKJ07n+E2SzgWOJuqzuNfMns17ZM41gOIVG5mwcA03n76vz91wUVOVTwDMKc4VB8D7QDnRDZwm5KjrXMEYM3ExLVLiPF9ixOETAOOKs1bVt4DxwNeJ1qn6UNI38x2Yc/m2bUcFT08q4eShPShq3zrpcFwj4PM44olzxXETcIiZrQKQ1JXoCsSHx7qC9vqc5azatI2LRninuIt4H0c8cUZHlQAbMp5v4IuLDjpXkJ6YsJieHdtw7GBf28xFfAJgPHGuOJYAH0l6jmhxwbOB8ZJ+CFBlORHnCsKStZ/z9rwyvnfiYNIp7xR3EZ8AGE+cxDE//FSqXFuqff2H41zDGDMhumi+4DDvFHe7pPx+HLHEGY77i4YIxLmGsr28gsfHL+K4IUW+oKH7gpTEjgofj5uLzwB3zc6rs5azYsNWrjiif9KhuEYmnfLO8Tg8cbhm56EPFtK3S1uOG9I9Z13XvKRSwvNGbnHmcRwVp8y5QjB32QY++nQ1l43s753i7ktSwvs4YohzxfHHmGXONXoPf7iQVi1SXDjc5264L0v7PI5Yauwcl3QEcCRQVDn0NugApPMdmHP1bcOW7Tw7eQlnHdSLLnu2Sjoc1wh5U1U82UZVtQLahTqZQ2/XEy094lxBeWbyEjZtK/dOcVejlKDCM0dONSYOM3sLeEvSA2b2WQPG5Fy9MzMe/vAzDu7TkYP7dko6HNdIpf1GTrHEmQD4gKQvvZNmdmIe4nEuLz5YsIriFRu544KDkw7FNWI+ATCeOInjxxmP2wDnATvyE45z+fHg+wvpvEdLzjzI7ynuapaSvKkqhpyjqsxsUsbPe2b2Q2Bkru0kjZa0QtLMjLJbJS2RNDX8nJHx2k8lFUuaK+m0jPJRoaxY0s11OEfXzC1cuYlXZy/nspH9adPSx3W4mnlTVTw5rzgkZd5PMwUcBnwlxr4fAP4EPFSl/G4zu6PKMfYDLgb2B3oBr0kaEl7+M3AK0Sq9EySNNbPZMY7vHAB/f+9TWqZSXHGkd4q77KIrjqSjaPziNFVNIloVV0RNVJ8C1+TayMzeljQgZhxnA0+Y2VbgU0nFwIjwWrGZLQCQ9ESo64nDxbJ28zbGTCzh34b1onv7NkmH4xo5nwAYT5xFDgfW8zFvkHQFMBH4kZmtAXoDH2bUKQll8MV7f5RQQzOZpGuBawH69etXzyG7QvXY+EV8vr2ca46u73/GrinytariibPkSBtJP5T0jKSnJf27pLr+6XYPMAgYBiwF7qw8TDV1LUv5lwvN7jWz4WY2vKjIb8zjolvDPvj+Qo4Z3I2hPTskHY4rAD4BMJ44S448RNT38EeiPouhwMN1OZiZLTezcjOrAO5jV3NUCZC5BkQfoDRLuXM5/XN6KcvXb/WrDRebN1XFE6ePYx8zyxz8/oakaXU5mKSeZrY0PP06UDniaizwmKS7iDrHBwPjia44BksaSHQnwouBS+tybNe8mBn3vfMpg7u347ghfgXq4kn7PI5Y4iSOKZION7MPASSNBN7LtZGkx4HjgW6SSoBbgOMlDSNqbloIXAdgZrMkjSHq9N4BXG9m5WE/NwCvEK2PNdrMZtXqDF2z9MH8VcxZup7fnXcgkq+C6+KRL3IYS5zEMRK4QtKi8LwfMEfSDMDM7KDqNjKzS6opvr+mg5jZbcBt1ZS/CLwYI07ndvq/by+g656tOHtY79yVnQvSKZ8AGEecxDEq71E4V49mLlnHW5+UcdNp+/iEP1crae8cjyVO4vi1mX0js0DSw1XLnGss/vxGMe3btOAbvgquqyUJnzkeQ5xRVftnPpHUgmj2uHONTvGKDbw8axlXHjGADm1aJh2OKzBpX6sqlhoTR1g7agNwkKT1kjaE58uB5xosQudq4S9vzKdNizTf9CG4rg6ipipPHLnUmDjM7Ddm1h74vZl1MLP24aermf20AWN0LpZFqzbz3LRSLhvZz+/w5+pEivo4zJNHVnH6OF6SdGzVQjN7Ow/xOFdnf317PmmJbx+7V9KhuAKVDkO3KwzSPoq7RnESx00Zj9sQzfaeBPiNnFyjsWzdFp6aWMIFw/vQo4MvZujqJh3aYCrMSFe74pGDeIscnpX5XFJf4Pa8ReRcHdz3zgLKzfjOcYOSDsUVsMrJouUVho/krlmcUVVVlQAH1HcgztXVivVbeOTDzzhnWG/6dtkj6XBcAUunKpuqvI8jmzg3cvoju1akTRGtbFuntaqcy4c/vVFMeYXxg5MGJx2KK3CZfRyuZnH6OCZmPN4BPG5mOdeqcq4hlKzZzOPjF3HhV/vSr6tfbbjdU7msma9XlV2cxPEksDfRVcd8M9uS35Cci+8Pr89DEt87ce+kQ3FNwM6mKk8cWWWbANhC0u1EfRoPAo8AiyXdLsmn5LrELSjbyNOTl3D5yP707Ng26XBcE+B9HPFk6xz/PdAFGGhmh5nZIUR37+sE3NEQwTmXzf+8No9W6RTfPd5HUrn6sXNUlSeOrLIljjOBb5vZhsoCM1sPfBc4I9+BOZfNx8vW8/z0Uq46agBF7VsnHY5rIio7xz1vZJctcZhVM+8+3GDJ31aXqLte/YR2rVpwnc8Sd/Uo5Z3jsWRLHLMlXVG1UNLlwMf5C8m57MZ/uppXZy/n28fuRac9fE0qV39SqV0TAF3Nso2quh54RtI3iZYYMeCrQFui+4U71+AqKozbXpjNVzq04dvH+NWGq1/eVBVPjYnDzJYAIyWdSHRPDgEvmdnrDRWcc1U9P72UaSXruOOCg2nbyteEcPWrYmjnAAASkklEQVQrFdpgvHM8uzhrVf0L+FcDxOJcVlu2l3P7y3PZv1cHzj3E7yXu6l9K3lQVR13WqopF0mhJKyTNzCjrImmcpHnhd+dQLkl/kFQsabqkQzO2uTLUnyfpynzF6xq/0e99ypK1n/OfXxu6sy3aufpUOY/D78eRXd4SB/AAMKpK2c3A62Y2GHg9PAc4HRgcfq4F7oEo0QC3ACOJlnO/pTLZuOZl5cat/OWN+Zw8tDtHDuqWdDiuiUr5PI5Y8pY4wo2eVlcpPptoFjrh9zkZ5Q9Z5EOgk6SewGnAODNbbWZrgHF8ORm5ZuB/XvuEz7eXc/PpQ5MOxTVh3lQVTz6vOKrTw8yWAoTf3UN5b2BxRr2SUFZT+ZdIulbSREkTy8rK6j1wl5zZpet57KNFXDayH3t3b5d0OK4J29VUlXAgjVxDJ46aVNdgbVnKv1xodq+ZDTez4UVFRfUanEtORYXx8+dm0nmPVvzolH2SDsc1cT4BMJ6GThzLQxMU4feKUF4C9M2o1wcozVLumomnJpcw6bM13Hz6vnTcw9fWdPm1cwKgX3Jk1dCJYyxQOTLqSuC5jPIrwuiqw4F1oSnrFeBUSZ1Dp/ipocw1A2s3b+O3L33M8P6dOe/QPkmH45qBXRMAPXFkE+d+HHUi6XHgeKCbpBKi0VG/BcZIugZYBFwQqr9ItHBiMbAZuBrAzFZL+hUwIdT7pZlV7XB3TdTvX5nLus+386tzDvDht65B7OocTziQRi5vicPMLqnhpZOqqWtES5xUt5/RwOh6DM0VgGmL1/LY+EVcfeRAhvbskHQ4rpmonDnu9+PIrrF0jju3U3noEC9q15p/P8XvI+4aTuUVh98BMDtPHK7R+ft7nzK9ZB3/deZ+tG/jHeKu4aS9czwWTxyuUfl05SZ+/8pcTh7ag7MO6pl0OK6Z2XnF4XkjK08crtGoqDB+8tR0WrdIcdvXD9h5G0/nGkrlGAxvqsrOE4drNB7+8DPGL1zNz8/cjx4d2iQdjmuG0n4jp1g8cbhGYfHqzfzu5Y85dkgR5x/mczZcMnY1VXniyMYTh0ucmXHzM9NJSfzm3AO9icolxhNHPJ44XOIeeH8h7xWv4qdn7EvvTm2TDsc1Y7uaqhIOpJHzxOESNWfpen7z4sectG93Lh3RL+lwXDOX9gmAsXjicInZsr2c7z8+hY57tOT28w/yJiqXOHlTVSx5W3LEuVxue2EO81Zs5KFvjqBru9ZJh+PczkUOfVRVdn7F4RIxbvZyHv7wM759zECOHeL3T3GNQ2Ufh+eN7DxxuAa3bN0W/uOpaezXswM/Ps1vzuQaD/kEwFg8cbgGtXVHOd99dBLbdlTwh0sOoXWLdNIhObeTr1UVj/dxuAb1q3/OZsqitfzlskP9/uGu0Ul753gsfsXhGsxTk0p45MNFXHfsXpxxoC9g6Bof+bLqsXjicA1i5pJ1/OezMzhir67c5P0arpHyzvF4PHG4vFuzaRvfeWQSXfdsxZ8uPYQWaf9n5xqnytVxfThudt7H4fJqy/Zyrn14Iis2bGXMdUf4fA3XqKVS3scRh//p5/KmosK46anpTFi4hrsuPJhhfTslHZJzWXnneDyJJA5JCyXNkDRV0sRQ1kXSOEnzwu/OoVyS/iCpWNJ0SYcmEbOrvTvHzeX5aaX8x6h9OPOgXkmH41xOKfkih3EkecVxgpkNM7Ph4fnNwOtmNhh4PTwHOB0YHH6uBe5p8EhdrY2ZsJg/vzGfS0b05bvHDUo6HOdiSfkih7E0pqaqs4EHw+MHgXMyyh+yyIdAJ0k+lrMRe3PuCn727AyOGdyNX57tt4B1hSPtw3FjSSpxGPCqpEmSrg1lPcxsKUD43T2U9wYWZ2xbEspcI/TB/FVc9/AkhvRoz18uO5SWPoLKFZCdTVV+xZFVUqOqjjKzUkndgXGSPs5St7o/V7/0qYYEdC1Av35+X4ckTF60hmsenEC/Lnvw8DUjaN+mZdIhOVcrO0dV+RVHVon8OWhmpeH3CuBZYASwvLIJKvxeEaqXAH0zNu8DlFazz3vNbLiZDS8q8tVWG9qs0nVcNXo8Re1b88i3RvqwW1ew0in5BMAcGjxxSNpTUvvKx8CpwExgLHBlqHYl8Fx4PBa4IoyuOhxYV9mk5RqHT5Zv4Bv3j6dd6xY8+q2R9OjQJumQnKuzlLypKpckmqp6AM+GDtMWwGNm9rKkCcAYSdcAi4ALQv0XgTOAYmAzcHXDh+xqMqNkHVeM/ogW6RSPfvtw+nTeI+mQnNstKcmbqnJo8MRhZguAg6spXwWcVE25Adc3QGiulsZ/uppvPjCBjm1b8ui3RjKg255Jh+TcbouaqjxxZONLjrg6eXPuCr7zyCR6d2rLI98aSc+ObZMOybl6kZJ8AmAOnjhcrY2dVsqPxkxlSI/2fr9w1+Sk5BMAc/HE4WIzM/7wejF3v/YJIwZ04b4rh9OxrQ+5dU1LypuqcvLE4WLZsr2cnzw9neemlnLuob35zbkH+m1fXZOUlnxZ9Rw8cbicyjZs5bqHJzJ50Vr+Y9Q+fPe4Qb6MiGuyUj6PIydPHC6rjxas4vtPTGHd59v56+WHMuoAXybMNW0p+czxXDxxuGpVVBj3vDWfO1+dS/+ue/L3q0awX68OSYflXN6lJZ8AmIMnDvclKzdu5d+fnMo781Zy1sG9+M25B9Kutf9Tcc2Dd47n5t8G7gtenLGU//PcTNZv2cF/f/1ALhnR1/szXLPiM8dz88ThgOgq45bnZvHCjKUc0LsDj3zrYPb9ijdNueYnnRLlnjey8sTRzJkZY6eV8ovnZ7Nxyw5uOm0frj12L7+Phmu2fAJgbp44mrEZJev4xfOzmPjZGg7u24nfn38QQ3q0Tzos5xLlTVW5eeJohso2bOWOV+YyZtJiuu7Zit+ddyDnH9aXdMr7MpxLp3wCYC6eOJqR1Zu2cd87C3jw/YVsL6/g28fsxQ0n7k0Hv1Ofczul5BMAc/HE0QysyUgYm7eXc9ZBvbjx5MHsVdQu6dCca3RSKe/jyMUTRxO3o7yCs/70LkvWfs7XDuzJD04azGDvx3CuRr5WVW6eOJq4maXrKVnzObeffxAXDu+bewPnmjmfAJibj7ls4t6fvxKAE/bpnnAkzhWGqI/DE0c2njiauA/mr2JIj3YUtfebLTkXR1qiwu8AmJUnjiZs244KJixczZGDuiUdinMFQ8IXOczBE0cTNnXxWrZsr+CIQV2TDsW5gpFO+QTAXAomcUgaJWmupGJJNycdTyF4f/5KJDh8oCcO5+JKe+d4TgUxqkpSGvgzcApQAkyQNNbMZicbWf0wM8yiseNG+G1EPxgVFuoAVvHFsopQp7L+zn1UGO/MW8kBvTrScQ+f4OdcXJIvcphLQSQOYARQbGYLACQ9AZwN1GviWLt5G+f/9YOdX+SZX+KVv6laRjVf/BWV5V/8cq8woJoEkU/XHbdXfg/gXBOTFsxZup5T7nor6VDqZN+eHfjjJYfk9RiFkjh6A4sznpcAIzMrSLoWuBagX79+dTpIOiWG9GiHJBTtk5RAREP0UPS78rkU1ZEI9UL9cP+KVOZrlfV31tm1r137qHz+5bIvHI9d+0yFDSqPnxlLOiVOHtqjTu+Fc83VpSP707ZVOukw6qxv57Z5P0ahJI7qVt/7wt/qZnYvcC/A8OHD6/R3fPs2LfnLZYfVZVPnXBNxyn49OGU//4Mrm0LpHC8BMqc99wFKE4rFOeeatUJJHBOAwZIGSmoFXAyMTTgm55xrlgqiqcrMdki6AXgFSAOjzWxWwmE551yzVBCJA8DMXgReTDoO55xr7gqlqco551wj4YnDOedcrXjicM45VyueOJxzztWKrAku5iWpDPhsN3bRDVhZT+Ekzc+lcfJzaZya+7n0N7OiXJWaZOLYXZImmtnwpOOoD34ujZOfS+Pk5xKPN1U555yrFU8czjnnasUTR/XuTTqAeuTn0jj5uTROfi4xeB+Hc865WvErDuecc7XiicM551yteOLIIGmUpLmSiiXdnHQ8u0PSQkkzJE2VNDHpeGpL0mhJKyTNzCjrImmcpHnhd+ckY4yrhnO5VdKS8PlMlXRGkjHGIamvpDckzZE0S9IPQnnBfS5ZzqXgPhcASW0kjZc0LZzPL0L5QEkfhc/myXBbit0/nvdxRCSlgU+AU4huHDUBuMTM6vW+5g1F0kJguJkV5GQmSccCG4GHzOyAUHY7sNrMfhsSe2cz+0mSccZRw7ncCmw0szuSjK02JPUEeprZZEntgUnAOcBVFNjnkuVcLqTAPhcARfer3tPMNkpqCbwL/AD4IfCMmT0h6a/ANDO7Z3eP51ccu4wAis1sgZltA54Azk44pmbLzN4GVlcpPht4MDx+kOg/eqNXw7kUHDNbamaTw+MNwBygNwX4uWQ5l4JkkY3hacvwY8CJwFOhvN4+G08cu/QGFmc8L6GA/yER/aN5VdIkSdcmHUw96WFmSyH6jw90Tzie3XWDpOmhKavRN+9kkjQAOAT4iAL/XKqcCxTo5yIpLWkqsAIYB8wH1prZjlCl3r7TPHHsomrKCrkd7ygzOxQ4Hbg+NJe4xuMeYBAwDFgK3JlsOPFJagc8DdxoZuuTjmd3VHMuBfu5mFm5mQ0D+hC1oAytrlp9HMsTxy4lQN+M532A0oRi2W1mVhp+rwCeJfqHVOiWh7bpyjbqFQnHU2dmtjz8R68A7qNAPp/Qfv408KiZPROKC/Jzqe5cCvVzyWRma4E3gcOBTpIq7/Rab99pnjh2mQAMDqMQWgEXA2MTjqlOJO0ZOvyQtCdwKjAz+1YFYSxwZXh8JfBcgrHslsov2uDrFMDnEzpg7wfmmNldGS8V3OdS07kU4ucCIKlIUqfwuC1wMlG/zRvA+aFavX02PqoqQxh69z9AGhhtZrclHFKdSNqL6CoDovvKP1Zo5yLpceB4oqWhlwO3AP8AxgD9gEXABWbW6DudaziX44maQwxYCFxX2U/QWEk6GngHmAFUhOKfEfUNFNTnkuVcLqHAPhcASQcRdX6niS4IxpjZL8N3wRNAF2AKcLmZbd3t43nicM45VxveVOWcc65WPHE455yrFU8czjnnasUTh3POuVrxxOGcc65WPHG4JkfS3ZJuzHj+iqS/ZTy/U9IPJfWS9FT1e6lx31dJ+lN9xltbkgZkrrRbx330lPRqfcXkmhdPHK4peh84EkBSimj+xP4Zrx8JvGdmpWZ2fjXbNwejgFeSDsIVJk8cril6j5A4iBLGTGCDpM6SWhOt4TMl8y/3cCXxjKSXw70Lbq/cmaSrJX0i6S3gqOoOKOm4jHs4TJHUXtLxkt6W9Kyk2ZL+GhIZkk6V9IGkyZL+X1gzCUmHSXorLE75SsZSHoeFey18AFxfQwzHh23HhHh/K+kyRfdpmCFpUEb1UcBL4crj7RD3TEnH1P1td82FJw7X5IR1unZI6keUQD4gmt18BDAcmB6Wzq9qGHARcCBwkaKb/fQEfkGUME4B9qvhsD8Grg+LzB0DfB7KRwA/CvscBJwrqRvwX8DJYSHKicAPw9pJfwTON7PDgNFA5Yz/vwPfN7Mjcpz+wUT3YTgQ+AYwxMxGAH8Dvgc77z2zT7jXzKXAKyHug4GpOfbvHC1yV3GuIFVedRwJ3EW0nPSRwDqipqzqvG5m6wAkzQb6EzVzvWlmZaH8SWBIDce7S9KjRDfOKYmWQ2K8mS0I2z4OHA1sIUpA74U6rYiS2z7AAcC4UJ4GlkrqCHQys7fCsR4mWvW4OhMql8iQNB+o7MeYAZwQHo9k1xLiE4DRIWn9w8w8cbic/IrDNVWV/RwHEjVVfUh0xXEk0Zd8dTLX8Cln1x9WOdflMbPfAt8C2gIfStq3hm2NaAn/cWY2LPzsZ2bXhPJZGeUHmtmpoTzu2kCZ51CR8bwi43xOB14Ocb8NHAssAR6WdEXM47hmzBOHa6reA84kuqVpeVh0rxNR8vigFvv5CDheUtfwV/kF1VWSNMjMZpjZ74ianioTx4iw4nKKqBnsXaIkdpSkvcO2e0gaAswFiiQdEcpbSto/LJO9LizMB3BZLeKvzknA6+EY/YEVZnYf0Wqxh+7mvl0z4InDNVUziJqZPqxStq4292EPzT63EiWb14DJNVS9MXQuTyPq33gplH8A/JboqudT4NnQ7HUV8Lik6SHGfUO/y/nA78J+prKrk/9q4M+hc7yy/6TWJBUBWzJuwHQ8MFXSFOA84H/rum/XfPjquM7liaTjgR+b2ZlJx1JJ0uVAn9C05lydeOe4c82ImT2SdAyu8PkVh3POuVrxPg7nnHO14onDOedcrXjicM45VyueOJxzztWKJw7nnHO18v8Bxh4AVZF3I+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12b31fa90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def production_curve(wind_speed):\n",
    "    \"\"\" Production curve of a turbine/farm can be simplified as a \n",
    "    cube of wind speed, with production \n",
    "         starting from a particular wind speed        'v_cut_in'\n",
    "         getting maximum output power at wind speed   'v_rated'\n",
    "         production halted at very strong wind speeds 'v_cut_off'\n",
    "    Note that this simplified approach is useful for comparing places only.\n",
    "    \"\"\"\n",
    "    v_cut_in = 5\n",
    "    v_rated = 15\n",
    "    v_cut_off = 25\n",
    "    max_production = 15**3\n",
    "    rt = np.zeros_like(wind_speed)\n",
    "    rt = np.where(wind_speed>v_cut_in, wind_speed**3,wind_speed)\n",
    "    rt = np.where(wind_speed>v_rated,max_production,rt)\n",
    "    rt = np.where(wind_speed>v_cut_off,0,rt)\n",
    "    return rt\n",
    "    \n",
    "def wind_production_smooth(wspd):\n",
    "    \"\"\"Weekly average production as a moving average\"\"\"\n",
    "    def moving_average(a, n=3) :\n",
    "        ret = np.cumsum(a, dtype=float)\n",
    "        ret[n:] = ret[n:] - ret[:-n]\n",
    "        return ret[n - 1:] / n\n",
    "    return moving_average(np.array(wspd), n=7*24)\n",
    "\n",
    "def station_statistics(lonlats, count = 100):\n",
    "    \"\"\"\n",
    "    Compute min, max, mean, 5'th and 95'th percentiles for both \n",
    "    hourly values and weekly averages\n",
    "    \"\"\"\n",
    "    ddwd = [(name,ds.get_json_data_in_pandas(count=count, \n",
    "                **{'vars':'100_metre_U_wind_component_surface,100_metre_V_wind_component_surface',\n",
    "                'lon':lon,\n",
    "                'lat':lat})) for name,lon,lat in lonlats]\n",
    "    wspds = [(name,np.sqrt(dd['100_metre_U_wind_component_surface']**2 + dd['100_metre_V_wind_component_surface']**2))\n",
    "              for name,dd in ddwd]\n",
    "    pcurves = [(name,production_curve(dd)) for name, dd in wspds]\n",
    "    weekly_ave = [wind_production_smooth(dd) for name,dd in pcurves]    \n",
    "    retdic = {}\n",
    "    for dd, wa in zip(pcurves,weekly_ave):\n",
    "        retdic[dd[0]] = (np.amin(dd[1]),np.amax(dd[1]),np.mean(dd[1]),np.percentile(dd[1],5),np.percentile(dd[1],95),\n",
    "                         np.amin(wa), np.amax(wa), np.mean(wa), np.percentile(wa,5), np.percentile(wa,95),\n",
    "                        np.sum(dd[1]))                    \n",
    "        \n",
    "    return pd.DataFrame(retdic, index=['min','max','mean','5 percentile', '95 percentile',\n",
    "                                       'wa min','wa max','wa mean','wa 5 percentile', 'wa 95 percentile', 'sum'])\n",
    "\n",
    "## Power curve demo\n",
    "x = np.arange(0,30, 0.2)\n",
    "fig = plt.figure()\n",
    "plt.plot(x,production_curve(x))\n",
    "plt.xlabel('Wind speed m/s')\n",
    "plt.ylabel('Output power')\n",
    "plt.title(\"Idealized turbine production curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stations = [('Hiiumaa',22.1, 59),            ## planned off-shore wind farm in North-West Estonia\n",
    "            ('Tõravere',26+28/60, 58+16/60), ## Inland climate station in South-East Estonia\n",
    "            ('GYM',-3 -35/60, 53+27/60)]     ## existing off-shore wind farm west to Scotland\n",
    "abc = station_statistics(stations,count=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>min</th>\n",
       "      <th>max</th>\n",
       "      <th>mean</th>\n",
       "      <th>5 percentile</th>\n",
       "      <th>95 percentile</th>\n",
       "      <th>wa min</th>\n",
       "      <th>wa max</th>\n",
       "      <th>wa mean</th>\n",
       "      <th>wa 5 percentile</th>\n",
       "      <th>wa 95 percentile</th>\n",
       "      <th>sum</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>GYM</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>3375.0</td>\n",
       "      <td>781.445151</td>\n",
       "      <td>2.359436</td>\n",
       "      <td>3375.000000</td>\n",
       "      <td>61.956120</td>\n",
       "      <td>2398.221162</td>\n",
       "      <td>780.378387</td>\n",
       "      <td>202.650405</td>\n",
       "      <td>1737.967417</td>\n",
       "      <td>7.814452e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Hiiumaa</th>\n",
       "      <td>0.050631</td>\n",
       "      <td>3375.0</td>\n",
       "      <td>770.282982</td>\n",
       "      <td>2.309145</td>\n",
       "      <td>2944.666319</td>\n",
       "      <td>71.048125</td>\n",
       "      <td>2062.972979</td>\n",
       "      <td>777.684667</td>\n",
       "      <td>231.585537</td>\n",
       "      <td>1564.945427</td>\n",
       "      <td>7.702830e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tõravere</th>\n",
       "      <td>0.117344</td>\n",
       "      <td>3375.0</td>\n",
       "      <td>403.100573</td>\n",
       "      <td>1.921589</td>\n",
       "      <td>1532.588197</td>\n",
       "      <td>22.684435</td>\n",
       "      <td>1380.830043</td>\n",
       "      <td>405.582212</td>\n",
       "      <td>65.306773</td>\n",
       "      <td>835.503148</td>\n",
       "      <td>4.031006e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               min     max        mean  5 percentile  95 percentile  \\\n",
       "GYM       0.000000  3375.0  781.445151      2.359436    3375.000000   \n",
       "Hiiumaa   0.050631  3375.0  770.282982      2.309145    2944.666319   \n",
       "Tõravere  0.117344  3375.0  403.100573      1.921589    1532.588197   \n",
       "\n",
       "             wa min       wa max     wa mean  wa 5 percentile  \\\n",
       "GYM       61.956120  2398.221162  780.378387       202.650405   \n",
       "Hiiumaa   71.048125  2062.972979  777.684667       231.585537   \n",
       "Tõravere  22.684435  1380.830043  405.582212        65.306773   \n",
       "\n",
       "          wa 95 percentile           sum  \n",
       "GYM            1737.967417  7.814452e+06  \n",
       "Hiiumaa        1564.945427  7.702830e+06  \n",
       "Tõravere        835.503148  4.031006e+06  "
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abc.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}