web_services_data_access.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "R-M5_FCS_ZS0"
   },
   "source": [
    "# Accessing Data through ORNL DAAC Web Services\n",
    "## Exploring Forest Disturbance Caused by Bark Beetles at Fairview Curve, Rocky Mountain National Park"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1vYq8f_eFsOb"
   },
   "source": [
    "*Author: ORNL DAAC*  \n",
    "*Date: April 21, 2019*  \n",
    "*Contact for ORNL DAAC: uso@daac.ornl.gov*  \n",
    "\n",
    "**Keywords:** Web service, Python, Open Geospatial Consortium (OGC), REST API, MODIS, Daymet, Spatial Data Access Tool (SDAT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UNAWghRR_koy"
   },
   "source": [
    "## Web Services at the ORNL DAAC \n",
    "The [ORNL DAAC](https://daac.ornl.gov/) (Oak Ridge National Laboratory Distributed Active Archive Center) for Biogeochemical Dynamics is a National Aeronautics and Space Administration (NASA) Earth Observing System Data and Information System (EOSDIS) data center managed by the Earth Science Data and Information System  (ESDIS) Project.\n",
    "\n",
    "The ORNL DAAC offers many [tools and services](https://daac.ornl.gov/tools/) for visualization and access to data, which will be demonstrated in this tutorial.\n",
    "\n",
    "The Web services used in this tutorial include:\n",
    "* [MODIS Web Service](https://modis.ornl.gov/data/modis_webservice.html), a [REST](https://en.wikipedia.org/wiki/Representational_state_transfer) Web service, to access [MODIS Land Products](https://modis.ornl.gov/documentation.html)\n",
    "* [Daymet Single Pixel Tool](https://daymet.ornl.gov/web_services.html), a [REST](https://en.wikipedia.org/wiki/Representational_state_transfer) web service to access daily weather parameters offered in the [Daymet data product](https://daymet.ornl.gov/) at a single geographic point\n",
    "* [Spatial Data Access Tools (SDAT)](https://webmap.ornl.gov/ogc) offers [OGC](http://www.opengeospatial.org/)-based [REST](https://en.wikipedia.org/wiki/Representational_state_transfer) Web services to visualize and download spatial data in various user-selected spatial/temporal extents and formats. The SDAT is a general-purpose Web application. It contains 187 archived data products (as of April 2019) in many science disciplines, including agriculture, biosphere, climate, land surface, etc.\n",
    "\t* [Web Map Service (WMS)](https://en.wikipedia.org/wiki/Web_Map_Service): Get maps (i.e., visualization) of geospatial data from a remote server\n",
    "\t* [Web Coverage Service (WCS)](https://en.wikipedia.org/wiki/Web_Coverage_Service): Download geospatial data from a remote server\n",
    "\n",
    "With both OGC WMS and WCS, you can customize maps for the data you need. For example, you can choose the spatial region, the data and map format (e.g., GeoTIFF for WCS; PNG for WMS), the coordinate reference system or projection, and the resolution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jqJM86mfFjq8"
   },
   "source": [
    "## Overview\n",
    "In this tutorial, we will use several Web services to access and combine datasets near Fairview Curve at Rocky Mountain National Park in Colorado (latitude = 40.399391, longitude = -105.835901) to investigate tree mortality owing to bark beetle infestation from 2006-2010.\n",
    "\n",
    "![Image_Bark_Beetle_Farview_Curve](https://www.nps.gov/romo/learn/nature/images/ForestHealth_MPB_Comparison_688x300.jpg)\n",
    "\n",
    "\n",
    "The photos above were taken in late summer in 2005 (left) and 2009 (right), and illustrate the effects of bark beetle infestation. One way to quantitatively access the impacts of bark beetles on tree mortality is to use a vegetation index, which indicated the health of vegetation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A4BlYkYEGBWs"
   },
   "source": [
    "### Data Products\n",
    "Data products used in theis tutorial include:\n",
    "* [MOD13Q1](https://lpdaac.usgs.gov/products/mod13q1v006/): MODIS (Collection 6) 16 day Enhanced Vegetation Index (EVI) at 250 m spatial resolution.\n",
    "\t* Accessed through MODIS Web service\n",
    "* [Daymet](https://daymet.ornl.gov/overview): Daily total precipitation and daily maximum/minimum temperature at 1 km spatial resolution.\n",
    "\t* Accessed through Daymet single pixel Web service\n",
    "* [Tree Mortality from Bark Beetles](https://doi.org/10.3334/ORNLDAAC/1512): yearly total tree motality (in Mg C km2) at 1 km Resolution in the western USA for 2003-2012\n",
    "\t* Accessed through SDAT OGC WMS and WCS\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PvzrUCiOGIK9"
   },
   "source": [
    "### Prerequisites\n",
    "Python 2.7 or later and modules: numpy, pandas, matplotlib, urllib, zeep, owslib, datetime, modisViirsClient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bz4LraMfG-iu"
   },
   "source": [
    "## Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IgmV7y7G6pJO"
   },
   "source": [
    "We will use [modisViirsClient](https://github.com/tquaife/modisViirsClient)  to call ORNL DAAC's MODIS REST service. The client has been cloned into this repository for use in the tutorial. \n",
    "\n",
    "Let's import some python modules to be used in this Jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-0GyPbzTApM6"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done!\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "# above generates plots in line within this page\n",
    "\n",
    "import numpy as np # numpy module for array operations\n",
    "import pandas as pd # pandas module for tabular data\n",
    "import matplotlib.pyplot as plt # matplotlib module for plotting\n",
    "import datetime as dt # datetime module for converting dates\n",
    "import json  # json lets you produce REST output in readable format\n",
    "import urllib # urllib allows you to access legend images from SDAT\n",
    "from owslib.wms import WebMapService # OWSlib module to access WMS services from SDAT\n",
    "from owslib.wcs import WebCoverageService  # OWSlib module to access WCS services from SDAT\n",
    "from IPython.display import Image # Displaying image\n",
    "\n",
    "# import modisViirsClient\n",
    "import sys\n",
    "sys.path.insert(0, 'modisViirsClient')\n",
    "from modisViirsClient import *\n",
    "\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lTIy3LeQ6y4_"
   },
   "source": [
    "Then define the location of the study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "k5m8qIc264Ek"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Latitude: 40.399391\n",
      "Longitude: -105.835901\n"
     ]
    }
   ],
   "source": [
    "# Fairview Curve, Rocky Mountain National Park, Colorado\n",
    "latitude = 40.399391\n",
    "longitude = -105.835901\n",
    "\n",
    "# Deep Creek Ridge in Central Idaho \n",
    "# latitude = 45.0949 \n",
    "# longitude = -114.1086\n",
    "\n",
    "print('Latitude: ' + str(latitude))\n",
    "print('Longitude: ' + str(longitude))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VP3SUiLV8Ugj"
   },
   "source": [
    "### 1. Download Weather Data using the Daymet Single Pixel Tool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Np6Ua_ycCWoI"
   },
   "source": [
    "The Daymet Single Pixel Tool requires a url to be formated like: \n",
    "> https://<i></i>daymet.ornl.gov/data/send/saveData?lat=**lat**&lon=**lon**&measuredParams=**params**&year=**years**\n",
    "\n",
    "where **lat** and **lon** are the latitude and longitude values in decimal degrees, **params** is a comma-separated list of parameters (tmax, tmin, dayl, prcp, srad, swe, vp) and **years** is a comma-separated list of years.\n",
    "\n",
    "Here, we will download daily precipitation data (**prcp**), daily minimum temperature data (**tmin**), and daily maximum temperature data (**tmax**) from 2003-2011 at our location of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "pJxKZFRoBizJ"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Firt 5 rows of daymetdf:\n",
      "              year  yday  prcp (mm/day)  tmax (deg c)  tmin (deg c)\n",
      "2003-01-01  2003.0   1.0            6.0          -5.0         -15.0\n",
      "2003-01-02  2003.0   2.0            3.0          -4.5         -16.0\n",
      "2003-01-03  2003.0   3.0            3.0          -0.5          -7.0\n",
      "2003-01-04  2003.0   4.0            3.0          -0.5          -5.5\n",
      "2003-01-05  2003.0   5.0            0.0          -1.5          -7.0\n",
      "\n",
      "Firt 5 rows of daymetdf_30yrs:\n",
      "              year  yday  prcp (mm/day)  tmax (deg c)  tmin (deg c)\n",
      "1981-01-01  1981.0   1.0            0.0           5.0         -13.5\n",
      "1981-01-02  1981.0   2.0            0.0           6.0         -10.5\n",
      "1981-01-03  1981.0   3.0            0.0           3.5         -11.0\n",
      "1981-01-04  1981.0   4.0            0.0           3.5         -10.5\n",
      "1981-01-05  1981.0   5.0            0.0           2.5         -12.0\n"
     ]
    }
   ],
   "source": [
    "# list of years to retrieve data\n",
    "yrs = \"2003,2004,2005,2006,2007,2008,2009,2010,2011\"\n",
    "\n",
    "# Daymet URL. It is important that the URL is structured correctly\n",
    "daymeturl = \"https://daymet.ornl.gov/data/send/saveData?lat={lat}&lon={lon}&measuredParams=prcp,tmax,tmin&year={year}\"\n",
    "\n",
    "# read daymet_file, first 6 lines are headers\n",
    "daymetdf = pd.read_csv(daymeturl.format(lat=latitude, lon=longitude, year=yrs), header=6)\n",
    "\n",
    "# convert year and day of year into python datetime\n",
    "daymetdf.index = pd.to_datetime(daymetdf.year.astype(int).astype(str) + '-' + daymetdf.yday.astype(int).astype(str), format=\"%Y-%j\")\n",
    "\n",
    "# preview the first five rows\n",
    "print('Firt 5 rows of daymetdf:')\n",
    "print(daymetdf[:5])\n",
    "\n",
    "# retrieve Daymet data in 1981-2010 to calculate 30-yr normal\n",
    "yrs = \"1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010\"\n",
    "daymetdf_30yrs = pd.read_csv(daymeturl.format(lat=latitude, lon=longitude, year=yrs), header=6)\n",
    "daymetdf_30yrs.index = pd.to_datetime(daymetdf_30yrs.year.astype(int).astype(str) + '-' + daymetdf_30yrs.yday.astype(int).astype(str), format=\"%Y-%j\")\n",
    "\n",
    "# preview the first five rows\n",
    "print('\\nFirt 5 rows of daymetdf_30yrs:')\n",
    "print(daymetdf_30yrs[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MAS9QXu8Ee3r"
   },
   "source": [
    "Now, rename the columns into something shorter, create annual precipitation, and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1vNbJZCIChiA"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plots:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1YAAAEICAYAAABcR2XFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XmczfX+wPHXe+ykZEkiBu0iN0rdEsUt7ctNFBW3aC8tt1S3mxbS9iuUiohka1HpRosb3VQSly7tyZpdSGR///54f485xjkzZ8w55zsz5/18PM5j5rt+3ufMme/n+/l+NlFVnHPOOeecc87tvaywA3DOOeecc8654s4LVs4555xzzjlXSF6wcs4555xzzrlC8oKVc84555xzzhWSF6ycc84555xzrpC8YOWcc84555xzheQFK1dsiIiKyCEhpd1LRF4JI+1ccUwUkSv38tivRaR1kkNyzhVBIvK7iDQIO46iTkSGicjDCe676/orIl1EZGpqowuX5xnh8f/f4ssLVq5EEJEpInJ1qvYvKlT1TFUdvpfHNlLVKQU9TkQOEJHRIrJURNaLyKci0iLXPpeJyEIR2Sgib4lI1aht5URkqIj8JiLLReS2qG3Vg/OtEZF1IvK5iJyURyyHicjbIrJKRH4VkfdF5PBc+9wapLM+SLdc1LYbRWSGiGwRkWG5jisrIq+LyIKgEN+6oJ+Vc0WFqu6jqj8nsm+YD62Kk8Jcf8MkIj+IyGEFPW5v84wgzRLznSpI4TtZCvL/Gy0JeWRVEXkzyMsXishlUduOCvLPtcFrkogclUcsJ4jIh0Ecq0TkNRGpFbVdROTRIP9fIyKPiYhEbX9IROaIyHYR6ZXr3LVEZHxwX6Iikl3QzypVvGDlXJqJSKmwYyigfYAvgWZAVWA48K6I7AMgIo2AF4DLgZrAJmBg1PG9gEOBesCpwJ0i0i7Y9jvwN6AGsD/wKPCOiJSOE0sVYDxweJDWdODtyEYROQPoCbQBsoEGwANRxy8FHgaGxjn/VKAzsDzOdudclDz+V10RICINgSxV/SHsWIqidHx/0/w/Utg88llga3BsJ+C5II8Hyz8vxu4DqgfpjMkjlv2BQUE69YANwEtR27sDFwDHAE2Ac4Brorb/BNwJvBvj3DuB94C/5pF+OFTVX/7a7YX9083D/gm+AS6M2tYFu/l8AlgLzAfOjNo+BXgI+DQ4/gOgerCtNbAkV1oLgLbB78cDnwPrgGXAM0DZqH0VOCRGvL2BHcBm7Eb9mWD9n7ECwfrg55/z2b8fsBj4DZgJtIxKoxfwSpzPqzWwBLgHWB28p05R24cBzwETgI1AW6Bc8BkuAlYAzwMVoo45H5gdxDIPaBf1+V4d9bf4FBgQvMfvgDZ5/F2jP+tewKvAy8Hf6WugeQG+I78BzYLf+wCjorY1xC7MlYPlX4DTo7Y/BIyJcc4s4Nzg73xAgnFUDfavFiyPAvpEbW8DLI9x3MPAsDzOuwRoHfb/or/8Ff0CugLvRC3/BLwatbwYaBr8vut6GVyDnsVuUDYAXwANg23/CfbdGFwPOwTrzwmuQeuAz4AmUeksAO4C/gdsAUrHiFWB64EfgzQfCq4NnwfXj1fZ/fqeV3p7nSfFiOtPwH+Dc43FbgwfDrbtD/wLWBWc619Anahjp7D79Xdq8PuzwJO50nkH6BEnhoQ/m7xiwq5/S4Bzg+V9gu/EFVFp3Qz0j/oeDAQmBn/rT4EDgaeDc38H/CnX37nAeQZ7/536O/ad2ggMwW7uJwbpTQL2D/bNDs7fHbvZXwbcHnWurKjvzJog7qq5jr0Ky3//E6x/DXugtj6Iv1GwvjuwDcvTfif4/yPX/Ujw2Ua+R62Dv8tdwTlH5Pf+43xH8v3/TeCakXAeCVQK3udhUdtHAH1jnLc0cAOwqQDXr2OBDVHLnwHdo5avAqbFOO4VoFecc5YO3l92Qa+nqXqFHoC/it4LaA8cFFycOgQXuVrBti7BRaYbUAq4LriwSbB9SnAxOwyoECz3Dba1Ju+CVTPghOAfJRv4lqiMKfeFLNd5phBkeMFyVSyjuDw436XBcrVY+wfrOgPVgv1vDy6I5YNtvci7YLUd+D+swNQq+MwOD7YPwy7WJwWfaXksIxsfxFkZy4QfCfY/Ptj/L8H+tYEjcscd/C22A7cCZYK/1XqCDCRGnNGfdS+sYHlW8Hd8hBgXtDjnaRocu1+w/DZwV659fg/+nvsHf7eaUdsuBubk2v9/2AVdgcEF+K5eACyLWv6KIBMPlqsTlalErfeClb+K3Qt7urwuuC7UAhYCv0RtW4vVTsCeN2a/BteW0sBIoh5usOdN4rHASqBFcH24Mrh+lAu2L8BuEA8m6oFQrlg1uMbtCzTCCmD/DuLcDysgXZlgenudJ+WKqWzwmUWumRcHx0ZuiKthT8ArYtfl14C3oo6fQuyC1fFBmpHPvjpWc18zCZ9NfjGdjuVVBwCDgddzpfUecEbU92A1dm0uD3yEFUSvCD67h4HJUccuYC/zDPbuOzUNK0zVDvb9L1YQLhfEen+wb3Zw/tFYYaAxVvCMxNojOFed4NgXgNG5jn05OLZCsP5vwedbDsufZ0fFPozgO5LH+9u1Dzn3BI8G56uQ3/vP6/Mjn//ffK4ZCeeRwWf9R67j7yDqYU6wbl3w/nYC/yjA9atH9HcGu19pEbXcnKiCV9T6YlWw8qaAbg+q+pqqLlXVnao6FnuqdnzULgtVdbCq7sCahdXCLoYRL6nqD6r6B/akqGmC6c5U1Wmqul1VF2AXw1Z7+TbOBn5U1RHB+UZjT+POzSP9V1R1TbD/k9gF8fB4+8dwn6puUdWPsSdLl0Rte1tVP1XVnVgm2g24VVV/VdUNWK1Px2Dfq4Chqvph8Df4RVW/i5PmSuBpVd0W/K2+D957Iqaq6oTg7zgCq47Pk4jsG+z7gKquD1bvg10go63HMql9opZzb9tFVZtgNxmXYU+f8yUidbCneLdFrc4dS+T33dJzrjhS63OxAbumtgLeB34RkSOC5U+Ca0ws41R1uqpux27M8roudwNeUNUvVHWHWr+iLdiDr4j+qro4uM7H86iq/qaqXwNzgQ9U9efg2jERu5HLN70k5EkRJ2AFqsg183WsNQNBOmtU9Q1V3RRcl3uTQB6kqtOxa02bYFVHYIqqrijsZ5NfTKr6AVbY+jd27d/VlEpEKgLHAR9HpftmkNduBt4ENqvqy8FnN5acv0ksBc4zoiTynRqgqitU9RfgE+ALVZ2lqluCWHPH9oCqblTVOVgTs0uD9dcA96rqkuDYXsDFuZrk9QqO/QNAVYeq6oao/Y8Rkf0K8P5y24kVBLcEaSTy/vNSkP9fYK/yyLzy8l1UtQr2AOBGYFYiwYtIE+CfWK1kXrHsE93PqjjygpXbg4hcISKzg8EE1gFHY081Inb1P1HVTcGv+8Tajj21i96WV7qHici/gk6Vv2GFjer5HRfHQdiTyWgLsSdh8dK/XUS+DTp0rsMuHImmv1ZVN+ZK66Co5cVRv9fAnj7OjPqM3wvWgz0Fnpdgur+o2mObOOnmJfffqXxebcFFpAJWszZNVR+J2vQ7ViiKti92A/h71HLubbtR1c1BAbiniBwTpPl71KtuVCw1sGamA4Nj4sUS+X2P9Jwrpj7GnoifEvw+BbvRbsXuN9C5FeS6XA+4PXJ9Cq5RBxP/mhZPdMHijxjLkRjyTC8JeVLEQcS+ZhKkU1FEXgg67f+GNQmrkmC/2OFYqweCnyPy2T+hzybBmAZhn8lLqroman0b4LOgEFWgdOMoUJ6RSyLfqYLGFv0djM776gFvRqXzLdb8v2asY0WklIj0FZF5wWe8INi0t/cfAKtyfe6JvP+8FOi+ai/zyLzy8t0E9zvPAy+LDXBVNzq/zhXLIdjDgltU9ZN8Yvk91/9nseMFK7cbEamHNSe4EWs+VQV7mpaMJwgbsQJFJK1S5BQmwPohfQccqqr7Yn2WEk039z/iUuxCFq0u1t9nj/1FpCXWHvoSrB13FezpSaLp7y8ilXKltTROfKuxTKKRqlYJXvupauRCuRhrb5+I2rme7uRONynERg16C/v8rsm1+WuinlyKDRFbDvhBVddi7d+jn2weExwTTxmsSQxqIyNFXouC8++PZRjjVbV3XrEEv6/IdbPhXHEWKVi1DH7/mMQKVgWxGOgddX2qoqoVc92gJfPmJ256Sc6TlhH7mhlxO9ZKoUWQB50SrE8krVeA84OHQkdi18tkyDOmIB99AWvadp3sPhLfWcTu+B+GRL5TBXVw1O/Red9irJ9ddFrlg5qwiOjv72VYv+a22APV7GC9xNg3YhNR9zNYX7VouY9JxfuPqRB55A9AaRE5NNf2ePl1FvYZ1FbVRdH5dVQs9bD+cQ+pau6HDbFiyeveoFjwgpXLrRJ2QVgFICJdsSdhyfAD9oTrbBEpA/wDuwGPqIx13P09aNpyXQHOvYLgZjwwAThMbBjw0iLSATgK6/gba//KWJvhVdiF5Z/s+eQmPw+IDdndEuuk+lqsndSa6gwGnhKRAwBEpHYwWg9Yp92uItJGRLKCbUfESfMA4GYRKSMi7bEMfUIB485T8Ld6HSsMXqF7NjUaCZwrIi2DwuWDWLOFyFOul4F/iMj+wfvohrUZjwzHenLwuVUQkbuwp4pfxIllX6z506eq2jPGLi8DV4kNC7s/9h0bFnV8aREpj7VxLyUiuz1xFRsavnywWDbYXqybJbgS52NsdM0KqroEazLVDusjkVCznBhyXw8HA9eKSAsxlYLrdqqa1OaVXjLzpM+x6/zNwbXgInZvUlgZu86tE5sy4v5ETxz8Lb7Eaqre0LybSBZEfjHdE/z8GzaAx8tRtVlnkuT8oADS8Z26L6jRa4QN7DI2WP880Du4qUdEaojI+XmcpzLWLG8NVlDok897AetjeFlQ29WO/JuMpuV/qjB5ZFALNQ54MIjvJKzAOSI4919E5E/Be94X61e+FqsRjBVLbaxv3LOq+nycWG4L7nEOwh4iDIs6vkyQH2dh92Xlo77bBNsi95DReXeovGDldqOq3wBPYhnQCqxT6KdJOvd6bCSkF7Gaj43YIAERd2BPjjZgF6Gxuc+Rh35YG+q1ItI/ePpyDvaPugYbsvMcVV0da3/sQjQRK/wtxDrpJtLUJWI5doFZihU0rtX4/aLAasd+AqaJNT2YRNCfK2iv3xV4Cqs1+5g9a98ivsCGMl+Ntb2/OAW1M3/GPsvTscw9Ut3fMoj3a+Ba7H2vxDKp66OOvx9r2rgweC+Pq+p7wbZyWBvwNdh34izgbFWNV+t2IdZnoKvEaCYYnPcxYHKQ3kJ2vxH5B3aT0hNrrvNHsC7i+2Bdbew78QfxP3vn0k5t2OzfsQIVqvob8DN2I7VjL0/bCxgu1kTpElWdgT0AeQa7rv2EDdaQEnmll8w8SVW3AhcF516LDYQxLmqXp7GBBlZjgx+8R8EMD+LLrxlgQcSNSUSaYf1nrgj+9o9ihdCeInI01qxqURJjKYhepP479XFwnn8DT6j1NwPL38cDH4jIBuxzaxH7FIDd4C/E8qBvgv2jDQGOCt5LpCbyFqzP9jpsWPI8ayjT+D9V2Dzyeuz7thIbHOS6II8HG8p9NHZfMg84BBuxOLrJY7SrsQLp/RK7meALWPeCOVgt9LvBuojBWB58KXBv8PvlUdv/IKe7wXfBcugiI7k55/aS2ESyr6hqnTSn2wUboerkdKbrnHNuTyJyCtYkMDtGzX66Y7kTm+rkzjDjSAWxyWDnA2WCwRycKzJ8Yj/nnHPOuUIImkzfArwYdqEqsACrDXDOpZEXrJxzzjnn9pKIHAnMwOYI6hpyOACo6qthx+BcJvKmgM4555xzzjlXSD54hXPOOeecc84VUoltCli9enXNzs4OOwznnMtoM2fOXK2qNfLfM/N4PuWcc+FLZj5VYgtW2dnZzJgxI+wwnHMuo4nIwrBjKKo8n3LOufAlM5/ypoDOOeecc845V0hesHLOOeecc865QvKClXPOOeecc84VUontYxXLtm3bWLJkCZs3bw47lJQpX748derUoUyZMmGH4pxzzjnnXMbIqILVkiVLqFy5MtnZ2YhI2OEknaqyZs0alixZQv369cMOxznnnHPOuYyRUU0BN2/eTLVq1UpkoQpARKhWrVqJrpFzzhUPI0eCjSTerFnIobhMFPkCZmXZz5Ejw47IOZcBMqrGCiixhaqIkv7+nHNF38iR0L07bNoUdiQuI+X+Ai5caMsAnTqFF5dzrsTLqBor55xzqXfvvV6ociGK9QXctMnWO+dcCnnByjnnXFItWhR2BC5jTZtmNVSxLFxotVnxtjvnXCF5wSoPYTXR3rFjR3oScs65FKhdO+wIXEbZvh3GjoUTToATT4R4TeJFoHNny9Dr1rVmgc8/D19/DTt3pjVk51zJ5AWrOCJNtBcuBNWcJtqFLVwtWLCAI444giuvvJImTZpw8cUXs2nTJrKzs3nwwQc5+eSTee211/jpp59o27YtxxxzDMceeyzz5s1jypQpnHLKKVx44YUcddRRXHvttez0zMA5V8T4oKQuLdatg8cfhwYNoGNHWLMGnnkGXnwRKlbcfd+KFWH4cJg1C/r3t0LYRx/BddfB0UfDAQfA+efDE0/AF1/Atm3hvCfnXLGWcYNXRPToAbNnx98+bRps2bL7uk2b4KqrYPDg2Mc0bQpPP51/2t9//z1DhgzhpJNO4m9/+xsDBw4EbA6qqVOnAtCiRQt69uzJhRdeyObNm9m5cyeLFy9m+vTpfPPNN9SrV4927doxbtw4Lr744kTesnPOpdy778Inn9g96uzZ3urKpcBPP0G/fvDSS7BxI5x6Kjz7LJx9tjUxAShXzvpULVpktVO9e+cMXNG0Kdx0kz01/fln+8JGXuPH2z4VK1rhq2VLe51wAlSqFM77dc4VGxlbsMpP7kJVfusL4uCDD+akk04CoHPnzvTv3x+ADh06ALBhwwZ++eUXLrzwQsAKXBHHH388DRo0AODSSy9l6tSpXrByzhUJa9dCt27QuLG1zCpXDkRmzgw7LlcCqMKUKfDUU/Cvf0GZMnDppfaUtGnTPffv1Cn/EQBFoGFDe3XpYuuWL4epU3MKWg89ZM0ES5eGY4/NKWidfDJUq5bsd+mcK+YytmCVX81SdnbsJ6316tm1vTByD4keWa4UPA1T1QIf65xzYbvlFli50u57y5ULOxpXImzZAmPGWKY9ezbUqAH33WdN+A48MPnpHXggXHyxvQDWr4fPP88paD3zDDz5pG076qjdC1r16iU/HudcseJ9rOLo3Tt2E+3evQt/7kWLFvH5558DMHr0aE4++eTdtu+7777UqVOHt956C4AtW7awKRg6dvr06cyfP5+dO3cyduzYPY51zrkwvPMOjBgB99xjD/adK5RVq+DBB62w0qWL9Xl68UVr2vfAA6kpVMWy337Qrp1l/v/5j/Xr+uQT6NPHYhs92gfEcM7t4gWrODp1gkGD7LopYj8HDUrO3IJHHnkkw4cPp0mTJvz6669cd911e+wzYsQI+vfvT5MmTfjzn//M8uXLATjxxBPp2bMnRx99NPXr19/VXNA558Ly6682uE+TJvCPf4QdjSvW5s6Fq6+Ggw+G+++HZs3ggw9gzhzr5BzVND4U5ctb7dTdd8OECfbl9wExnHOBjG0KmIhEmmjvjaysLJ5//vnd1i1YsGC35UMPPZSPPvpot3WLFi2iYsWKjB07NvlBOefcXrr5Zli9GiZOhLJlw47GFTs7d8J771lzvw8/hAoVoGtXa1t6xBFhR5e3UqWsj1f0gBjz5hV8QIyRI+MPtuGcKza8xso559xee/vtnHvCWGMIuDzMnJneSRKLmk2brNlco0Y2ot/XX1sTu8WL4bnnin6hKhYROOQQKxgOHQo//gjLlsGrr1pN3Nq1NiBG27ZQpQq0aGHv/aqrkj+/i3PpFtYEsEWI5DVQQnHWvHlznTFjxm7rvv32W4488siQIkqfTHmfzrlwrVlj98S1allLp1i1VSIyU1Wbpz+6oq+5iM4Aq8lIVlvz4uCXX2x49BdesKZ0zZrBrbdC+/aZUeWZe0CMTz6JvV/duj5fgSs+IhPABmMCAMXm2pbMfCplNVYiMlREVorI3Kh1VUXkQxH5Mfi5f9S2u0XkJxH5XkTOiFrfTETmBNv6SyGHwSupBcmIkv7+nHNFx003WeFq2LDMuB9OmU2brMqvpJsxI2egh0cfhdatrVDx5Zd245UpX6LcA2LEu61ZtMhGJxw6FJYuTW+MzhWEKvz977sXqiBzrm1RUtkUcBjQLte6nsC/VfVQ4N/BMiJyFNARaBQcM1BESgXHPAd0Bw4NXrnPmbDy5cuzZs2aElv4UFXWrFmz27xXzjmXCm++aQOi3XcfHHNM2NGUAAsXWkn1+eetsLFmTdgRJceOHTBunPUpOu446290443WRO6NN2wgiEyfNqRu3djr99kHpk2zZoK1a1tb27vvtu/H9u3pjdG53BYtskm6O3e27+eyZfH3yyApbQooItnAv1T16GD5e6C1qi4TkVrAFFU9XETuBlDVR4L93gd6AQuAyap6RLD+0uD4a/JLO1ZTwG3btrFkyRI2b96cnDdYBJUvX546depQpkyZsENxzpVQq1dbE8Data0JYF6Xm0xoCigijwPnAluBeUBXVV2X33G7mgKC1daUKwcbNuTsULOmfdBHHWU/I78Xh4lpf/vNalr694f5862W6uabrZCw775hR1e05NWE6rLLbKTECRNsdJhPP7VC1X77wV/+AmeeabVfBx0UXvwuM6xZA5Mnw7//ba8ff7T1NWpAmzY2euevv+55nAhcc409ODrqqPTGnKBk5lPpHhWwpqouAwgKVwcE62sD06L2WxKs2xb8nnt9TCLSHavdom6MJ0BlypShfv36hYnfOecy3o03Wh/8SZPyLlRlkA+Bu1V1u4g8CtwN3JXw0dE30YsXwzff2EAOX39tvw8fXnwKXPPnW2FqyBCL+eSTbbjx88+3EfTcniL9T+KNCti4sb3uusv6Z02aZIWsiRPh9ddtn6ZNrZB11lk22mBpH/TZFdLGjTB1qhWiJk2yCbpVrSa1VSubVqBNG5taICsr9gOC8uVtgJaXXrLa+LZt7QHLWWeV2OtBumus1qlqlajta1V1fxF5FvhcVV8J1g8BJgCLgEdUtW2wviVwp6qem1/asWqsnHPOFc4bb1i3j4ceSmzOqkyosYomIhcCF6tqvr21m4vojHr18h9aWzV2gevrr+H333P2ixS4che6qlZNwjvLg6rVpDz1FLz1lt1kXXKJDUjRPGP+9OmnCv/7X04h69NPrelllSpWm3XWWVabla7JlF3xtm2b9XecNMkKU59/buvKlIETT7RCVNu21qQ33hO1eNMGrF4NgwfboDW//AINGtgTuq5d7fsasmTmUxnVFNA559zeW7XK7tPr1rU8N5HaqgwsWL0DjI08KIyxPbplRbOFhRn1LVLgii5oRX5PR4Fr61Z47TWbf2rGDNh/f2vyc8MNUKdO4c7tCm7dut1rsyJ9Xv70J6vNOvNMr81yOVRt4u1I076PP7brhojVgLZta4Wpk0/OmW+tsLZts4cv/ftbbVilSnDFFdZMMMTRrItzwepxYI2q9hWRnkBVVb1TRBoBo4DjgYOwgS0OVdUdIvIlcBPwBVaLNUBVJ+SXthesnHMuuS65xOatmjnTWn8koqQUrERkEhDr0f+9qvp2sM+9QHPgIk0gc01ZPpXsAlfup9D33GP9LZ55xkarO/xw6NEDLr88eTdgrnBU4auvcgpZn32WU5t1+uk5fbO8NiuzzJ+fU5D66CNYudLWH3qoFaLatIFTT01Ps+JZs2DAABg1CrZssVrWSDPBrPROs5vUfEpVU/ICRgPLyOkndRVQDSs0/Rj8rBq1/71Yp9/vgTOj1jcH5gbbniEoDOb3atasmTrnnEuOV19VBdXevQt2HDBDU5TPFKUXcCXwOVAx0WP22aeZvvSSfU5bt6q2aqU6YoQtb9xoy2PG2PK6dbb8xhu2vGqVLY8fb8vLltnyxIm2vGiRLX/4oS3Pm6faqtVOnTJ2ueqECfrdnUO01YHf6qdHXqW6zz46h0baisk6neaqNWvqrOO6aavaP+is8/6pWq6cTqe5tmKyzqGRKuinnKitqszW716Yorpjh06ZYunNm2fpffihLS9aZMsTJ9rysmW2PH68La9aZctvvGHL69bZ8pgxtrxxoy2PGGHLW7fa8ksv2XLEoEGqbdrkLD/7rGq7djnLTz+teu65OcuPP6560UU5y488otqhQ87ygw+qduqUs3zffapduuQs9+yp2q1bzvLtt6tef33O8i232Cvi+uttn4hu3ewcEV26WBoRnTpZDBEdOliMERddZO8h4txz7T1GtGtnn0FEm1ZbdVD3L1W7dlU98EBtxWR9iStVjz1Wt/a8T1s1Xacjhm1X1VR991SnTLHl776z5U8/teU5c2x5+nRbnjXLlmfNsuXp0215zhxb/vRTW/7uO1v2714e372VK/X603/Q2xu/r9qggSpoN17QnpX6WyJDh2qX9htS+91rY59RRKtWuvt1789bdET7t1Rr19aNVNBW5afpmMv/pbp2bdq+e8nMp1JWH6yql8bZ1CbO/r2B3jHWzwASfDbqnHMu2VauhOuvt+4yd94ZdjRFj4i0wwaraKWqm/LbPzxitVStzoQGWDuQPi/CiYNh0nK4vSycegts+Aimb4Jly2H828CWPU9VrTocfQy0IrUTt7jkKF3G/oG7DYWdO+H4jZC9AVbMsznF9DS4fjRMXA+nnQNbLwEyZF6xkuT332HBKli0GI65yfrg8QyU3QntjoZbboFP2kPDKtD3JjvmP6FGDGXKwnnnw8izYMx4uLEsjHgZxnWAjtfApgeB4lMTntKmgGHypoDOOVd4qtC+PbzzjrXcKOhouSWlKWBeROQnoBwQmXxqmqpem99xRT6fUrWRu2LdJ4jYDbor/tauhQ8/zGk2uGKFrW/WLKdvVosWJXYUt2Jt61ab8yIyct8XX9hw/OXKwUkn5TTva9asePWtmznTmgmOHm3ApqAIAAAgAElEQVTv8fTTrZngmWempJlgMvMpf87knHMurldftZEAH3igyE5BEjpVPURVD1bVpsEr30JVsSASf/LaeOtd8bP//taB8qWXrM/czJnw8MN2c96nj92gH3AAXHopvPyyVWGPHGlzk2Vl2c+RI8N+FyVPrM94504b9vyJJ6yQUbUqnHIKPPigFUDuuMMKyWvXWmHrnnusUFycClVgBcFhw6yv6MMP21xu55wDhx1mg+WsXx92hHF5jZVzzrmYVqywsQwaNrSRnPcmb86EGqu9VSzyqbwmr81riHhXMvz66+61WZHBDrKydq+x9O9EcsX6vytVCipUyBmA5ogjcmqkWre2AnJJtW0bjBtnowl+9pkNktOliw3ZfsQRhT59sRkVMEzFIsNyzrkiShX++leYMMGaAO7tSLhesIqv2ORT8eamcZll5067GLRpE7vGoFYtm6NIJP2xlTQHHZQzXH60SpVg4ED7G9Sunf64ioLczQTPOMOaCbZrt9fNBL0poHPOuZQaMwbefNNamIQ4vYgrCjp1ggUL7MZ6wQIvVGWqrCxrovXbb7G3L1tm7YUfeAB++CG9sZUE8+ZZ08smTWIXqsBqsK64InMLVbB7M8GHHrIBOs4+26Z96N8//vczTbxg5ZxzbjfLl1sLixNOgNtvDzsa51yREq9/XdWq1hfrgQfsJrdZM+sLtHhxeuMrThYtss/ouOPgkEOsZrhy5fjN+rxvY44DDoB//MMe9oweDTVq2KiHtWvbhMPffx9KWF6wcs45t4sqXHstbNxofdl9IDDn3G5697Y+VdEqVrTago8/tsLCk09aDdff/26FgZYtrQnbqlXhxFyULFtmn9VJJ0G9evYZATz+OCxcaB1aBwyI/Rn33mNWIle2LHTsaH2vvvwSLrrI+vsdcYQN8DFxYlpHMPWClXPOuV1GjYK337aBmJLQJ9g5V9J06mQ3rvXqWX+qevV2H7iiTh247Ta7yf3hB2uu9euvcMMN1g/rjDOsKVcRHtkt6VavhhdegFNPtRqVW26xQSh694Yff7TP6o47cmqk8vuMXWzNm8Pw4Va4f/BB+OorOOssy8zS1EzQB69wzjkH2IPURo0sD/rkk+TUVvngFfF5PuUyhirMmWOdN0ePtuZbZcvaTe+ll9pQ2rlraIq7deuso+rYsTbH1I4d1kSyY0fo0ME7r6bD1q02X0j//jBtGuyzD3Ttam3dDzts18A8zRcuZIZqUkZd8YKVc845VOH8821k5a++sjwnGbxgFZ/nUy4jqdpEtmPGWKFj+XIb7e78862QdfrpVugqjjZsgPHj7X29954NE16/vhWkOna0gSl81MRwTJ9uTSzHjrW/S5Mm1g9ryxaaQ9IKVt4U0DnnHK+8Au+8Yy1TklWocs65PYjYyDhPPw1LlthEtpddZn1hzj0XDjwQunWDjz6yWp6ibtMmeP11uPhiG1Chc2cblv6mm+xmft48eOQROOYYL1SF6fjjYcQIayb4wAM26fCWLUlPxmusnHMuwy1dak0AGzWyvufJHLDCa6zi83zKuShbt8IHH1hN1ltv2Qg6Bx4Il1xiNVktWhSdgsmWLfD++1b78fbbFmvNmtC+vdVO/fnPez2nkkuTrCyrPYWk1liVTsZJnHPOFU+q0L273Sf4KIDOudCULWt9rc45x2qB/vUvK2Q9/7z1kalf35rTdewIjRunv5C1bZvVro0da32n1q+3IeYvu8xiatXKL6DFSd26Ngpjknlx2jnnMtjLL8O779q8lIceGnY0zjmHDWRxySUwbhysXGlPfQ47DB57zJrUHX20DV3600+pjWPHDpg8Ga65xkY0PPNMi+mCC2DCBOsfNmgQnHaaF6qKm1jTBiSBNwV0zrkM9csv1vyvSROYMiU1LVe8KWB8nk85V0ArV1p/ptGjYepUW9e8uTUV7NDBhjIvrJ074fPPrbbs9ddzBtc47zxL44wzoHz5wqfjwuejAibOMyznnItPFc4+2wpU//sfHHJIatLxglV8nk85VwiLF1uzvNGj4b//taaBLVtaIevii6F69cTPpQozZtj5xo61QTXKlbOLZMeO9rOkDQfvdklmPuVNAZ1zLgMNG2aDcPXtm7pClXPOpczBB9ukujNn2rDZvXpZjdZ119mgF2eeaW2dI5PCjhwJ2dlWNZ+dbctffQX33GMXweOPt75cTZva6HErV9ocSO3be6HKJcxrrJxzLsMsWWJNAJs2te4DqRy8ymus4vN8yrkkU7XC0pgx9lq40GqeGje2qvmtW3P2FbH9S5WCNm2smd+FF8L++4cXvwtFMvMpHxXQOecyiCpcfTVs3279wX1EYOdciSFiT4yaNrW5oyJ9pZ591vpORVO1Uf2++w5q1AgnXlfieJbqnHMZZOhQm37l0UehQYOwo3HOuRQRsfmk+vffNV/RHtau9UKVSyovWDnnXIZYtAhuuw1at4brrw87GuecS5O6dQu23rm95AUr55zLAJEmgDt2WK2VNwFMPhG5Q0RURAowHJlzLuVizVlUsaKtdy6JPGt1zrkM8OKL8OGHNr9m/fphR1PyiMjBwF+ARWHH4pzLpVMnm8i3Xj1rIlivni136hR2ZK6ECaVgJSK3iMhcEflaRHoE63qJyC8iMjt4nRW1/90i8pOIfC8iZ4QRs3POFVcLF8Ltt8Npp8G114YdTYn1FHAnUDKH2nWuuOvUCRYssEEsFizwQpVLibSPCigiRwPdgOOBrcB7IvJusPkpVX0i1/5HAR2BRsBBwCQROUxVd6QxbOecK5YiTQBVYcgQbwKYCiJyHvCLqn4lIvnt2x3oDlDX+3c451yJEsZw60cC01R1E4CIfAxcmMf+5wNjVHULMF9EfsIKZZ+nPFLnnCvmBg2CSZPguedsTky3d0RkEnBgjE33AvcApydyHlUdBAwCm8cqaQE655wLXRjPLucCp4hINRGpCJwFHBxsu1FE/iciQ0UkMkNbbWBx1PFLgnV7EJHuIjJDRGasWrUqVfE751yxsGAB3HEHtG0L11wTdjTFm6q2VdWjc7+An4H6wFcisgCoA/xXRGIVwpxzzpVgaS9Yqeq3wKPAh8B7wFfAduA5oCHQFFgGPBkcEqtdRcynfKo6SFWbq2rzGj4vgXMug+3cCVddZf20X3zRfrrkU9U5qnqAqmarajb28O9YVV0ecmjOOefSLJTW9qo6RFWPVdVTgF+BH1V1haruUNWdwGCsuR9YJnVw1OF1gKXpjdg554qXF16Ajz6CJ56wAbCcc845l1phjQp4QPCzLnARMFpEakXtciHWZBBgPNBRRMqJSH3gUGB6OuN1zrniZP58+Pvf4S9/gW7dwo4mswQ1V6vDjsM551z6hTF4BcAbIlIN2AbcoKprRWSEiDTFmvktAK4BUNWvReRV4BusyeANPiKgc87FtnMn/O1vNvqfNwF0zjnn0ieUgpWqtoyx7vI89u8N+PTYzjmXj+eegylTYPBg8NG8nXPOufTxGU2cc66E+PlnuPNOOOMMG7jCOeecc+njBSvnnCsBdu6Erl2hdGlvAuicc86FIaGmgMGcUgcBfwALgpH7nHPOFRHPPgv/+Q8MGQJ16oQdjXPOOZd54hasRGQ/4AbgUqAssAooD9QUkWnAQFWdnJYonXPOxfXTT9CzJ5x5ptVaOeeccy798qqxeh14GWipquuiN4hIM+ByEWmgqkNSGaBzzrn4IqMAliljA1Z4E0DnnHMuHHELVqr6lzy2zQRmpiQi55xzCRswAD75BF56CWrXDjsa55xzLnMl2seqCZAdvb+qjktRTM455xLw449w991w9tlw5ZVhR+Occ85ltnwLViIyFGgCfA1EBq1QwAtWzjkXkh07rD9VuXIwaJA3AXTOOefClkiN1QmqelTKI3HOOZew/v3h009h+HA46KCwoylaRKQW0AFoSc6ItnOBd4EPVFVDDM8551wJlcg8Vp+LiBesnHOuiPjhB7jnHjj3XLj88rCjKVpEZDDwCpa/9QO6ArcBU4ELgE9F5OTwInTOOVdSJVJjNRwrXC0HtgACqKo2SWlkzjnn9rBjB3TpAhUqwAsveBPAGJ5R1a9irJ8NvCoi5YG6aY7JOedcBkikYDUUuByYQ04fK+eccyF4+mn4/HMYMQJq1Qo7miLpFxE5XFW/j14pIkcAq1R1DfBDOKE555wryRJpCrhIVcer6nxVXRh5pTwy55xzAIwcCdnZkJUFd9wBxx4LnTqFHVWR1R/rV5VbfaxpoHPOOZcSiRSsvhORUSJyqYhcFHmlPDLnnHOMHAndu8PChRAZcuHbb2HUqHDjKsKOUdXJuVeq6kSgaQjxOOecyxCJFKwqYH2rTgfODV7npDIo55xz5t57YdOm3df98YetdzHl1cS9TNqicM45l3Hy7WOlql3TEYhzzrk9LVpUsPWOeSJyhqq+H71SRE4H5ocUk3POuQyQyATB9YGbgOzo/VX1vNSF5ZxzDqBuXWsGGGu9i+k24B0R+RiYGaxrDpyCtbhwzjnnUiKRpoBvAQuAAcCTUS/nnHMp9vDDew6pXrEi9O4dTjxFnap+BzQGvgCOCF5fAE2CbSkhIjeJyPci8rWIPJaqdJxzzhVdiQy3vllV+6c8Euecc3uoVMkGraheHdassZqq3r19VMC8qOpmYHC60hORU4HzscLbFhE5IF1pO+ecKzoSKVj1E5H7gQ+wQSwAUNX/piwq55xzqFohqkED+P57KJ3IFTvDichk4FXgbVVdGrW+NPBn4Epgqqq+lMRkrwP6quoWAFVdmcRzO+ecKyYSyaYbYxMEn0bOBMEaLDvnnEuRDz+EmTNh0CAvVBXA2cDVwJsiUhv4FRvdthzwb+BZVZ2R5DQPA1qKSG9gM3CHqn4Za0cR6Q50B6jrHeWcc65ESSSrvhBooKpbUx2Mc865HL17Q+3acMUVYUdSfKjqJmyS4P4iUg44APhDVVcX5rwiMgk4MMame7G8dH/gBOA44FURaaAamXlst/gGAYMAmjdvvsd255xzxVciBauvgCqAN21wzrk0mToV/vMfeOopKFcu7GiKp6Bp3uIknattvG0ich0wLihITReRnUB1YFUy0nbOOVc8JDIqYE3gOxF5X0TGR16FSVREbhGRucHoST2CdVVF5EMR+TH4uX/U/neLyE/BiEtnFCZt55wrDvr0sQErunULOxKXgLcImseLyGFAWaBQNWTOOeeKn0RqrO5PZoIicjTQDTge2Aq8JyLvBuv+rap9RaQn0BO4S0SOAjoCjYCDgEkicpiq7khmXM45V1TMmgUTJ9pQ65UqhR2NS8BQYKiIzMXytStjNQN0zjlXssUtWImIqPk4v30KmOaRwLSgHTzBJI4XYkPVtg72GQ5MAe4K1o8JmnTMF5GfsELZ5wVM1znnioU+fWDffeGGG8KOpHgTkTrAoao6OehvVVpVNyY7naAPcudkn9c551zxkldTwMnBhIe7DVskImVF5DQRGY4NW1tQc4FTRKSaiFQEzgIOBmqq6jKA4GdkHpDa7N5Gfkmwbg8i0l1EZojIjFWrvGm7c674+e47eOMNK1RVqRJ2NMWXiPwNGA+8GKyqB7wdXkTOOedKurwKVu2AHcBoEVkqIt+IyM/Aj8ClwFOqOqygCarqt8CjwIfAe9jgGNvzOERinSbOuQepanNVbV6jRo2Chuacc6Hr2xfKl4cePcKOpNi7GRul7zcAVf2BnAd2zjnnXNLFbQoYzFw/EBgoImWwEY7+UNV1hU1UVYcAQwBEpA9WC7VCRGqp6jIRqUXOKIRLsBqtiDrAUpxzroRZsABeecVqqw7wIkBhbVbVrSL2bE5EShH7QZ1zzjmXFImMCoiqblPVZckoVAGIyAHBz7rARcBorMlGpGnhleQ02RgPdBSRciJSHzgUmJ6MOJxzrih5/HHIyoI77gg7khLhUxG5EygvIqcCY4F/hRyTc865EiyRUQFT4Q0RqQZsA25Q1bUi0hebVPEqYBHQHkBVvxaRV4FvsCaDN/iIgM65kmbZMhgyxCYDPvjg/Pd3+boT6A58B9wCvA+8EGpEzjnnSrRQClaq2jLGujVAmzj79wZ6pzou55wLy1NPwbZt0LNn2JEUf0Gzv6GqeiXwXNjxOOdcUbNt2zaWLFnC5s2bww4lbcqXL0+dOnUoU6ZMytIIq8bKOedc4Ndf4bnnoEMHOOSQsKMp/lR1h4jUEpEyqrot7Hicc66oWbJkCZUrVyY7O5tIX9SSTFVZs2YNS5YsoX79+ilLJ9+ClYhchI3idwDW8VcsPt03ZVE551wGGTAAfv8d7r477EhKlJ+BT0TkbWDX3FWq2j+8kJxzrmjYvHlzxhSqAESEatWqkerpmBKpsXoMODcYJt0551wSbdgA/frBeedB48ZhR1OirMKm9agYvJxzzkXJlEJVRDrebyIFqxVeqHLOudR44QVYuxbuuSfsSEoWVb0v7Bicc87Ft27dOkaNGsX111+f8DHPP/88FStW5IorrkhhZHsvkYLVDBEZC7wFbImsVNVxKYvKOecywObN8OST0KYNtGgRdjQli4h8SIzJ5FX19BDCcc45l8u6desYOHBggQpW1157bQojKrxEClb7ApuA6MxIAS9YOedcIQwdCsuXw8iRYUdSIv0j6vfywF+JejjonHMuXD179mTevHk0bdqUMmXKULFiRWrWrMns2bO56KKLaNy4Mf369eOPP/7grbfeomHDhvTq1Yt99tmHO+64g9atW9OiRQsmT57MunXrGDJkCC1b7jHweFrlW7BS1a7pCMQ55zLJtm3w2GNwwglw6qlhR1PyqOoXuVZ9LCIfhxKMc84VYT16wOzZyT1n06bw9NN579O3b1/mzp3L7NmzmTJlChdccAHffvstVatWpUGDBlx99dVMnz6dfv36MWDAAJ6OccLt27czffp0JkyYwAMPPMCkSZOS+0YKKCu/HUSkjoi8KSIrRWSFiLwhInXSEZxzzpVUo0bBwoXWtyrD+g+nhYjsG/WqIiJtgFphx+Wccy624447jlq1alGuXDkaNmzI6adbY7nGjRuzYMGCmMdcdNFFADRr1izuPumUSFPAl4BRQPtguXOw7i+pCso550qyHTvgkUegSRM455ywoymxvsaarQuwHZgPdAs1IuecK4Lyq1lKl3Llyu36PSsra9dyVlYW27dvz/OYUqVKxd0nnRIpWNVQ1ZeiloeJSI9UBeSccyXdm2/C99/DmDFeW5VCDXJPDiwiieR5zjnn0qBy5cps2LAh7DCSKpFMZrWIdAZGB8uXAmtSF5JzzpVcqtCnDxx6KFx8cdjRlGhfAMfmWjc9xjrnnHMhqFatGieddBJHH300FSpUoGbNmmGHVGiJFKz+BjwDPIU1q/gsWOecc66A3nsPZs2CIUOgVKmwoyl5ROQArC9VBRFpjDUFBBvh1icKds65ImTUqFEx10+ZMmXX761bt6Z169YA9OrVK+Y+1atXLx59rFR1EXBeGmJxzrkSr08fOPhg6Nw57EhKrLOxh391gIFR6zcAPmmwc865lIlbsBKRO1X1MREZQOxJFm9OaWTOOVfC/Oc/MHUq9O8PZcuGHU3JFPQJfklELlHVV9ORpog0BZ7H5svaDlyvqtPTkbZzzrmiI68aq2+DnzPSEYhzzpV0ffpAjRpw1VVhR1LyqeqrInIG0Agr8ETW90lBco8BD6jqRBE5K1hunYJ0nHPOFWFxC1aq+k7w6yZVfS16m4i0j3GIc865OGbMgPfft2HWK3pPn5QTkYFAFeAUbIqQvwLTUpScYn24APYDlqYoHeecc0VYvhMEA3cnuM4551wcffrAfvvB9deHHUnGOFlVLwPWqOp9QAus31Uq9AAeF5HFwBPkkUeKSHcRmSEiM1atWpWicJxzzoUhrz5WZwJnAbVFpH/Upn2xNuTOOecS8M03NnfVP/4B++6b//4uKTZHforIgdg0Idl7ezIRmQQcGGPTvUAb4FZVfUNELgGGAG1jnUdVBwGDAJo3b75H/2XnnHPFV159rJZi/avOA2ZGrd8A3JrKoJxzriSJNP+75ZawI8koE0SkClaDNBvYAQzf25OpasyCEoCIvAxE/rqvAS/ubTrOOZcp1q1bx6hRo7i+BDXliNsUUFW/UtXhwCGqOjzqNU5V16YxRuecK7Z+/hlGj4ZrroHq1cOOJjOISBYwUVXXBX2E6wONVfWeFCW5FGgV/H4a8GOK0nHOuXCMHAnZ2ZCVZT9Hjiz0KdetW8fAgQPz37EYSaSPVbaIvC4i34jIz5FXyiNzzrkS4LHHbCLg228PO5LMoao7gX5Ry3+o6q8pTLIb8KSIfAX0AbqnMC3nnEuvkSOhe3dYuBBU7Wf37oUuXPXs2ZN58+bRtGlTjjvuOFq1asUll1zCYYcdRs+ePRk5ciTHH388jRs3Zt68eQC88847tGjRgj/96U+0bduWFStWAHDzzTfz4IMPAvD+++9zyimnsHPnzsK9772Q7wTB2GhK9wNPAacCXcmZyd4551wcS5fCSy9Bly5Qu3bY0WScD0XkfFV9O9UJqepUoFmq03HOuZTo0QNmz46/fdo02LJl93WbNtncIYMHxz6maVN4+uk8k+3bty9z585l9uzZTJkyhQsuuIBvv/2WqlWr0qBBA66++mqmT59Ov379GDBgAE8//TQnn3wy06ZNQ0R48cUXeeyxx3jyySfp27cvxx13HC1btuTmm29mwoQJZGUlUn+UXIkUrCqo6r9FRFR1IdBLRD7BClvOOefiePJJ2L4d7ror7Egy0o3AfiKyBfgDeyCoqlo13LCcc66YyV2oym/9XjruuOOoVasWAA0bNuT0008HoHHjxkyePBmAJUuW0KFDB5YtW8bWrVupX78+ABUrVmTw4MGccsopPPXUUzRs2DCpsSUqkYLV5qC9+o8iciPwC3BAYRIVkVuBq7G5P+ZgtWA9seYUkfFn71HVCcH+dwNXYZ2Pb1bV9wuTvnPOpdqaNfD883DppdCgQdjRZCTv0eacc4nIp2aJ7Gxr/pdbvXowZUrSwihXrtyu37OysnYtZ2VlsX27DUh+0003cdttt3HeeecxZcoUevXqteuYOXPmUK1aNZYuDW8qwUTqyHoAFYGbsaYOnYEr9jZBEakdnKu5qh4NlAI6BpufUtWmwStSqDoq2N4IaAcMFJFSe5u+c86lQ79+1lLibp/1LxSqugNoD9wV/F4LaBpuVM45Vwz17r3nzPYVK9r6QqhcuTIbNmwo0DHr16+ndtC2fvjwnIFeFy5cyJNPPsmsWbOYOHEiX3zxRaFi21sJDV6hqr+r6hJV7aqqfwXqFjLd0kAFESmNFdryKlqeD4xR1S2qOh/4CTi+kOk751zK/PYbDBgAF14IjRqFHU1mEpFnsH7BlwerNgHPhxeRc84VU506waBBVkMlYj8HDbL1hVCtWjVOOukkjj76aP7+978ndEyvXr1o3749LVu2pHow1K6qctVVV/HEE09w0EEHMWTIEK6++mo2b96cz9mST1Tznp9QRP6rqsfmt65AiYrcAvTG2r1/oKqdRKQX0AX4DZs/63ZVXRtkjtNU9ZXg2CHYMLqvxzhvd4LRmOrWrdtsYaxqS+ecS7FHH4WePeHLL6F587CjCZeIzFTVtH8KkXxKRGap6p+CdV+p6jHpjiWe5s2b64wZM8IOwzmXgb799luOPPLIsMNIu1jvO5n5VNwaKxE5U0QGALVFpH/UaxiwfW8TFJH9sVqo+sBBQCUR6Qw8BzTEmmosA56MHBLjNDFLg6o6SFWbq2rzRYtqJGuYfedcPlIwvUWx9ccf8H//B6ef7oWqkG0L+gcrgIhUA9I/9q5zzrmMkVdTwKVYzdFmYGbUazxwRiHSbAvMV9VVqroNGAf8WVVXqOqOYP6RweQ091sCHBx1fB3ybjq4S5KG2XfO5SFF01sUW0OGwMqVcE+qpqJ1iXoWeAOoISIPAFOBR8MNyTnnXEkWd1RAVf0K+EpERqrqXtdQxbAIOEFEKmJNAdsAM0SklqouC/a5EJgb/D4eGCUi/4fVcB0KTE80sU2b4N57C90M1DkXx7332v9ZtEz9v9u61SYEPukkOOWUsKPJbKr6sojMxB7mAbRX1bl5HeOcc84VRtyClYi8qqqXALNEZI+md6raZG8SVNUvROR14L9Yk8JZwCDgRRFpijXbWABcE+z/tYi8CnwT7H9DMMJTwhYt2ptInXOJiPf/lYn/dyNHwuLFNsy6+DTqRUEpYBuWr6R/pkjnnCvCVBXJoMwqv3ElkiHu4BWRGiQRqRdrezBZcJEl0lytJaO5/np44AGo7jObOJdUBx4IK1bsub5u3djTXpRUO3bAkUdCpUrw3/96wSoixMEr7gUuA97E+uqeD4xU1UfSHUs8PniFcy4s8+fPp3LlylSrVi0jCleqypo1a9iwYcOuSYUjkplP5dUUcFnwc6GIHIj1eVLgS1VdnozE06FCBTj5ZHjhBXuafN99cOONEDUHmXNuL6nCfvvFLlg1bWrbM+B6DcAbb8CPP8Krr2bOey7iOgPNVHUTgIj0xvoJF5mClXPOhaVOnTosWbKEVatWhR1K2pQvX546deqkNI24BasIEbka+CfwEfbUb4CIPKiqQ1MaWRLUq2dzl3XqBN98A3fcYa/nnoPHH4cLLvAbIOcKY+JE+OEH6NoVPvrImv8dfDA0bAjjx8NDD8E//xl2lKmnCn36wOGHw0UXhR2NCyxk9zyuNPBzSLE451yRUqZMmT1qblzhJdLm/O/An1S1i6peCTQD7kptWIXXrBksWJDTef6oo2DCBHjvPShf3m5+Tj3Vmuy4xPmw2i5i504boKJBA+tTtGCBrVu4ECZNgi5d4P77bU6nkm7CBPjqK7j7bihVKuxoXGAT8LWIvCgig4E5wDoR+b9gMCTnnHMuqfKtscKGO98QtbwBWJyacFLvjDOgTRt48eR4MrcAACAASURBVEVrFti8OVx5pdVsHXRQ2NEVbZFhtSMjwEWG1YbMG/3NwWuvwezZ8MorULbs7tuysux/bMsWmyi3XDno0SOcOFNN1a4f9erBZZeFHY2L8m7wipgWViDOOecyQ9zBK3btIPIy0Bh4G+tjdT423PkPAKpaJJ/8JdIpeP16uyHq1w/KlIG77oLbb4eKFdMUZDGTnR17MIJ69ay2wmWObdugUSOr/Z092wpSsWzfDh07Wv+jgQPhuuvSG2c6TJ4Mp50Gzz5rg+S43YU1eEVx4INXOOdc+JKZTyXSFHAe8BbB7PVYAWsZUDl4FVv77WdzznzzDbRrZ31BDj/cnsDv3Bl2dEXL5s3xR3jLxGG1M93w4TZQw8MPxy9UAZQuDaNGwbnnWqFjaJHvmVlwffpAzZrWz8wVHSLSTkS+FJGVIvKriKwVkV/Djss551zJlW/BSlUfyOuVjiD3yvffw7Bh9vu2bdC6tZWYwNqytW4NY8cC0LD6el5f3Zq5D42jZk3ocflqZlVpzTePvmP7L19u+7/3ni0vXmzLkybZ8s8/2/LHH+ek3bo1fPaZLc+da8tffmnLs2fb8uzZtvzll7Y8N5i78rPPbPn77235449t+eeg3/WkSba8OGiR+d57trw8GKzxnXdsefVqWx43zpbXr7flsWNtOdKm75VXbHnbNlseNsyWgTlzYEybwUyrHJljE65jIBM4c9fy/VX6wXnn5Xz2TzwBf/1rznLfvlZtEfHQQ9C5c87yP/+5+13p3XfntDEEG3Hkhhtylnv02L1d2Q032D4R3bvbOSK6dt19BIXOnS2GiI4dLcaIv/7V3kPEeedZtWbEmWda9UtE27YweHDOcuvWCX/3WL/elseNs+XVq235naL73dt5SmsG3beYE0+Ec8vk/90re3prXntxPWecAR9cNZYVR7ZO6LsH2OfaNue7x8CB9vlH9Av3u7fi4hs4Y9Id3H67jUDq373We373wvMMNh9ibaAGUD346ZxzzqVE3IKViDwd/HxHRMbnfqUvxPRpdBRMnw7PPmN9Q+7qCR065JRfMsUfm2HpMmjRApo0sfumatXsfjF3M8msLDj77HDidOH4ZSksW241NYmOqlmuHLz5pg0i8+13Offyxd3MmVC+HFx7bdiRuBiWALNVdZuq7oi8wg7KOedcyZXXBMHNVHWmiLSKtV1VP05pZIVU2LbrGzfakOyPPWbNAm+91QoW++6bxCCLEFX44gsbcGDMGHv/jRpBt272kL1aNdtv5EgbCW7RIqhSBdautaZel14abvwuPX77zUYBPPZY+OCDgh+/caM1u502DV5/Hc4/P/kxpsucOfbg4f77oVevsKMpukKcIPh44H5gCrAlsl5V+6c7lni8j5VzzoUvLX2sVHVm8OsM4BNV/TgoTE0FvkxG4kVZpUp2s/TDD1Zr1bcvHHooDBoEO0rQM881a6ylUZMmcOKJVqjq2BE+/9xuHG+5JadQBTb6X2RY7VWrbFTFW2+FdetCewsujZ56yr4zffrs3fGVKsG779p0CO3b2zDlxVXfvvZ+brop7Ejc/7d37/FWzekDxz9Pd6cokkSdYuSSmOhIjNG4dGGQaUJNohkqZISfSwklTrmEjHLJpXHPpZnRMGmIXBqkEHKrqHSZFFJJ1/P8/nj2nrM77XOps/Zea+3zvF+v/Tp7rX1ZT7t9znc96/v9Pt9S3ABsARpgQwCTtx0mImeIyBwRKRKRghKPDRaReSLyhYh0rsxxnHPOxVNFildMBVIHgO0EvJKZcKKnaVObqP/ee7D//tC/Pxx2WPE0gzgqKrLFXHv2tBLzl15qQ/zGjYNly6zXqn378od5Va8O999vCVbqtBKXm1auhNtvt2lABZW4rrPLLjZt55BDbD25OP4uzZtnFyEuvHDrCw8uUvZQ1dNUdYiqXpe8VfI9PwG6AW+k7hSRVkAP4GCgC3CPiPiKZs45V8VUJLGqo6prkxuJ+1WuIHlBAbzxhq3ds3YtdOxolc6Sc/zjYOlS62lo2dLW8nrpJUsUZ8+2YYB9+8LO21nn8fDD4ZJLLMF6x1eJyWkjR9pQvtTaCzuqQQMbSnjAAVaf4Y03yn9NlNx6qy3RcPnlYUfiyjBVRI4P8g1V9TNVTfdXvyswQVU3qOrXwDygXZDHds45F30VSax+EpHDkxsi0hb4OXMhRZcIdO9u5dlvucWKOrRubcPlvo9oEd/Nm63IV9eukJ9v86Py822u1NKl8Je/2DDAyhg+HPbe25K0ZHE3l1sWL7Z1ms45Bw46KJj3bNgQXn7Z1kf77W9t+GkcLF5shff+9Cdo0iTsaFwZ+gKviMjaLJRb3xtILXO0OLFvGyLST0RmisjMFStWZCgc55xzYahIYnUp8KyIvCkibwJPAxdnNqxoq1MHrrrKhgOdfz6MGQP77QejR8PGjWFHZ776Cq691hbvPe0065G68kqbM/baa/CHPyTKQwdg550tQfvoo60rQ7vcMXy4FTgJukjDHnvA1KmWoHTpUlwVPspGjbLhtFddFXYkrhy7AzWB+mxHuXUReUVEPklzK6vUSrqB02krQ6nqOFUtUNWCRo28+rtzzuWSiqxj9R5wIHAhcBFwUEphiyptjz3g3nttKN0RR1gRh9atYdIkOwnNtg0bbN7HiSfCL35hQ7cOO8zKXH/zjW23bJmZY59+ug2NHDrUFwzONV9+aQv7XnCBJepBa9LE5vw1bAidOxcvsRVFK1bYXMRevaynzUVXorT6GcDViftNgDYVeN2Jqto6ze35Ml62GGiWst0UWFqZ+J1zzsVPuYmViOQBVwMDVfVjoIWInJLxyGKkdWubr/Tii1bQoWtXS25mz87O8efMsaRur72sIMW8edbDsHAhvPCCJT01a2Y2BhG4+26771XScsv111sv7TXXZO4YTZtaclWvns1fTK5XHDWjR8P69TBoUNiRuPKIyBjgOKB3Ytc64L4MHW4S0ENEaovIPkBLYEaGjuWccy6iKjIUcDywETgqsb0YuCljEcWUCJx8sg2Hu/tuu+p+2GE2VPC//w3+eGvXWi/C0UdbYjd2rCVz//63DQO87jo7Wc2m5s3hhhusx+4f/8jusV1mfPghPP20VY5s3Dizx2rRwpKrWrXsuxy1wjA//mjDfrt1C26emcuoo1W1P7AeQFW/B2pV5g1F5HcishhrD18UkSmJ954DPAN8CrwEDPDFiJ1zruqpSGL1C1W9FdgEoKo/k348ucN6hi6+2HqNLrsMHn3Uht+NGAE/V7Lkh6rNQenXz4ZPnXeeLdB7++2wZImdAHfsCNUq8r+aIQMHWjGMP/8Z1qwJLw4XjCFDYNdd4YorsnO8/fazOVeqcPzxMH9+do5bEWPH2gLJmey5c4HaJCLVSMx1EpGGQFFl3lBV/66qTVW1tqo2VtXOKY8VquovVPUAVZ1cudCdc87FUUVOwTeKyE4UN06/IGUVe5ferrtawjNnjl19HzIEDjzQ5kBt7/yrH36wXrA2baBdO3j8catO+NZbVqHw8sshKnOga9a00utLlth8Kxdfb71lC/gOGmTl0bPlwAMtudqwwZKrhQuzd+zSrFtniyN36WJLDLjoEpEaibtjgYlAIxG5AVvc/pbQAnPOOZfzKpJYDcWGNjQTkSewBYO9HlYFtWxpxSNefRV2283mQB19dPlrPqnCtGlw9tnWO3XJJZa03HuvLeI7fjz86lflL+IbhvbtrfT6XXfBBx+EHY3bEaq26HOTJtYDm22tW1sp9tWrLblavDj7MaR64AFbIHnIkHDjcBUyA0BVHwWuBUYBPwBnqOqEMANzzjmX20TL6D4REcGqG60D2mNDAN9R1ZXZCW/HFRQU6MyZM8MOYytbttjQwGuusXlXPXtaEnLHHVZJLz/fSjivWQMPPmjDCevXt+TqvPNszlZcrFplPQ/5+bY+UfXqYUfktsfkyTZn8J574MILw4tjxgzr8W3SxNaN23PP7MewcSPsu69V2nz99ewfP+5EZJaqFmTxeB+oaiz+WkaxnXLOuaomyHaqzMQq5WBtgzhYNkW5wVq71hYYvvlmW8A3nWOPtcIXv/895OVlN76gPPWUrZc1ZgwMGBB2NK6iioqgbVvrLfrsMysmEabp060Me/Pm1oub7WGvDz4Iffta5c/Onct/vttaCInVYuCO0h5X1VIfy7Yot1POOVdVBNlOVWQo4DsickQQB3OmXj248cbSq6zttZddGe/dO75JFUCPHlZM45prYKmv6BIbzz5r1QBvuCH8pApsyOsLL8DXX9v36fvvs3fszZvtAkjbttCpU/aO6yqlOlAP2LmUm3POOZcRNcp/CscBF4jIAuAnbDigquqhO3pQEbkMOB8riPEx8EcgD3gaaAEsAM5U1R8Szx8MnAdsAS5R1Sk7euwoKS3ZWLYsu3FkiogNJWvd2iokPv102BG58mzebKX6W7e2oapR8ZvfwPPP2yLUnTpZcYv69TN/3GeftcqEEydGcz6jS2uZqg4POwjnnHNVT0V6rE4C9gWOB04FTkn83CEisjdwCVCgqq2xq4s9gEHAVFVtiRXIGJR4fqvE4wcDXYB7RCQnZuzk52/f/jjabz+49lp45hkbSuWi7a9/hblzobAwevPiOna0BOejj6w6X6bL+RcV2TIJBx1ki2y72PAU2DnnXChKTaxEpI6IXApciSU0S1R1YfJWyePWAHZKlMXNA5YCXYFHEo8/AiRPZboCE1R1g6p+DcwD2lXy+JFQWLjtUL+8PNufS6680gpZXHSRla120bR+vQ3/a9/eeoai6Le/tZ7P996z+z/9lLljvfACfPKJVUcMc204t91OCDsA55xzVVNZpwuPAAXYUL2TgNuDOKCqLsHK3y4ClgE/quq/gcaquizxnGXAHomX7A18k/IWixP7Yq9XLxg3zibli9jPceNsfy6pXRvuu8/myNx0U9jRuNLcc4+VNR8xItrD3n73O3jiCStq0bVr5RfeTkfVLnC0aBGtIZGufKqaxVl4zjnnXLGyEqtWqnq2qt4PdAd+HcQBRWRXrBdqH2AvoK6InF3WS9LsS1vKUET6ichMEZm5YsWKygebBb16wYIFNuxowYLcS6qSOnSAPn3gttts0WQXLatXW0LVsSMcd1zY0ZTvrLNs2OKrr1rlzA0BL1n+6qtW6v3qq6FGRWaiOuecc67KKyux2pS8o6qlFAXfIScCX6vqClXdBPwNOBpYLiJNABI/v008fzHQLOX1TbGhg9tQ1XGqWqCqBY2yXZPZleu222CXXeCCCyyRdNFx553w3XeWXMVF797Wwzt5siVamzaV/5qKKiy0tbP69AnuPZ1zzjmX28pKrH4pIqsTtzXAocn7IrK6EsdcBLQXkbzEAsQnAJ8Bk4BzE885F3g+cX8S0ENEaovIPkBLYEYlju9Csvvully99RaMHx92NC5p5Uq4/Xbo1g0KsrbaUDDOP9/WSXv+eVszrbR14bbH22/Da6/BFVdAnTqVfz/nnHPOVQ2lDnJR1YzUBFPVd0XkOeB9YDPwATAOW3fkGRE5D0u+zkg8f46IPAN8mnj+AFXdkonYXOb98Y82hOuqq+C007K/2Kvb1s03WxGIuM5/GzAANm6Eyy+3dbcefbRyFQ1HjIDddoN+/YKL0TnnnHO5L5TZA6o6FBhaYvcGSqnmpKqFQI7VyquaRKyQRZs21iPwyCPlv8ZlzuLF1uNzzjlWVjyuLrvM5lkNHmzFUh58cMcq+c2ebdUAhw+3hbydc8455yrKiwi7rGvVykqwP/qoDbly4Rk+3Oa7DS15mSOGBg2yf8f48daLpWlL3JRt5EjYeWe4+OLg43POOedcbvPEyoXi2mth332tkEXQFd1cxXz5JTz8sP0ftGgRdjTBGDrUEqz77rNerO1Jrr780hayvugi2HXXzMXonHPOudzkiZULxU472bpJX34Jt9wSdjRV09ChVpxhyJCwIwmOiM2RuvRSuOsuS7IqmlzdcosNI7zssszG6Jxzzrnc5ImVC03nzlYme8QImDs37Giqlg8/hAkTLAFp3DjsaIIlAnfcARdeCLfeCsOGlf+aRYtsaOr55+fe5+Gcc8657PDEyoXqzjutl+DCC3dsTozbMUOG2HC3K64IO5LMELGiHOedZ/PIylufa9Qo+3nllZmPzTnnnHO5yRMrF6omTaxgwNSp8OSTYUdTNbz1FvzrX3D11dCgQdjRZE61anD//baQ8JAhtlZXOsuXwwMP2PPy87Mbo4sXETlDROaISJGIFKTs7ygis0Tk48TP48OM0znnXDg8sXKh698f2rWzdYh++CHsaHKbKlxzDey5J/z5z2FHk3nVq1uBjjPPtN65MWO2fc7o0VZAZdCg7MfnYucToBvwRon9K4FTVfUQbIH7x7IdmHPOufB5YuVCV7269Sx8952f3GbaSy/Bm2/CdddBXl7Y0WRHjRrw+ONw+umWTD7wQPFjP/wAY8fCGWfA/vuHF6OLB1X9TFW/SLP/A1VdmticA9QRkdrZjc4551zYPLFykdCmDQwcCOPGwX/+E3Y0uamoyHqr9tnHijRUJTVrWrGOk0+2HtL+/a3E/G67wZo19v1zLiC/Bz5QVV9IwjnnqhhPrFxk3HADNGtmJ72bNoUdTe557jmrBjh8ONSqFXY02Ve7NkycaAtUjxsHCxcWP3bTTfDEE+HF5qJDRF4RkU/S3LpW4LUHA7cA/ct4Tj8RmSkiM1esWBFk6M4550LmiZWLjHr14O674ZNPrFqgC87mzTb87+CDoWfPsKMJT506sHr1tvvXrcut9bzcjlPVE1W1dZrb82W9TkSaAn8HzlHV+WW8/zhVLVDVgkaNGgUdvnPOuRB5YuUipWtXuw0bBgsWhB1N7vjrX20x5sJCm9NWlS1enH7/okXZjcPlDhFpALwIDFbV6WHH45xzLhyeWLnIuftuK5V98cW+tlUQ1q+3YZbt28Npp4UdTfhKK6nupdZdeUTkdyKyGDgKeFFEpiQeuhjYD7hORD5M3PYILVDnnHOh8MTKRU6zZjYP6MUX4W9/Czua+Lv3XuulGTHCFs6t6goLt62ImJdn+50ri6r+XVWbqmptVW2sqp0T+29S1bqq2ibl9m3Y8TrnnMsuT6xcJF1yiVVqu+SS9HNiXMWsWWMJ1YknwnHHhR1NNPTqZcUrmje3RLN5c9vu1SvsyJxzzjkXZ55YuUiqUcPWtlq2zIouuB1zxx2wcqUlV65Yr142h6+oyH56UuWcc865yvLEykVWu3Zw4YUwZgzMmhV2NPGzciXcfjt06wZHHBF2NM4555xzuc0TKxdpI0bAHnvY2lZbtoQdTbzcfDP89BPceGPYkTjnnHPO5T5PrFyk1a8Po0dbj9XYsWFHEx+LF1tPX+/etiCucy56Zs2CFi18cWrnnMsVnli5yDvzTOjcGa69FpYsCTuaeLjxRps/NGxY2JE458qycCH06+fJlXPO5QJPrFzkicA998CmTTBwYNjRRN/cufDQQ3DBBXY13DkXbevWwZAhYUfhnHOusjyxcrGw775WHXDiRFvfypXu+uuhdm0/UXMuThYutItHzjkXV088YRd0q1WrusOcPbFysXHFFTZfaMAAK8rgtvXhhzBhAlx6KTRuHHY0zrntcfDB9vtbVBR2JPHnJ3jOZdcTT9iw5oULQbXqDnMWVQ07howoKCjQmTNnhh2GC9ibb8Kxx8JVV8Ett4QdTfSccgpMnw5ffw0NGoQdjXMgIrNUtSDsOKJIpEBhJnl5NnT33/+GTz6BQw+FwkL47W9tKLTbPskTvHXrivfl5flC4C5ennjCRp4sWgT5+fY3IdPf340b4ccf7bZqVfH9iuz76qv0F4WaN7f1IqMsyHYq64mViBwAPJ2ya1/geqAB0BdYkdh/jar+K/GawcB5wBbgElWdUt5xPLHKXeedB48+Cu+/D4ccEnY00TF9OhxzDIwcCYMGhR2Nc8YTq9KJFGjz5jP/d8JUVGQ9VtdfD/Pnw9FH25ITHTqEHWm8tGhhV8tLisMJnnOwYxcHNm/eNunZ3gRp/fryY6tXzyo2N2hgP5O3CRPSP18k+r3wsU6stjq4SHVgCXAk8EdgraqOKvGcVsBTQDtgL+AVYH9VLXNVI0+sctd338GBB0LLlvDWWzbUo6pTtZOvuXPthCwvL+yInDOeWJWutHZq0yYYPx5uuAGWLoVOnexqdYF/imXavBkmTYLf/77051x0EbRta7dWraBmzezF51xFlXZxoF4968lOlyBVZIpEXt7WyVDJ5Ki8fbvsAtWrb1/MtWrBG2/AkUduzyeQXUG2UzWCeJNKOAGYr6oLpfTxDl2BCaq6AfhaROZhSdbbWYrRRUzDhjBqFPTpAw8+aFd1qropU2yY5NixnlQ5F3c1a9rftd69rSLqyJFwxBHQrZstpeBr021t+XJ44AG4/35bw6969fQLyteuDY89Zp9pcvvQQ4sTrbZtbZ5brVrZjd+5pEWL4KWX0icoAGvXwgcfFCc8e+1V8eSofv3MXkgoLNy2l61WLdhpJ2jf3s7Zbr459+d/h91j9TDwvqqOEZFhQB9gNTAT+D9V/UFExgDvqOrjidc8BExW1efSvF8/oB9Afn5+24WlfTNd7KnC8cdbsYbPP8/9X9SyFBXZlexVq+yz8JMCFyXeY1W6io6sWL0a7rwTbr/drkr37g1Dh8I++2QhyIhShbfftotJzz5rvXwnnggXX2yf1wUXpB9G1bMnzJtnizMnb++/b68B+/uZTLYOP9x+tm5tSZhzQduwwS6KvvQSTJ4Mn35q+0u7OBD14azp5oWddhrcdJP9DdtpJ+uJHzAgWr3FOTEUUERqAUuBg1V1uYg0BlYCCtwINFHVP4nIWODtEonVv1R1Ylnv70MBc9/nn1sDeOaZ8PjjYUcTnmeegbPOsnlnvXuHHY1zW/PEqnTb206tXGlFe8aMsZOufv3sJKZJkwwGGTHr1sGTT1pC9eGHNjSpTx8b4nfAAcXP256J/0VFNvG+ZLK1apU9XrOmzedN9modfri1PZ5suR3x1VeWRL30Erz6qn2na9WywlwnnQRdulivVK4VYPniC1uLdMoU63W/+267QB4FuZJYdQUGqGqnNI+1AF5Q1daJwhWo6sjEY1OAYapa5lBAT6yqhuuvt6ExL79sVyurms2bbehKzZowe3bpY5+dC4snVqXb0XZqyRL7u/fQQ/a7P3AgXHkl7LZbBoKMiHnz4N574eGHLeE55BC76t2rl807CZqqVVdNTbZmzYIffrDHa9SwnqzUYYSHHgp16gQfi4u3n3+GadOKe6XmzrX9++5bnEgddxzUrbv168KoCphpqjYP8rLL7Pere3fric/PDzeuXEmsJgBTVHV8YruJqi5L3L8MOFJVe4jIwcCTFBevmAq09OIVDqyCzSGHWNWZjz6qeo3aQw/B+efDP/4BXbuGHY1z2/LEqnSVbafmzYNhw6wHZ5ddLLkaODAziUYYtmyxE9GxY+2ktEYNm2d28cVWATXbpehVbRhWskcrmWx99509XqOGXehKDiFs2xZ++Usb/uSqDlX48sviROr11+1cpU4dS6CSyVTLlmFHGp6ff7a58iNH2vbgwfb3K6xzuNgnViKSB3wD7KuqPyb2PQa0wYYCLgD6pyRaQ4A/AZuBS1V1cnnH8MSq6njlFejY0Xqvbrgh7GiyZ/16+8O8994218DXu3FR5IlV6YJqpz7+GK67Dp5/Hho1sqvc/fvH90LTd99Zz9S999pV7SZN7N/Tt69N1o8SVetRSB1COGsWrEgsHFO9ug17Sh1G2KbNtkWGcrF3oipZuxZee614iN/XX9v+Aw4oTqSOPdaT7JIWLoQrroDnnrM5o3feaXOysn0+E/vEKhs8sapaevWyX8yPPtp6nH0uu/NOuPxymDo1OuOUnSvJE6vSBd1OvfMOXHONneA1a2a9WeecYz0pcTBzpvVOTZhgF446dLDhfqefHq2J7uVRteqEJYcRfvutPV6tGhx0UHGy9d13cNttdhU/Ke7zaXKdqhWaSCZSb75pi+vWrQsnnGCJVJcuVbvAzPaYOhUuucQ+086d4a67snsu54lVBXhiVbUsX25rW7VpY5NBc733Zs0aG5/dpo3NL3MuqjyxKl2m2qmpUy3BmjED9t/f5mN17x7NNf/Wr7cCPGPHWrx161oRngEDbA5TrlC1NclKJlv//W/pr2nUyE7a8/NtmZFcb9eibvVqGyHz0kt2++Yb29+6tSVRJ50Ev/qVFzXZUZs2WWGeYcPsIsOll1pP/M47Z/7YnlhVgCdWVc/991uJ3Ucesau0uWz4cCu3/O670K5d2NE4V7pcSqxE5AxgGHAQ0E5VZ5Z4PB/4FCuwNGrbd9haJtspVRsaeO21MGeOXYQpLLSTvyicoC9cCPfdZ2sRrlxpV6cHDLC/3fXrhx1d9ixdCk2b2v9XWXbayRKs/HwruZ28n7w1beon9EFTtaJQyV6p//zHCkbtsotNP+jSxXpXmjULO9Lcsnw5DBoEf/2rDQO+7Tb4wx8y+3fLE6sK8MSq6ikqsgnNc+daKfaGDcOOKDNWrrTeqhNPhL/9LexonCtbjiVWBwFFwP3AFWkSq4mJx98NO7FK2rIFnnrK5qB+/bX9jRwxAn7964weNq2iIrviP3YsvPCC7TvtNCtGcfzx0Uj4wtCiRfoFYffc0+aZLVpkt4ULi++X7OkSseeXTLhSE7Hddqu6n3FFff+9jQJJ9kolP+fDDivulWrfPl5DU+PqnXfgz3+2IcLHHGPl2du0ycyxgmynYjLy2rnyVatmV0APPxyuvtquhOaim2+2RUJvuinsSJyrWlT1MwBJc3YqIqcDXwE/ZTmsMlWvDmefbev9Pfyw9XYfe6ydJBYW2t/LTFu1ykYSjB1rF74aNbIr0v37h19mOQoKC9OvWTRqlM0vS2fDBpvH682MdwAAEKNJREFUlZpsJW+zZ8M//2nDLFPl5aVPuFJ7vSq6wHwci22ki7lnTxuSmazg9+67dgFg112hUydLpDp1qlprxUVF+/b2//Hww1Y1sG1b+5tx443RvnDuPVYu51x1lXUdv/FGOFdlM2nxYthvP+jRw7rJnYu6XOqxShKRaaT0WIlIXeAVoCNwBbC2tB4rEekH9APIz89vuzBdV0UGrVtnCc7NN9vV+e7d7UTlwAODP9ZHH9mxHn/cjnvUUTbcr3t3H7ZWUtCJiqqNbkhNuEomYcuXb/0aEUsgyur12nVXK+8ft8Vrn3hi25irV7chlmvX2r+9oKC4gl+7dr4uZJT88INNfxg7Fho0sN+Pvn2D+z/yoYAV4IlV1fXTT7aWSN26tnp5Ra/AxUH//jB+vK2R0aJF2NE4V764JVYi8gqwZ5qHhqjq84nnTGPrxGoUMENVnxGRYZSRWKUKs5368Ue44w67rVsH555rJy7Nm1fufTdutCHKY8fCW29Zyfc//MESqmz0jrmKW7/eCjCU7PFKTcA2bNj6NXXr2r7Nm7d9v112sf/noqKK31S37/k7+ppPPrHiCCXVrWvzszt1sp5UF20ff2zDA19/3YZn3n23FQypLE+sKsATq6rthRfg1FNtLsHgwWFHE4y5c61E74UX2h8T5+IgbolVRaRJrN4EklPYG2DzrK5X1TFlvU8U2qkVK6z3auxYO2Ht3996Tho33r73WbLEeizGjbN5KfvuCxddBH/8o83tcfGjat+PksnW6NGlv6ZmTRuWX62a9QIl72fqVtFj/POf6eMVscTLxYeqVRK94gobxXP22XDrrZUbrhloO6WqOXlr27atuqqtWzfVOnVU588PO5Jg9OihmpenumxZ2JE4V3HATI1AmxDkDZgGFJTy2DAs6Sq/napXT3X8ePugNm5U7dBB9bHHbPunn2x7wgTbXrXKtidOtO0VK2x70iTbXrbMtidPtu1Fi2z75Zdte/582542zbY//9y2p0+3l7/8sX7epIMeWW2G5uWp3n3eB7rxVx1UP/jAnj9jhj3/449te/p0LerQQd955HPt3l31uGrT9DU66J9+M19ffFF1y5SX7fmLFtnzJ0+27eQfsEmTbHvFCtueONG2V62y7QkTbPunn2z7scdse+NG2x4/3raTxo1TPeGE4u2xY1W7dCneHj1a9dRTi7dvu80aiaSRI1XPOqt4e/hw1V69irevu061T5/i7UGDVPv2Ld7+v/9Tveii4u2BA+2WdNFF9pykvn3tPZL69LFjJPXqZTEknXWWxZjUrZv9G5JOPdX+jUlduthnkHTCCfYZJXXosMPfvTZNV+hrdNBTmKSg2phl+hod9Jw9duy7px9/bNszZtj2Bx/YdhnfPe3Qwd5H1d63Q4fixv7lbb97b9fuoI1ZpqB6CpP0NTpoQ1Zo8+bq370YffdS/+6tXas64pJlOk066Ol1Juutt6pumLdj370g26kIrmrhXDDuussWxhwwoPxStlE3e7YtmjlwoFV+cs5ln4j8TkQWA0cBL4rIlLBjCsqee8IB+9uV4K5d4cGH4N134KGHbPL4aafBtNdtyNRDD9lwv5nvwTnn2tqBZ5wBRx5pj518cjTXzHLBGDJk2//fatVsOGlU7bMP5O209b68nWyujoununVtRFK7I+DQQ21+fadONnc0TD4U0OW0u+6yReZ2391Wt49L9aKSTjkFpk+Hr76yycPOxUUuDgUMSpTbqdmzbQ2sZFn0dA4/3Eql9+hhRQBc1ZErVQGjHrOruBdftPO9efOsmuYdd1hCXRE+x6oCotxguex57DG7ipb6NY969aKSpk+3NRxGjrQSxc7FiSdWpYtDO7XnnttWj0vuX7rU10VyzkXHhg2WUN10k82du+oqW34nL6/s1wXZTnlnvctp11237TDAdeusAMSoUTbs5Z13bOL1li3hxFgWVevqbtzYKuE451w2fftt+v3Ll3tS5ZyLltq17Zzpiy+s12r4cCv6NXFi9qaE+ALBLqctWpR+/5o1cOWVW++rUcMWSGzWrHjtjtT7+flQv37mY041ZQq8+SaMGWPjiZ1zLpvy860iXLr9zjkXRU2bwlNPwQUX2EXp7t3hhBPgL3+BVq0ye2xPrFxOK+2koHlzm0OwaNG263h8840Nv3v66W3X6th5560TrZJJ2PasXF+eoiK45hobI9y3bzDv6Zxz26OwMP1isD7p3zkXdR06wPvvw3332QimX/7SEq2hQzN3odwTK5fTyjopqF8fDjnEbuls2WLDXVITrtT7M2faGh+pRGzuQbqkK3m/UaOyh9CkTrBVtSsuubTIsXMuPpJzUX3Sv3MujmrUsCI7Z51lF6tHj4Ynn4RbboHq1a1ID7RtG9TxvHiFy3mZrAT088+WZJVMulJ7wH7+eevX1K5detI1e7ZdSUl9TdyKbTiXyotXlM7bKeecy66ZMy3RevddWybAFoguQHVmILNGPbFyLoNUbU2F0nq9Fi2yylrlrfzevDksWJCVkJ0LlCdWpfN2yjnnsq+oCPbYw5bhMcElVj4U0LkMEoGGDe122GHpn7NpEyxbZknWr3+d/jmlFeFwzjnnnHMVV61a5hYS9nLrzoWsZk0bBnjMMdYzlY5X4HLOOeecC0amzqs8sXIuQgoLt13IzitwOeecc84FJ935VhA8sXIuQnr1skIVzZvbMMLmzb1whXPOOedckFLPt4Lkc6yci5hevTyRcs4555zLpOT5lsisWUG9p/dYOeecc84551wlZT2xEpEDROTDlNtqEblURHYTkZdFZG7i564prxksIvNE5AsR6ZztmJ1zzjnnnHOuLFlPrFT1C1Vto6ptgLbAOuDvwCBgqqq2BKYmthGRVkAP4GCgC3CPiFTPdtzOOeecc845V5qwhwKeAMxX1YVAV+CRxP5HgNMT97sCE1R1g6p+DcwD2mU9Uuecc84555wrRdiJVQ/gqcT9xqq6DCDxc4/E/r2Bb1Jeszixbxsi0k9EZorIzBUrVmQoZOecc84555zbWmhVAUWkFnAaMLi8p6bZp+meqKrjgHGJ9/9ZROZUKsjsygcWhR3EdohbvBC/mOMWL8QvZo838w4OO4ComjVr1npvpzIubjF7vJkXt5jjFi/EL+bA2qkwy62fBLyvqssT28tFpImqLhORJsC3if2LgWYpr2sKLK3A+69V1YLgws0sEVnh8WZW3GKOW7wQv5g93swTER8+UDpvpzIsbjF7vJkXt5jjFi/EL+Yg26kwhwL2pHgYIMAk4NzE/XOB51P29xCR2iKyD9ASmFGB918VVKBZ4vFmXtxijlu8EL+YPd7Mi2PM2RK3zyZu8UL8YvZ4My9uMcctXohfzIHFG0qPlYjkAR2B/im7bwaeEZHzsO7DMwBUdY6IPAN8CmwGBqjqlgoc5sdgo844jzfz4hZz3OKF+MXs8WZeHGPOlrh9NnGLF+IXs8ebeXGLOW7xQvxiDizeUBIrVV0HNCyx7zusSmC65xcChdt5mHE7Fl1oPN7Mi1vMcYsX4hezx5t5cYw5W+L22cQtXohfzB5v5sUt5rjFC/GLObB4RTVtHQjnnHPOOeeccxUUdrl155xzzjnnnIs9T6ycc84555xzrpJik1iJSDMReU1EPhOROSIyMLF/NxF5WUTmJn7umvKawSIyT0S+EJHOKftfEpHZife5T0SqRznelMcnicgnQccadLwiMi2x78PEbY90x4xYzLVEZJyIfCkin4vI76Mar4jsnPLZfigiK0VkdNDxBhlzYn9PEflYRD5K/A7uHvF4z0rEOkdEbg061h2JV0QaJp6/VkTGlHivtonPd56I/EVE0q0BGLWYC0XkGxFZm4lYsy3g75+3UxmOWbLQVgUcr7dTGYw5sd/bqQBilpDbqoDj3b52SlVjcQOaAIcn7u8MfAm0Am4FBiX2DwJuSdxvBcwGagP7APOB6onHdkn8FGAi0CPK8SYe7wY8CXwSg893GlAQs+/EDcBNifvVgN2jHG+J950FHBvlzxgrlPNt8nNNvH5YhONtiFUnbZR43iPACRGIty5wDHABMKbEe80AjsL+rk0GTorId6KsmNsn3m9tJmLN9i3I33G8ncrGZzyNDLdVAcfr7VQGY8bbqSBjDrWtCjje7WqnAv/ws3XD1rnqCHwBNEn5IL9I3B8MDE55/hTgqBLvURP4J3BWlOMF6gFvJb4UGWuwAox3GllIrAKO+RugblziTdnXMhG7RDnmxO/aCqA59sf0PqBfhOM9AnglZX9v4J6w4015Xh9S/vgnnvN5ynZP4P4ofCdKi7nEYzmRWAX1/SvxHt5OZS7maWS5rapkvN5OZTBmvJ0KLOaU5/UhAm3VjsZb4rEKtVOxGQqYSkRaAIcB7wKNVXUZQOJnsit/b+wXOWlxYl/yPaZgVybWAM9FPN4bgduBdZmMMymIzxcYn+j+vy4T3bxBxiwiDRLbN4rI+yLyrIg0jmq8Jd6qJ/C0Jn7rM6kyMavqJuBC4GNgKXby9VBU4wXmAQeKSAsRqQGcDjSLQLyl2RuLPSnddyVwlYw5p3k7lXlxa6u8nfJ2Ksh4CaGd2o6YS5P1tirb7VTsEisRqYcNi7hUVVeX9dQ0+/73S62qnbFstTZwfKBBpgZRyXhFpA2wn6r+PSMBlgwimM+3l6oeAvw6cesdbJQlAql8zDWApsB0VT0ceBsYFXigySAC+g4n9ACeCiq2UgOp/Pe4JtZgHQbsBXyEXYXLiMrGq6o/YPE+DbwJLMAWKM+I7Yi31LdIsy+jJzEBxJyzvJ3KvLi1Vd5OeTu1TRAxa6cgfm1VGO1UrBKrxJd+IvCEqv4tsXu5iDRJPN4Eu7oHlgWnZu5NsSsQ/6Oq64FJQNcIx3sU0FZEFmDDLPYXkWkRjhdVXZL4uQYbb98uE/EGGPN32FXW5EnBs8DhEY43+V6/BGqo6qxMxBpwzG0AVHV+4qrlM8DREY4XVf2nqh6pqkdhwwfmRiDe0ixOxJ60zd+7IAUUc07ydiqz7VSAMWetrfJ2ytupDMWbtXZqB2IuTdbaqrDaqdgkVoku+oeAz1T1jpSHJgHnJu6fi42jTO7vISK1RWQfbIzvDBGpl/Kh1gBOBj6Paryqeq+q7qWqLbCJdV+q6m+iGq+I1JBEFZ3El/oUIFOVDIP6jBWbw/CbxPNOAD6Narwpr+tJhq8CBhjzEqCViDRKPK8j8FmE40USFcLEqgZdBDwYgXjTSgxpWCMi7RPveU55r9lRQcWci7ydymw7FWTM2WqrvJ3ydiqD8WalndrBmNPKVlsVajulWZjkFsQN+2OtWNfsh4nbyVhVlKlYlj4V2C3lNUOw6ilfkKg6AjQG3ku8zxzgbuxqSiTjLfGeLchcVcCgPt+6WPWf5Od7F2kqBEUp5sT+5sAbifeaCuRHOd7EY18BB8bh9y6x/wKskfoIO0FoGPF4n8JOXD4lAxXZKhHvAuB7YC129a9VYn8BdmI4HxhDhiaKBxzzrYntosTPYZn8Pmf6FtT3D2+nsvEZZ6WtCvIzxtupbHzG3k4FF/MCQmqrAo53u9opSbzIOeecc84559wOis1QQOecc84555yLKk+snHPOOeecc66SPLFyzjnnnHPOuUryxMo555xzzjnnKskTK+ecc84555yrJE+snHPOOeecc66SPLFyzjnnnHPOuUr6f98QjGP8q8acAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# rename columns to make them shorter and no whitespace\n",
    "daymetdf=daymetdf.rename(columns = {'prcp (mm/day)':'prcp', 'tmax (deg c)':'tmax', 'tmin (deg c)':'tmin'})\n",
    "daymetdf_30yrs=daymetdf_30yrs.rename(columns = {'prcp (mm/day)':'prcp', 'tmax (deg c)':'tmax', 'tmin (deg c)':'tmin'})\n",
    "\n",
    "# compute annual total precipitation\n",
    "daymetdf_ann = daymetdf.groupby(pd.Grouper(freq='1Y')).aggregate(np.sum)\n",
    "daymetdf_30yrs_ann = daymetdf_30yrs.groupby(pd.Grouper(freq='1Y')).aggregate(np.sum)\n",
    "\n",
    "# compute monthly mean temperature in winter season (Dec, Jan, and Feb)\n",
    "daymetdf_mon = daymetdf.groupby(pd.Grouper(freq='1M')).aggregate(np.mean)\n",
    "daymetdf_wtr = daymetdf_mon[(daymetdf_mon['yday']<59) | (daymetdf_mon['yday']>334)]\n",
    "daymetdf_wtr = daymetdf_wtr.groupby(pd.Grouper(freq='1Y')).aggregate(np.mean)\n",
    "\n",
    "daymetdf_30yrs_mon = daymetdf_30yrs.groupby(pd.Grouper(freq='1M')).aggregate(np.mean)\n",
    "daymetdf_30yrs_wtr = daymetdf_30yrs_mon[(daymetdf_30yrs_mon['yday']<59) | (daymetdf_30yrs_mon['yday']>334)]\n",
    "daymetdf_30yrs_wtr = daymetdf_30yrs_wtr.groupby(pd.Grouper(freq='1Y')).aggregate(np.mean)\n",
    "\n",
    "# make two plots\n",
    "# 1st plot shows time series of annual total precip in 2003-2011, comparing with 30-yr normal\n",
    "# 2nd plot shows time series of winter mean daily max/min temperature in 2003-2011, comparing with 30-yr winter normal\n",
    "fig = plt.figure(figsize=(14,4))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "\n",
    "daymetdf_ann.prcp.plot(legend=True, marker='o', style='b', ax=ax1)\n",
    "ax1.axhline(y=daymetdf_ann.prcp.mean(), color='r', linestyle=':')\n",
    "ax1.set(title='annual total precip in 2003-2011')\n",
    "ax1.set(ylabel='Precipitation (mm)')\n",
    "\n",
    "daymetdf_wtr.tmin.plot(legend=True, marker='o', style='b', ax=ax2)\n",
    "daymetdf_wtr.tmax.plot(legend=True, marker='o', style='r', ax=ax2)\n",
    "ax2.axhline(y=daymetdf_30yrs_wtr.tmin.mean(), color='r', linestyle=':')\n",
    "ax2.axhline(y=daymetdf_30yrs_wtr.tmax.mean(), color='b', linestyle=':')\n",
    "ax2.set(title='winter mean daily max/min temperature in 2003-2011')\n",
    "ax2.set(ylabel='Temperature (C)')\n",
    "\n",
    "print('Plots:')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "93Uh5TakE0sI"
   },
   "source": [
    "In the chart above, we can see that years 2004-2009 were dry (received less precipitation) compared with the 30-yr normal. During a drought period, the water-stressed trees are more susceptible to the bark beetles. Also, daily minimum temperature during winter time in years 2004-2009 was higher than normal (calculated from 1981-2010 Daymet data). Warmer winters favor the growth of bark beetle infestations.\n",
    "\n",
    "### 2. Download Vegetation Index with the MODIS Web Service\n",
    "\n",
    "Now, we will use ORNL DAAC's MODIS Web service to download the MODIS Vegetation Index (EVI from [MOD13Q1](https://doi.org/10.5067/modis/mod13q1.006)). The MODIS product MOD13Q1 is a gridded product with a time frequency of 16 days and a spatial resolution of 250 m.\n",
    "\n",
    "We will use  [modisViirsClient](https://github.com/tquaife/modisViirsClient) to access the MODIS Web service, and then print the bands available for product MOD13Q1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "xOO5booqR4Mo"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[\n",
      "  {\n",
      "    \"units\": \"bit field\", \n",
      "    \"band\": \"250m_16_days_VI_Quality\", \n",
      "    \"fill_value\": \"-1\", \n",
      "    \"description\": \"VI quality indicators\", \n",
      "    \"valid_range\": \"0 to 65534\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"Surface Reflectance Band 3\", \n",
      "    \"fill_value\": \"-1000\", \n",
      "    \"valid_range\": \"0 to 10000\", \n",
      "    \"band\": \"250m_16_days_blue_reflectance\", \n",
      "    \"units\": \"reflectance\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"units\": \"rank\", \n",
      "    \"band\": \"250m_16_days_pixel_reliability\", \n",
      "    \"fill_value\": \"-1\", \n",
      "    \"description\": \"Quality reliability of VI pixel\", \n",
      "    \"valid_range\": \"0 to 3\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"Surface Reflectance Band 2\", \n",
      "    \"fill_value\": \"-1000\", \n",
      "    \"valid_range\": \"0 to 10000\", \n",
      "    \"band\": \"250m_16_days_NIR_reflectance\", \n",
      "    \"units\": \"reflectance\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"Surface Reflectance Band 7\", \n",
      "    \"fill_value\": \"-1000\", \n",
      "    \"valid_range\": \"0 to 10000\", \n",
      "    \"band\": \"250m_16_days_MIR_reflectance\", \n",
      "    \"units\": \"reflectance\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"16 day NDVI average\", \n",
      "    \"fill_value\": \"-3000\", \n",
      "    \"valid_range\": \"-2000 to 10000\", \n",
      "    \"band\": \"250m_16_days_NDVI\", \n",
      "    \"units\": \"NDVI ratio - No units\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"Surface Reflectance Band 1\", \n",
      "    \"fill_value\": \"-1000\", \n",
      "    \"valid_range\": \"0 to 10000\", \n",
      "    \"band\": \"250m_16_days_red_reflectance\", \n",
      "    \"units\": \"reflectance\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"units\": \"Julian day of the year\", \n",
      "    \"band\": \"250m_16_days_composite_day_of_the_year\", \n",
      "    \"fill_value\": \"-1\", \n",
      "    \"description\": \"Day of year VI pixel\", \n",
      "    \"valid_range\": \"1 to 366\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.1\", \n",
      "    \"description\": \"Relative azimuth angle of VI pixel\", \n",
      "    \"fill_value\": \"-4000\", \n",
      "    \"valid_range\": \"-3600 to 3600\", \n",
      "    \"band\": \"250m_16_days_relative_azimuth_angle\", \n",
      "    \"units\": \"degrees\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.01\", \n",
      "    \"description\": \"Sun zenith angle of VI pixel\", \n",
      "    \"fill_value\": \"-10000\", \n",
      "    \"valid_range\": \"-9000 to 9000\", \n",
      "    \"band\": \"250m_16_days_sun_zenith_angle\", \n",
      "    \"units\": \"degrees\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.0001\", \n",
      "    \"description\": \"16 day EVI average\", \n",
      "    \"fill_value\": \"-3000\", \n",
      "    \"valid_range\": \"-2000 to 10000\", \n",
      "    \"band\": \"250m_16_days_EVI\", \n",
      "    \"units\": \"EVI ratio - No units\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }, \n",
      "  {\n",
      "    \"scale_factor\": \"0.01\", \n",
      "    \"description\": \"View zenith angle of VI Pixel\", \n",
      "    \"fill_value\": \"-10000\", \n",
      "    \"valid_range\": \"-9000 to 9000\", \n",
      "    \"band\": \"250m_16_days_view_zenith_angle\", \n",
      "    \"units\": \"degrees\", \n",
      "    \"add_offset\": \"0\"\n",
      "  }\n",
      "]\n"
     ]
    }
   ],
   "source": [
    "product = 'MOD13Q1' # MODIS Terra vegetation index\n",
    "\n",
    "response=modViirRequest(product)\n",
    "m=parseModViirJSON(response)\n",
    "print (json.dumps(m.bands, indent=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N1eQVdPGVQtB"
   },
   "source": [
    "The vegetation index we are interested is '250m_16_days_EVI' which contains metadata elements (e.g., scale_factor). Let's retrieve the `scale_factor` value and assign it to variable \"scale.\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Waxu2NmXMaOT"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scale factor: 0.0001\n"
     ]
    }
   ],
   "source": [
    "band = '250m_16_days_EVI' # EVI band\n",
    "\n",
    "evi_band = [item for item in m.bands if item['band'] == band]\n",
    "scale = float(evi_band[0]['scale_factor'])\n",
    "print(\"Scale factor: \" + str(scale))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "I-C1iTlcsgBc"
   },
   "source": [
    "Also note the '250m_16_days_VI_Quality' band which contains quality information about the band. For the sake of simplicity, we will not apply quality filters in this tutorial. Read more about the quality band [here](https://lpdaac.usgs.gov/documents/103/MOD13_User_Guide_V6.pdf). \n",
    "\n",
    "We will now download MODIS EVI data from the start of 2003 to the end of 2012. This might take a while to run, so sit back and relax while it finishes downloading the data. Once we download the data, we apply the scale factor and arrange data into a Pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "QqBVDhLmVRGg"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start: 2019-06-04 12:52:57.966000\n",
      "End: 2019-06-04 13:34:54.556000\n"
     ]
    }
   ],
   "source": [
    "#start and end times\n",
    "startdate=dt.datetime(year=2003,month=1,day=1)\n",
    "enddate=dt.datetime(year=2012,month=12,day=31)\n",
    "\n",
    "# count how long it takes to download all data\n",
    "print('Start: ' + str(dt.datetime.now()))\n",
    "\n",
    "# download and parse data\n",
    "response=modViirRequest(product,band=band,latitude=latitude,longitude=longitude,start_date=startdate, end_date=enddate)\n",
    "modEvi=parseModViirJSON(response)\n",
    "\n",
    "# apply scale factor and put data into a Pandas dataframe\n",
    "evidf = pd.DataFrame(data=scale*modEvi.data[band][:,0,0], index=modEvi.dates, columns=['evi'])\n",
    "\n",
    "print('End: ' + str(dt.datetime.now()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_vo8ym6hS4cq"
   },
   "source": [
    "Save the Pandas dataframe as a CSV file so that we don't have to download the data again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "soXACDw2Najy"
   },
   "outputs": [],
   "source": [
    "evifile = 'modis_evi_fairview_curve.csv'\n",
    "evidf.to_csv(evifile, header=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "TR59EDlBzIFb"
   },
   "source": [
    "Print the head of data frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "kFt743oizdDg"
   },
   "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>evi</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2003-01-01</th>\n",
       "      <td>0.1654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003-01-17</th>\n",
       "      <td>0.1694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003-02-02</th>\n",
       "      <td>0.2017</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003-02-18</th>\n",
       "      <td>0.2767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003-03-06</th>\n",
       "      <td>0.1556</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               evi\n",
       "2003-01-01  0.1654\n",
       "2003-01-17  0.1694\n",
       "2003-02-02  0.2017\n",
       "2003-02-18  0.2767\n",
       "2003-03-06  0.1556"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evidf[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9gI7fGSOTRrG"
   },
   "source": [
    "Let's compute and plot the annual mean EVI values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Xf_ohEO0TFjT"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plot: \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEICAYAAACXo2mmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8VfX9x/HXJwMSIBDCJgkkYYpsQxQRtQ4QqwUHOFBcCNhqrbPS9tdqa6t1tVongqJ1IlLcxV32CFtGmAGSsCEQIAkZn98f9wQvMSEJ3Jtz783n+Xjch/eecc/nfg33fb/ne4aoKsYYY4wvhLldgDHGmNBhoWKMMcZnLFSMMcb4jIWKMcYYn7FQMcYY4zMWKsYYY3zGQsUELBE5X0Sy3K4jUIhIfRFZLSKt3a4lUIhITxGZ63Yd5kcWKuYYEfleRPaLSH23a6kOEVER2SkiEV7TIkRkl4jU+glYInKziJSIyKFyj7YiMkNE/lzBOkNFZIdT92QRefQEmxgDzFTVHc66D4jIDyKSJyKbReSBcu+dJCLficgREVkrIheVm3+9iGwRkcMiMl1E4rzmPSEi20TkoLPM76v47E+JyHqnlrUiMqrc/N4istipZbGI9Paa191pnz3l/785QTrJqSFPRJaKyJCy+aq6AsgVkctPVJ+pPRYqBvB8AQEDAQV+4WoxNZMLDPF6fSmw36VaAOapaqNyjxxgMnCjiEi55W8E3lbV4mq891jg316vBRgFNAUuAe4UkWu95r8LLAWaAb8HpopICwAROR14xdl+K+AI8KLXupOArqraGDgbuF5ErjxBbYeBy4EmwE3AsyJytrOtesBHwFtOrW8AHznTAYqAKcBtFbxvBLANOM957/8Dpjh/r2XedtrGBAJVtYc9AP4IzAGeAT4tN28y8ALwGZAHLAA6eM1XYBywHs8X+guAOPMeBt7yWjbJWT7CeX0LsMZ5303AWK9lzweyTlCzAn8APvCaNhXPF6h6TWuC50tyO5ANPAqEO/M6AN8Ce4E9eL6gYr3WzQTuB1YAB4D3gahK6rkZmF3JvGhn/XO9pjUFCoBeXu38aCXrtwPyy9qtkmWeA/7lPO8MFAIxXvNnAeOc538D3vGa1wE46r2817x4YCXwYA3+nj4G7nOeD3LaXbzmbwUuKbdOR+//byd47xXAVeXqywfqu/3vyB5qPRVzzCg8X6hvA4NFpFW5+dcBj+D5ItwA/LXc/MuAfkAvYAQwuJrb3eWs2xhPwPxDRPrWoO7pwLkiEisisXh6Wx+VW+YNoBjPl1YfPF9yo515AjwGtAVOAxLxBKG3EXh6AslATzzhUSOqmo/n17j3bqERwFpVXV6Nt+gBbNJKejROD2ggsMqZdLqzfJ7XYsud6WXzj21XVTfiCZXOXu/5kIgcArKAhsA71agTEYnG87fgXcsKdRLAscKrlmpz/i47e703qpqNp7fTpabvZ3zPQsUgIucA7YEpqroY2AhcX26xaaq60PlSexvoXW7+46qaq6pbge8qmF8hVf1MVTeqx/+AL/F8OVZXAfAJcA1wLZ5fyAVen60Vnt1jv1HVw6q6C/iHsyyqukFVv1LVQlXdjaendl65bTynqjmqus/Z1ok+21kikuv12Og17w1guPOlC56AeaOanzMWT2+uMg/j+ff8uvO6EZ6ekbcDQEw156Oqjzuv++LZ7VZ++cq8jCewZlR3W9UhIpF4/vbeUNW15Wbn4Wkj4zILFQOefeBfquoe5/U7zjRvO7yeH8HzRVGT+RUSkSEiMl9E9olILp4xkebVrtzjTTxf0KOc597aA5HA9rIvejxjCS2d7bcUkfdEJFtEDuLZ719++zX5bPNVNdbr0aFshqrOBnYDQ0UkBc+v+Wr9+sezW7HCL2ERuRPPZ/+5qhY6kw/h6f15a8yPwVTV/LKaVVWX4tm99IizvZe9DkL4XblangS6AyO8eibV2taJiEgYnmA7CtxZwSIxeMbXjMsiql7EhDLnV/MIIFxEyr486wOxItKrmrtmTuQw0MDr9bHDYZ2jzD7E84X4kaoWich0PLukamIW0AbPGMtsPOMDZbbhGVtoXsmuo8ec9Xqq6l4RGQY8X8Pt10RZAHbBE+Q7q7neCiBFRCK8P4eI3Ao8hGesxvvw61XO8jFeu8B68WOIrXJel71PCp7/7+sq2X4ETruq6jg8Y2jHEZFH8PQKz1PVg+VquU9ExCtoeuIZe6uSs2tvEp4DCi5V1aJy89sC9YCM6ryf8S/rqZhhQAnQDc9und54xhZmcfz+/5O1DM+YRzsRaQKM95pXD88X2W6g2DlUdFBNN+B8UV0O/KLcfntUdTueXWpPi0hjEQkTkQ4iUraLKwbPL+lcEYkHjjss1w/eBC4Cbqf6u75wAmM9kFY2TURG4hlwv1hVN5Vbfh2etv+TiESJyBV4vsg/dBZ5G7hcRAaKSEPgz3h2ceY5bTRWRJqKRxrwK+CbyuoTkfF4dplerKp7y83+Hs/f2K+dQ4TLehrfOuuKiETh+XvAqdf7sPaX8PxNXu6MTZV3PvCtVy/NuMhCxdwEvK6qW1V1R9kDz6/1keJ1DsjJUNWv8BwxtQJYDHzqNS8P+DWeAez9eL6UPj7J7axS1VWVzB6F5wtrtbOdqXh6NuDZpdMXzz7+z4BpJ7N9L/3lp+ep9POqMxOYi2fgu6aftewQ4DKP4jlceJHXtl72mn8tkIrnMz8OXO2MG+G01Tg84bILT7j+0mvdK/CMreXh2SX4L+dRmb/hOUJtffldY6p6FM+Pl1F4dlHdCgxzpoNnF2U+Pw6+5+P0OkSkPZ7DhXsDO7zee6TXtkfiGccxAUDK/bAzxgQo59f7UuBCpwdW54lID2CCqvZ3uxbjYaFijDHGZ2z3lzHGGJ+xUDHGGOMzFirGGGN8JqTOU2nevLkmJSW5XYYxxgSNxYsX71HVFr56v5AKlaSkJNLT090uwxhjgoaIbPHl+9nuL2OMMT5joWKMMcZnLFSMMcb4TEiNqRhjjD8UFRWRlZVFQUFB1QsHqKioKBISEoiMjPTrdixUjDGmCllZWcTExJCUlMRP7wgd+FSVvXv3kpWVRXJysl+3ZaFiTDnTl2bz5IwMcnLzaRsbzQODuzCsT7zbZRkXFRQUBG2gAIgIzZo1Y/fu3X7floWKMV6mL81m/LSV5BeVAJCdm8/4aSsBLFjquGANlDK1Vb8N1Bvj5ckZGccCpUx+UQlPzrD7PxlTHRYqxnjJya3oHlCVTzcmEJ199tmubdt2fxnj2JVXQP2IMAqKS38yr21stAsVmWDl9rjc3Llza21b5fm1pyIil4hIhohsEJGHKpg/UkRWOI+5IuJ9z+xYEZkqImtFZI2I2E14jN98sXI7g/8xk6KSUiLDj9/3XD8ijAcGd3GpMhNsysblsnPzUX4cl5u+NPuU3vett94iLS2N3r17M3bsWF544QUefPDBY/MnT57MXXfdBUCjRo1OaVunwm89FREJB14ALgay8Nzy9GNVXe212GbgPFXd79yffAJwpjPvWeC/qnq1iNQDGvirVlN3Hcgv4uGPV/Gfpdn0iG/CM2N7sSrn4LFfmQBdWzeyQXpzzCOfrGJ1zsFK5y/dmsvRkuN7u/lFJTw4dQXvLtxa4Trd2jbmT5efXul7rlmzhvfff585c+YQGRnJL3/5Sxo1asS0adN44oknAHj//ff5/e9/fxKfyLf8ufsrDdigqpsAROQ9YCie+4QDoKrefbT5QIKzbGPgXOBmZ7mjwFGM8aFZ63fz4NQV7Mor5O4LO3HnBR2JDA+jU6uYYyHy+BdrmTBzI1v3HqFdM/tdY6pWPlCqml4d33zzDYsXL6Zfv34A5Ofn07JlS1JSUpg/fz6dOnUiIyODAQMGnPQ2fMWfoRIPbPN6ncWPvZCK3AZ84TxPAXYDrzu7xBYDd6vq4fIricgYYAxAu3btfFC2CXX5R0t4/Is1vDFvCx1aNGTaHWfTKzG2wmVvPjuJibM28dqczTz8i8p/SZq640Q9CoABj39LdgUHdsTHRvP+2JPbi6+q3HTTTTz22GPHTZ80aRJTpkyha9euXHHFFQFx2LM/x1Qq+nRa4YIiP8MTKr91JkUAfYGXVLUPcBj4yZgMgKpOUNVUVU1t0cJntwQwIWrp1v38/LlZvDFvC7cOSOazXw+sNFAAWjeJ4he92zIlfRu5R6yzbKr2wOAuREeGHzctOjL8lMblLrzwQqZOncquXbsA2LdvH1u2bOHKK69k+vTpvPvuu1xzzTWnVLev+DNUsoBEr9cJQE75hUSkJzARGKqqe73WzVLVBc7rqXhCxpiTcrS4lKdmZHDVS3MpLC7lndFn8sfLuxFV7h9/RUafk8KRoyW8vaDi/eHGeBvWJ57HruxBfGw0gqeH8tiVPU5pXK5bt248+uijDBo0iJ49e3LxxRezfft2mjZtSrdu3diyZQtpaWm++xCnwJ+7vxYBnUQkGcgGrgWu915ARNoB04AbVXVd2XRV3SEi20Ski6pmABfiNRZjTE1k7MjjnveXsXr7Qa4+I4E/Xt6NxlHVv6het7aNGdipOW/MzWT0wGTqR1QdRKZuG9Yn3ucHd1xzzTUV9kY+/fTTn0w7dOiQT7ddE37rqahqMXAnMANYA0xR1VUiMk5ExjmL/RFoBrwoIstExPu2jXcBb4vICqA38Dd/1WpCU0mpMmHmRi7/12x2HizglRvP4KnhvWoUKGVGD0xhV14hHy/7SWfbGOPFryc/qurnwOflpr3s9Xw0MLqSdZcBqf6sz4SurXuPcP8Hy1mYuY9B3Vrxtyt70LxR/ZN+v3M7NadLqxgmztrM1WckBMSAqDGByC7TYkKKqvLewq0MeXYma7Yf5OnhvXjlxjNOKVDAczG+0QOTydiZx8z1e3xUrQkmqhUeZxQ0aqt+CxUTMnYdLOC2N9J5aNpKeiXG8t97zuUqH/YqftG7LS1j6jNx1iafvJ8JHlFRUezduzdog6XsfipRUVF+35Zd+8uEhM9WbOf301eSf7SEP13ejZv6JxEW5ttdVPUjwrnp7CSenJHB6pyDdGvb2KfvbwJXQkICWVlZtXI/En8pu/Ojv1momKB24EgRf/z4Bz5alkOvhCY8PaI3HVv677pHI89sxwvfbWDi7E08M6K337ZjAktkZKTf75gYKmz3lwlaM9ftZvA/Z/LZiu3cc1FnPrzjbL8GCkBsg3qMSE3k42U57DgQvPcrN8ZfLFRM0DlytJg/TF/JqNcW0igqgv/8cgB3X9SJiPDa+XO+dUAypapMnptZK9szJphYqJigsnjLfi59dhZvL9jK6HOS+fSuc+iR0KRWa2jXrAGXdG/N2wu2cKiwuFa3bUygs1AxQeFocSlP/Hctw1+eS1GJ8s7os/jDZdW7zIo/3D4whbyCYqYs2lb1wsbUITZQbwLemu0HuXfKctZsP8iI1AT+77JuxJzEWfG+1KddU1LbN2XS7M2M6t++1na9GRPo7F+CCVglpcrL/9vI0OfnsDuvgImjUnni6l6uB0qZ289NITs3n/+u2uF2KcYEDOupmIC0Ze9h7puynPQt+7nk9Nb89YruNDvFs+J97aLTWpHUrAGvztzEz3u0sUu3GIP1VEyAUVXeXrCFIc/OImNnHv+4phcv3dA34AIFIDxMuG1gCsuzDrAoc7/b5RgTEKynYlw1fWn2sfvBt2ocRWyDSNbuyOOcjs154uqetI2NdrvEE7q6bwLPfJnBhJmbSEuOc7scY1xnPRXjmulLsxk/bSXZufkosONgAWt35HFV33jevDUt4AMFILpeODee1Z5v1u5k42737mFhTKCwUDGuyMnN55FPVpFfVPKTefM37fP5dbv86cb+SUSGhzFp9ma3SzHGdbb7y/hdaamyftchFmXuIz1zH4sy95Odm1/p8jknmBeIWsTU58o+8Xy4OIv7Lu4ckOM/xtQWCxXjc4XFJax0Bq/TM/eRvmU/B/KLAM8XcFpSHKMHJvPi9xvZnVf4k/WDYbdXeaMHJvPeom28NX8rd1/Uye1yjHGNhYo5ZQfyi1iyZb/TE9nPsqxcjhaXAtChRUOGdG9NalIc/ZKa0i6uwbFDb5s2qMf4aSuP2wUWHRnOA4O7uPI5TkXHljFc0LUlb87LZOx5Ka6d6W+M2yxUTI3l5OYfC5BFmfvI2JmHKkSECd3jm3BT//akJsWR2r7pCXcFDesTD3Ds6K+2sdE8MLjLsenBZvTAZK5/dQH/WZrNdWnt3C7HGFdYqJgTOtF4SMN64fRt35RLe7QhNakpvRNjaVCvZn9Sw/rEB22IlNc/pRnd4xvz6qxNXJOaGFQHGxjjKxYq5jjVHQ/plxRH19Yxds0rLyLC7QNTuPu9ZXyXsYsLT2vldknG1DoLlTrE+0TDsl1NP+vaslrjIWlJcSTGRdulSKpwaY82/P2LtUyYuclCxdRJFip1RNmJhmWD4tm5+dzz/jLUmV/T8RBTscjwMG4ZkMxfP1/DiqxceibEul2SMbXKQiWEqSobdh1iYeY+Hv10zU9ONFQgJiqCCTem0jsxluh6dsSSL1yblshz36zn1Vmb+dd1fdwux5haZaESQopKSvkh+wCLnAH19Mx97D9SdMJ1DhUU079Ds1qqsG6IiYrk2rREXpuTyW8v6UJC0wZul2RMrbFQCWKHC4tZujWXhc6RWUu35h7rjSQ3b8hFp7WiX7JnPGTkxPlk5xb85D2C8UTDYHDLgGRen5PJ63My+b/LurldjjG1xkIliOw9VMgi59yQRZn7WJVzkJJSJUygW9vGXNMvkbTkOFKTmtIyJuq4dR8Y3DVkTjQMBm1jo/l5zza8t3Arv76wE02iA+PGYsb4m4VKgFJVsvbns3CzJ0AWZu5j0+7DANSPCKN3Yix3nNeBfslx9G0XW+XdEEPtRMNgcPvAFD5alsN7C7cy9rwObpdjTK2wUAkQJaVKxo480rfsOxYkOw96rovVOCqCfklxjEhNpF9SU7rHN6F+RM0H1UPpRMNg0D2+Cf1TmvH6nExuGZBMvQg7p8eEPr+GiohcAjwLhAMTVfXxcvNHAr91Xh4C7lDV5V7zw4F0IFtVL/Nnrb5W0Tkh3l/ohcUlrMhyBtU3e04yzCsoBqBNkyjOTG5Gv2TP9bI6t4yxs7OD1JhzU7hl8iI+W5nDFX0S3C7HGL/zW6g4gfACcDGQBSwSkY9VdbXXYpuB81R1v4gMASYAZ3rNvxtYAzT2V53+UNE5IQ9NW8Hq7QeIDA9j0ebjTzLs2LIRl/VsS1pyU/olxREfaycZhorzOregY8tGvDpzM8N6x9v/VxPy/NlTSQM2qOomABF5DxgKHAsVVZ3rtfx84NhPORFJAH4O/BW41491+tyTMzJ+ck5IQVEpE2ZuJiJMON05ybBfUhypSXHENaznUqXG38LChNsHJvPbD1cyd+NeBnRs7nZJxviVP0MlHtjm9TqL43sh5d0GfOH1+p/Ag0DMiTYiImOAMQDt2gXGlWEru8mUACseHlTjiy6a4Da0dzxPzsjg1VmbLFRMyPPnyGFF/XytYBoi8jM8ofJb5/VlwC5VXVzVRlR1gqqmqmpqixYtTqVen1BVGtWvODTaxkZboNRBUZHh3NQ/ie8zdrNuZ57b5RjjV/4MlSwg0et1ApBTfiER6QlMBIaq6l5n8gDgFyKSCbwHXCAib/mxVp9QVf7y6RryCosJLzewbueE1G03nNWeqMgwJs7a5HYpxviVP0NlEdBJRJJFpB5wLfCx9wIi0g6YBtyoquvKpqvqeFVNUNUkZ71vVfUGP9Z6ykpKlfHTVvLanM3cMiCJp67q6RlwB+Jjo3nsyh52OG8d1rRhPYafkcj0pTnsyvvplQ2MCRV+2xejqsUicicwA88hxa+p6ioRGefMfxn4I9AMeNE5KqZYVVP9VZO/FJWUcu+U5XyyPIe7LujIvRd3RkS44gw7hNT86LZzknlrwRbenLuF+63XakKUqFY4zBGUUlNTNT09vVa3WVBUwp3vLOHrNbt4aEhXxtmZ0+YExv47nfmb9jFv/AU2vmYCgogs9uWPeTvF9xQcLizmtjcW8fWaXfxl6OkWKKZKtw9M4UB+EVMXZ7ldijF+YaFykg7kFzHqtYXM27iXp4f34sb+SW6XZILAGe2b0qddLBNnbaakNHT2EhhTxkLlJOw9VMj1r85nRVYuL1zfl6ts7MRUU9l97LfuO8JXq3e4XY4xPmehUkM7DhRwzYT5bNh1iFdHpTKkRxu3SzJBZvDprUmMi2bCTDu82IQeC5Ua2LbvCMNfmcv23HzeuDWN87u0dLskE4TCw4TbBiSzZGsui7fsc7scY3zKQqWaNuw6xPCX53Ewv5i3bz+Ls1LsFrzm5A1PTaRJdCSvztzsdinG+JSFSjWsyjnANa/Mo7hUeX/sWfROjHW7JBPkGtaPYOSZ7Zixegdb9h52uxxjfMZCpQqLt+znugnzqR8RxpSxZ9G1dVBdhd8EsJvPTiIiTJg023orJnRYqJzA3A17uHHSAuIa1mPKuP6ktGjkdkkmhLRsHMXQ3vF8kJ7F/sNH3S7HGJ+wUKnEN2t2cvPkRSQ0jWbK2P4kNG3gdkkmBN0+MIX8ohLeXrDF7VKM8QkLlQp8uiKHsf9eTNfWMbw/pj8tG0e5XZIJUV1ax3Bu5xZMnruFwuKSqlcwJsBZqJQzZdE2fv3uUvq2a8rbo8+kqd2V0fjZmIEp7DlUyEdLf3JnCGOCjoWKl9fnbObBD1cwoGNz3rg1jZioSLdLMnXAgI7N6No6hldnbSKULvBq6iYLFccL323gkU9WM/j0Vky8KZXoeuFul2TqCBFhzLkprN91iO/X7Xa7HGNOSZ0PFVXl7/9dy5MzMriiTzwvXN+X+hEWKKZ2XdazLa0a17c7Q5qgV6dDpbRU+dPHq3jp+42MPLMdTw/vRUR4nW4S45J6EWHcMiCZORv2sirngNvlGHPS6uw3aHFJKQ9MXcGb87Yw5twUHh3WnbBy95U3pjZdl9aOhvXCmTjLToY0watOhsrR4lLuencpHy7J4p6LOjN+SFec2xkb45om0ZFc068dnyzPYfuBfLfLMeak1LlQKSgqYcy/0/nihx384eencfdFnSxQTMC4ZUASpapMnpPpdinGnJQ6FSqHCou56bWF/G/dbh67sgejB6a4XZIxx0mMa8ClPdrwzoKt5BUUuV2OMTVWZ0Il98hRRk5cQPqW/fzzmt5cl9bO7ZKMqdDtA1PIKyzm/UXb3C7FmBqrE6GyO6+QayfMZ03OQV4a2ZehvePdLsmYSvVKjCUtOY7X52RSXFLqdjnG1EjIh0pObj4jXpnHlr1HmHRzKoNOb+12ScZU6faBKWTn5vP5D3YfexNcQjpUMvccZvjL89iTV8i/b0tjYKcWbpdkTLVc2LUlKc0b8upMu3SLCS4hGyoZO/IY/so8jhwt5t0xZ5GaFOd2ScZUW1iYcNvAZFZmH2DBZruPvQkeIRkqK7JyuWbCPASYMrY/3eObuF2SMTV2Vd8E4hrW49WZdukWEzxCLlQWZe7j+lcX0LBeBB+M60+nVjFul2TMSYmKDOfGs9rzzdpdbNh1yO1yjKmWkAqVldkHGP7yPKLrhTH1jv60b9bQ7ZKMOSU39m9P/YgwJs223ooJDn4NFRG5REQyRGSDiDxUwfyRIrLCecwVkV7O9EQR+U5E1ojIKhG5uybbzSsoZsEm2w9tgl/zRvW5sm8CHy7JZs+hQrfLMaZKfgsVEQkHXgCGAN2A60SkW7nFNgPnqWpP4C/ABGd6MXCfqp4GnAX8qoJ1K1VQVMqTMzJO9SMYExBGD0zmaHEpb86z+9ibwOfPnkoasEFVN6nqUeA9YKj3Aqo6V1X3Oy/nAwnO9O2qusR5ngesAWp0xmJOrl2Qz4SGDi0acdFpLXlr/hbyj9p97E1g82eoxAPe15nI4sTBcBvwRfmJIpIE9AEWVLSSiIwRkXQRSfee3jY2uoblGhO4bh+Ywr7DRznrsW9IfugzBjz+LdOXZrtdljE/4c9QqejSvxWexSUiP8MTKr8tN70R8CHwG1U9WNG6qjpBVVNVNbVsWnRkOA8M7nLShRsTaHJy8xGBA/lFKJCdm8/4aSstWEzA8WeoZAGJXq8TgJzyC4lIT2AiMFRV93pNj8QTKG+r6rTqbjQ+NprHruzBsD52fS8TOp76ch3lT6zPLyqxsUMTcCL8+N6LgE4ikgxkA9cC13svICLtgGnAjaq6zmu6AJOANar6THU32CO+CXMeusAXtRsTUCobI7SxQxNo/NZTUdVi4E5gBp6B9imqukpExonIOGexPwLNgBdFZJnXuMgA4EbgAmf6MhG51F+1GhPoKhsjtLFDE2j82VNBVT8HPi837WWv56OB0RWsN5uKx2SMqZMeGNyF8dNWkl/049FfUZFhNnZoAo5fQ8UY4xtlY4RPzsgg29nlNbR3Wxs7NAHHQsWYIDGsTzzD+sSjqgx7cS6z1++lsLiE+hHhbpdmzDEhde0vY+oCEeH+QZ3Jzs3nvYV2y2ETWCxUjAlC53RsTlpyHM9/t8HOsjcBxULFmCDk6a10YXdeIf+en+l2OcYcY6FiTJBKS47j3M4teOn7jeQVFLldjjGAhYoxQe2+izuz/0gRr8/JdLsUYwALFWOCWq/EWAZ1a8WrMzeRe+So2+UYc+JQEZG4Ez1qq0hjTOXuHdSZQ0eLmWD3sjcBoKrzVBbjubJwZVccTvF5RcaYGunaujGX92zL63MyuWVAMi1i6rtdkqnDTthTUdVkVU1x/lv+YYFiTID4zUWdOFpSykvfb3S7FFPHVbX7a7WI/E5ELECMCWApLRpxVd943lqwhe0H7MrFxj1VDdRfB8QAX4nIAhH5jYi0rYW6jDE1dNcFnVBV/vXtBrdLMXVYVbu/lqvqeFXtANwNtAfmi8i3InJ7rVRojKmWxLgGXJfWjimLtrF17xG3yzF1VLUPKVbV+ap6DzAKaAo877eqjDEn5Vc/60h4mPDPb9ZVvbAxflCtUBGRfiLyjIhsAR4BJgB2zW1jAkyrxlGM6t+e6Uuz2bArz+1yTB1U1UD930RkI/ASnvvLD1DV81T1JVXdUysVGmNqZNx5HYiODOfvNRF6AAAWZ0lEQVQfX613uxRTB1XVUykEhqhqqqo+papZtVGUMebkNWtUn1vPSeazldtZlXPA7XJMHVNVqOSr6joAERnuPUNE/ua3qowxp2T0wBQaR0XwzJc2tmJqV1Whcq3X8/Hl5l3i41qMMT7SJDqSsed14Ju1u1iydb/b5Zg6pKpQkUqeV/TaGBNAbj47iWYN6/H0lxlul2LqkKpCRSt5XtFrY0wAaVg/gjvO78CcDXuZu9GOqzG1o6pQ6SUiB0UkD+jpPC973aMW6jPGnIIbzmpP68ZRPP3lOlTtd6Dxv6rOqA9X1caqGqOqEc7zsteRtVWkMebkREWGc+cFHVm8ZT/fr9vtdjmmDrCbdBkT4kakJpIYF83TX2ZYb8X4nYWKMSGuXkQYd1/YmR+yDzJj1Q63yzEhzkLFmDpgWO+2pLRoyDNfraOk1Horxn8sVIypAyLCw7j34s6s23mIT5bnuF2OCWEWKsbUEZd2b8NpbRrzz6/XUVRS6nY5JkT5NVRE5BIRyRCRDSLyUAXzR4rICucxV0R6VXddY0zNhIUJ913cmcy9R/hwsV3Gz/iH30JFRMKBF4AhQDfgOhHpVm6xzcB5qtoT+AueS+pXd11jTA1deFpLeiXG8tw36yksLnG7HBOC/NlTSQM2qOomVT0KvAcM9V5AVeeqatmFieYDCdVd1xhTcyLCA4O6kHOggHcXbHW7HBOC/Bkq8cA2r9dZnPjGXrcBX9R0XREZIyLpIpK+e7ed3GVMVQZ0bMaZyXE8/91G8o9ab8X4lj9DpaILTlZ4LKOI/AxPqPy2puuq6gTnfi+pLVq0OKlCjalLRIT7B3dhz6FC3piX6XY5JsT4M1SygESv1wl47h55HBHpCUwEhqrq3pqsa4w5Of2S4jivcwte/t9G8gqK3C7HhBB/hsoioJOIJItIPTz3ZvnYewERaQdMA24suxlYddc1xpya+wZ1JvdIEa/NznS7FBNC/BYqqloM3AnMANYAU1R1lYiME5FxzmJ/BJoBL4rIMhFJP9G6/qrVmLqoZ0Isg09vxcRZm8g9ctTtckyIkFC6wFxqaqqmp6e7XYYxQSNjRx6XPDuTced14LeXdHW7HOMCEVmsqqm+ej87o96YOqxL6xh+0astk+dksiuvwO1yTAiwUDGmjvvNRZ05WlLKS99vdLsUEwIsVIyp45KbN+Tqvgm8PX8rObn5bpdjgpyFijGGuy7siKL869sNbpdigpyFijGGhKYNuD6tHR+kb2PL3sNul2OCmIWKMQaAX/2sIxHhwrNfr3e7FBPELFSMMQC0bBzFTf2T+M+ybNbvzHO7HBOkLFSMMceMPa8DDSLD+cfX66pe2JgKWKgYY46Ja1iP285J5vOVO/gh+4Db5ZggZKFijDnObQNTaBIdyTNfWW/F1JyFijHmOE2iIxlzbgrfrt3F4i37q17BGC8WKsaYn7hlQBLNG9Xj6S8z3C7FBBkLFWPMTzSoF8Ed53dk7sa9zN2wx+1yTBCxUDHGVGjkme1o3TiKp77MIJSuZm78y0LFGFOhqMhw7rqwI0u25vJ9xm63yzFBwkLFGFOpEamJtItrwFNfZlBaar0VUzULFWNMpSLDw7j7wk6syjnIjFU73C7HBAELFWPMCQ3rE0+HFg15+qt1lFhvxVTBQsUYc0LhYcK9F3dhw65DfLw82+1yTICzUDHGVGlI99ac1qYx//x6PUUlpW6XYwKYhYoxpkphYcL9gzqzZe8Rpi7OcrscE8AsVIwx1XJB15b0TozluW/WU1BU4nY5JkBZqBhjqkVEeGBwF7YfKODdhVvdLscEKAsVY0y1nd2hGWelxPHCdxs5crTY7XJMALJQMcZUm4hw/6Au7DlUyBtzt7hdjglAFirGmBpJTYrj/C4tePl/GzlYUOR2OSbAWKgYY2rsvou7cCC/iEmzNrtdigkwFirGmBrrkdCES05vzaTZm9l/+Kjb5ZgA4tdQEZFLRCRDRDaIyEMVzO8qIvNEpFBE7i837x4RWSUiP4jIuyIS5c9ajTE1c++gzhwqLGbgE9+R/NBnDHj8W6YvtTPu6zq/hYqIhAMvAEOAbsB1ItKt3GL7gF8DT5VbN96Znqqq3YFw4Fp/1WqMqbnVOQcJF+FQYTEKZOfmM37aSguWOs6fPZU0YIOqblLVo8B7wFDvBVR1l6ouAioa7YsAokUkAmgA5PixVmNMDT05I4OScjfvyi8q4ckZdgviusyfoRIPbPN6neVMq5KqZuPpvWwFtgMHVPXLipYVkTEiki4i6bt3242EjKktObn5NZpu6gZ/hopUMK1a180WkaZ4ejXJQFugoYjcUNGyqjpBVVNVNbVFixYnXawxpmbaxkZXOL1xdITdfjgITF+azYDHv6Ve645n+PJ9/RkqWUCi1+sEqr8L6yJgs6ruVtUiYBpwto/rM8acggcGdyE6Mvy4aWECB/KLufOdpRwqtDPuA9X0pdmMn7aSbD/0KiN8/o4/WgR0EpFkIBvPQPv11Vx3K3CWiDQA8oELgXS/VGmMOSnD+nj2Zj85I4Oc3HzaxkZz/6DO7Mwr5In/riVjZx4v33AGHVs2crlSU96TMzLI99NFQf0WKqpaLCJ3AjPwHL31mqquEpFxzvyXRaQ1nrBoDJSKyG+Abqq6QESmAkuAYmApMMFftRpjTs6wPvHHwsVbz/gm3PXuUoY+P5unhvdiSI82LlRnKuPPcS8JpX2fqampmp5uHRpjAkFObj6/fHsJy7blMubcFB4c3IWIcDvfOhCc8Zev2OuctLr9jd9QuH19RWPgJ8X+Dxtj/KJtbDTvjz2LG89qz4SZm7hh0gJ25xW6XVadl565j4P5RRUeSeULFirGGL+pHxHOX4Z15+nhvVi6NZfL/jWLxVv2u11WnbUy6wC3vL6IxLgGPDL0dOIrOYLvVNjuL2NMrViVc4A73lrC9gP5/OHn3RjVvz0i/vq9bMpbtzOPa16ZR4N6EXwwrv+xQ8JFZLGqpvpqO9ZTMcbUitPbNuGTO89hYKcW/OnjVdzz/jK70VctydxzmJETFxAZHsY7t59Z6TlGvmChYoypNU0aRDJxVCr3XdyZj5bncOWLc8ncc9jtskJadm4+IycuoKRUeXv0mbRv1tCv27NQMcbUqrAw4a4LOzH5ljR2HCzg8udn8/XqnW6XFZJ2HSxg5KvzOVhQxJu3ptGpVYzft2mhYoxxxXmdW/DJnefQvlkDRr+ZzlMzMigpDZ0xXrftP3yUGyYtYFdeIZNv6Uf3+Ca1sl0LFWOMaxLjGjB13Nlck5rI899t4ObXF7LPbvp1yg4WFDHqtYVk7j3CxFGpnNE+rta2baFijHFVVGQ4f7+6J49f2YMFm/dx+b9msyIr1+2ygtaRo8XcNnkRa7Yf5OUb+nJ2x+a1un0LFWNMQLg2rR1Tx/UH4OqX5vHewq0uVxR8CopKGPvvxSzesp9nr+3DBV1b1XoNFirGmIDRMyGWT+46hzNT4nho2koenLqcAj9d+DDUFJWUcuc7S5m1fg9/v6onP+/pzvXWLFSMMQElrmE9Jt+Sxl0XdGRKehZXvzyXbfuOuF1WQCspVe6dspyv1+zkz0NPZ3hqYtUr+YmFijEm4ISHCfcN6sLEUals2XuEy5+fzfcZu9wuKyCVliq/m7aST5bn8NCQrozqn+RqPRYqxpiAdVG3Vnxy5zm0bhzFLZMX8dw36ym1w46PUVX+/Olq3k/fxq8v6Mi48zq4XZKFijEmsCU1b8h/fjmAYb3jeeardYx+M50DR4rcLisgPP3lOibPzeTWAcncc3Fnt8sBLFSMMUEgul44z4zoxZ+Hns6s9bu5/PnZrMo54HZZrnrhuw08/90Gru2XyP9ddlrAXJzTQsUYExREhFH9k3hvTH8Ki0u48sW5fLg4y+2yXDF5zmaenJHB0N5t+esVPQImUMBCxRgTZM5o35RP7xpIn3ax3PfBcv4wfSWFxXXnsOMpi7bx8CerGdStFU8N70V4WOAEClioGGOCUIuY+rx125mMPTeFt+Zv5ZpX5rP9gP/uux4oPlmew0PTVjCwU3P+dX0fIgPw9syBV5ExxlRDRHgY4y89jZdG9mX9zjwue242czfscbssv/l69U7ueX8Zqe3jmHBjKvUjwt0uqUIRbhdgjDGnYkiPNnRqFcO4txZzw6QF/LxHG5Zs3U9ObgFtY6N5YHAXhvWJd7vMUzJ7/R5++c4SurVtzKSbU4muF5iBAtZTMcaEgI4tG/HRrwbQM6EJn6zYTnZuAYrnBlXjp61k+tJst0s8aemZ+7j9zXRSmjfkzVvTiImKdLukE7JQMcaEhIb1I9idV/iT6flFJTw5Y60LFZ26lVkHuOX1RbRpEsW/bzuT2Ab13C6pShYqxpiQkZNbUOH07NwCXp+zmf1BdK+WdTvzGPXaAhpHR/LW6DNpEVPf7ZKqxULFGBMy2sZGVzg9Mlx45JPVnPm3b/jVO0v437rdAX2Xycw9hxk5cQGR4WG8PfrMSj9XILKBemNMyHhgcBfGT1tJvtfl8qMjw3nsyh50bhXDB4u3MX1pNp+t2E7bJlFcfUYCV5+RSLtmDVys+njZufmMnLiA4pJSpoztT1Lzhm6XVCOiGrhpXVOpqamanp7udhnGGBdNX5rNkzMyyMnNr/Dor8LiEr5Zs4v3F21j5vrdqEL/lGaM6JfAkO5tiIp078iqXQcLGPHKPPYePsq7t59VK/eVF5HFqprqs/ezUDHG1FU5uflMW5LFlPQstu47QkxUBL/o1ZYRqYn0TGhSq5c/2X/4KNdMmEfW/nz+fVtard1X3kLlBCxUjDEno7RUWbB5Hx+kb+PzH7ZTUFRKl1YxDE9N4Io+8TRr5N9B8oMFRYx8dQEZO/OYfHO/Wr2vfFCFiohcAjwLhAMTVfXxcvO7Aq8DfYHfq+pTXvNigYlAd0CBW1V13om2Z6FijDlVBwuK+HT5dt5P38bybblEhgsXndaKEamJDOzUnAgfXxrlyNFibnptIUu35vLKjWdw4Wm1e1/5oAkVEQkH1gEXA1nAIuA6VV3ttUxLoD0wDNhfLlTeAGap6kQRqQc0UNXcE23TQsUY40sZO/L4IH0b05Zms+/wUVo1rs/VZyQw/IxEnwygFxSVcPub6czZsIfnruvDZT3b+qDqmgmmUOkPPKyqg53X4wFU9bEKln0YOFQWKiLSGFgOpGgNCrRQMcb4w9HiUr5du5Mp6Vl8n7GLUoW05DhGpCZyaY/WNKhX8wNpi0pKueOtJXy9ZidPXt3TtfvK+zpU/HlIcTywzet1FnBmNddNAXYDr4tIL2AxcLeqHi6/oIiMAcYAtGvX7pQKNsaYitSLCOOS7m24pHsbdh4sYOriLD5I38b9Hyzn4Y9XcXmvNgxPTaRPYmy1BvdLSpV7pyzn6zU7+fPQ010LFH/w58mPFbVsdXsdEXjGWV5S1T7AYeChihZU1QmqmqqqqS1atDi5So0xpppaNY7iVz/ryHf3n8+Usf25pHtrpi/N4coX53LxP2YyYebGCi8XU6a0VPndtJV8sjyH317SlVH9k2qv+Frgz55KFuAdvwlATg3WzVLVBc7rqVQSKsYY4wYRIS05jrTkOP50eTc+W7GdKenb+Nvna3nivxlc0LUlI1ITOb9LCz5dsf3YuTMN6odzuLCEuy7oyB3nd3D7Y/icP0NlEdBJRJKBbOBa4PrqrKiqO0Rkm4h0UdUM4EJgdVXrGWOMG2KiIrk2rR3XprVjw648pqRnMW1JFl+u3klM/XDyi0opdi4Lc7iwhPAwISXIzpSvLn8fUnwp8E88hxS/pqp/FZFxAKr6soi0BtKBxkApcAjopqoHRaQ3nkOK6wGbgFtUdf+JtmcD9caYQFFUUsp3a3fx63eXUlBc+pP58bHRzHnoAhcqO14wDdSjqp8Dn5eb9rLX8x14dotVtO4ywGcf1BhjalNkeBiDTm9NYQWBAp6z+UORXaXYGGP8qLIrDAfTlYdrwkLFGGP86IHBXYgud5HK6MhwHhjcxaWK/MsufW+MMX5UdoXkE105OZRYqBhjjJ8N6xMfsiFSnu3+MsYY4zMWKsYYY3zGQsUYY4zPWKgYY4zxGQsVY4wxPhNStxMWkQPAerfrKKcdsNXtIsppAhxwu4hyrJ2qz9qqeqydqud0VfXZmZihFioTVHWM23V4E5HdqhpQ1+S3dqqeQGwnsLaqLmun6vF1O4Xa7q9P3C6gAie8BbJLrJ2qJxDbCaytqsvaqXp82k4hFSqqGoj/wwKtq2vtVE0B2k5gbVVd1k7V49N2CqlQCVAT3C4gSFg7VZ+1VfVYO1WPT9sppMZUjDHGuMt6KsYYY3zGQsUYY4zPWKjUkIgkish3IrJGRFaJyN3O9DgR+UpE1jv/beq1zngR2SAiGSIyuIL3/FhEfqjNz+FvvmwnEaknIhNEZJ2IrBWRq9z4TP7i47a6TkRWisgKEfmviDR34zP5Q03bSUSaOcsfEpHny73XGU47bRCR50RE3PhM/uCrdhKRBiLymfNvbpWIPF6tAlTVHjV4AG2Avs7zGGAd0A14AnjImf4Q8HfneTdgOVAfSAY2AuFe73cl8A7wg9ufLVDbCXgEeNR5HgY0d/vzBWJb4bmVxa6y9nHWf9jtz+diOzUEzgHGAc+Xe6+FQH9AgC+AIW5/vkBrJ6AB8DPneT1gVnXayXoqNaSq21V1ifM8D1gDxANDgTecxd4AhjnPhwLvqWqhqm4GNgBpACLSCLgXeLT2PkHt8GU7AbcCjznvVaqqe2rnU9QOH7aVOI+Gzi/vxkBOrX0QP6tpO6nqYVWdDRR4v4+ItAEaq+o89XxjvsmPbRv0fNVOqnpEVb9znh8FlgAJVW3fQuUUiEgS0AdYALRS1e3g+Z8KtHQWiwe2ea2W5UwD+AvwNHCkFsp1zam0k4jEOq//IiJLROQDEWlVK4W74FTaSlWLgDuAlXjCpBswqVYKr2XVbKfKxONpszLe/yZDyim2k/f7xAKXA99UtayFyklyehkfAr9R1YMnWrSCaSoivYGOqvofvxQYIE61nfDs0kkA5qhqX2Ae8JTPCw0APvibisQTKn2AtsAKYLzPC3VZDdqp0reoYFrInVvhg3Yqe58I4F3gOVXdVNXyFionwfnH+yHwtqpOcybvdLrVZd3rXc70LCDRa/UEPL8i+wNniEgmMBvoLCLf+7/62uOjdtqLpydXFr4fAH39XHqt81Fb9QZQ1Y3Obp0pwNm1UH6tqWE7VSaL43fjlLVfyPBRO5WZAKxX1X9WZ2ELlRpy9lVPAtao6jNesz4GbnKe3wR85DX9WhGpLyLJQCdgoaq+pKptVTUJzyDZOlU9vzY+Q23wYTspnuslne8sdyGw2s/l1ypftRWQDXQTkbKLA16MZ396SDiJdqqQs+snT0TOct5zVFXrBBNftZPzXo/iubLyb6pdgNtHKgTbA08AKJ5dC8ucx6VAMzz7G9c7/43zWuf3eI7QyaCCoyeAJELv6C+ftRPQHpjpvNc3QDu3P18At9U4PEGyAk8YN3P787ncTpnAPuAQnh5KN2d6KvCD04bP41xdJBQevmonPD04df6eyt5ndFXbt8u0GGOM8Rnb/WWMMcZnLFSMMcb4jIWKMcYYn7FQMcYY4zMWKsYYY3zGQsUYY4zPWKgYY4zxmf8HKD2wA0OXy74AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "eviyr = evidf.groupby(pd.Grouper(freq='1Y')).aggregate(np.mean)\n",
    "\n",
    "# plot\n",
    "ax = eviyr.plot(marker='o')\n",
    "ax.set(ylabel='EVI')\n",
    "ax.set(title='Annual Mean EVI (2003-2012)')\n",
    "\n",
    "print('Plot: ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_Yyg1DzbaALa"
   },
   "source": [
    "The trend of the annual mean EVI values decreased significantly from 2007 to 2011 owing to tree mortality caused by bark beetles. The recovery started in 2012.\n",
    "\n",
    "### 3. Download Tree Mortality from Bark Beetles with SDAT Web Service\n",
    "\n",
    "We will use the [Tree Mortality from Fires and Bark Beetles at 1-km Resolution, Western USA, 2003-2012](https://doi.org/10.3334/ORNLDAAC/1512) dataset. This dataset provides annual estimates of tree mortality (as amount of aboveground carbon in the trees killed) due to fires and bark beetles from 2003 to 2012 on forestland in the continental western United States. The data is also available via [SDAT](https://webmap.ornl.gov/ogc/dataset.jsp?ds_id=1512.). \n",
    "\n",
    "If you look up the URLs of each granules, you will notice that the bark beetle data from 2003-2012 have 'dg_id' of '1512_1', '1512_3', '1512_5',...,'1512_17', '1512_19'. These are the layer names we will need when accessing the SDAT service.\n",
    "\n",
    "Let's first connect to the ORNL DAAC WMS and retrieve its Capabilities document."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Lbl42yvETcPL"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18337 layers found from ORNL DAAC WMS Server\n"
     ]
    }
   ],
   "source": [
    "sdatwms = WebMapService('https://webmap.ornl.gov/ogcbroker/wms')\n",
    "print(str(len(sdatwms.contents)) + ' layers found from ' + sdatwms.identification.title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hrruWjsJC980"
   },
   "source": [
    "List WMS layers that belong to dataset number 1512: Tree Mortality from Fires and Bark Beetles at 1-km Resolution, Western USA, 2003-2012."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "BnAClc5I1IWY"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Layer Name: 1512_1 Title: Tree Mortality from Bark Beetles in the Western USA 2003\n",
      "Layer Name: 1512_2 Title: Tree Mortality from Bark Beetles in the Western USA 2003, Standard Deviation\n",
      "Layer Name: 1512_3 Title: Tree Mortality from Bark Beetles in the Western USA 2004\n",
      "Layer Name: 1512_4 Title: Tree Mortality from Bark Beetles in the Western USA 2004, Standard Deviation\n",
      "Layer Name: 1512_5 Title: Tree Mortality from Bark Beetles in the Western USA 2005\n",
      "Layer Name: 1512_6 Title: Tree Mortality from Bark Beetles in the Western USA 2005, Standard Deviation\n",
      "Layer Name: 1512_7 Title: Tree Mortality from Bark Beetles in the Western USA 2006\n",
      "Layer Name: 1512_8 Title: Tree Mortality from Bark Beetles in the Western USA 2006, Standard Deviation\n",
      "Layer Name: 1512_9 Title: Tree Mortality from Bark Beetles in the Western USA 2007\n",
      "Layer Name: 1512_10 Title: Tree Mortality from Bark Beetles in the Western USA 2007, Standard Deviation\n",
      "Layer Name: 1512_11 Title: Tree Mortality from Bark Beetles in the Western USA 2008\n",
      "Layer Name: 1512_12 Title: Tree Mortality from Bark Beetles in the Western USA 2008, Standard Deviation\n",
      "Layer Name: 1512_13 Title: Tree Mortality from Bark Beetles in the Western USA 2009\n",
      "Layer Name: 1512_14 Title: Tree Mortality from Bark Beetles in the Western USA 2009, Standard Deviation\n",
      "Layer Name: 1512_15 Title: Tree Mortality from Bark Beetles in the Western USA 2010\n",
      "Layer Name: 1512_16 Title: Tree Mortality from Bark Beetles in the Western USA 2010, Standard Deviation\n",
      "Layer Name: 1512_17 Title: Tree Mortality from Bark Beetles in the Western USA 2011\n",
      "Layer Name: 1512_18 Title: Tree Mortality from Bark Beetles in the Western USA 2011, Standard Deviation\n",
      "Layer Name: 1512_19 Title: Tree Mortality from Bark Beetles in the Western USA 2012\n",
      "Layer Name: 1512_20 Title: Tree Mortality from Bark Beetles in the Western USA 2012, Standard Deviation\n",
      "Layer Name: 1512_21 Title: Tree Mortality from Fires in the Western USA 2003\n",
      "Layer Name: 1512_22 Title: Tree Mortality from Fires in the Western USA 2003, Standard Deviation\n",
      "Layer Name: 1512_23 Title: Tree Mortality from Fires in the Western USA 2004\n",
      "Layer Name: 1512_24 Title: Tree Mortality from Fires in the Western USA 2004, Standard Deviation\n",
      "Layer Name: 1512_25 Title: Tree Mortality from Fires in the Western USA 2005\n",
      "Layer Name: 1512_26 Title: Tree Mortality from Fires in the Western USA 2005, Standard Deviation\n",
      "Layer Name: 1512_27 Title: Tree Mortality from Fires in the Western USA 2006\n",
      "Layer Name: 1512_28 Title: Tree Mortality from Fires in the Western USA 2006, Standard Deviation\n",
      "Layer Name: 1512_29 Title: Tree Mortality from Fires in the Western USA 2007\n",
      "Layer Name: 1512_30 Title: Tree Mortality from Fires in the Western USA 2007, Standard Deviation\n",
      "Layer Name: 1512_31 Title: Tree Mortality from Fires in the Western USA 2008\n",
      "Layer Name: 1512_32 Title: Tree Mortality from Fires in the Western USA 2008, Standard Deviation\n",
      "Layer Name: 1512_33 Title: Tree Mortality from Fires in the Western USA 2009\n",
      "Layer Name: 1512_34 Title: Tree Mortality from Fires in the Western USA 2009, Standard Deviation\n",
      "Layer Name: 1512_35 Title: Tree Mortality from Fires in the Western USA 2010\n",
      "Layer Name: 1512_36 Title: Tree Mortality from Fires in the Western USA 2010, Standard Deviation\n",
      "Layer Name: 1512_37 Title: Tree Mortality from Fires in the Western USA 2011\n",
      "Layer Name: 1512_38 Title: Tree Mortality from Fires in the Western USA 2011, Standard Deviation\n",
      "Layer Name: 1512_39 Title: Tree Mortality from Fires in the Western USA 2012\n",
      "Layer Name: 1512_40 Title: Tree Mortality from Fires in the Western USA 2012, Standard Deviation\n"
     ]
    }
   ],
   "source": [
    "# define a filter function to find layers whose name start with '1512_'\n",
    "def myFilter(x):\n",
    "  if x.startswith('1512_'):\n",
    "    return True\n",
    "  else:\n",
    "    return False\n",
    "\n",
    "# filter layers\n",
    "myLayers = filter(myFilter, sdatwms.contents)\n",
    "for l in myLayers:\n",
    "  print(sdatwms[l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "capJY1rgDyOn"
   },
   "source": [
    "Print out some metadata associated with layer 1512_1: Tree Mortality from Bark Beetles in the Western USA 2003."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ecuwwIKeD4Ma"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Layer:\n",
      "Tree Mortality from Bark Beetles in the Western USA 2003\n",
      "Bounding Box:\n",
      "(-2535990.0, 951253.0, -428992.0, 3239250.0, 'epsg:5071')\n",
      "Bounding Box (WGS84):\n",
      "(-130.435, 28.4705, -100.538, 52.0938)\n"
     ]
    }
   ],
   "source": [
    "print('Layer:')\n",
    "print(sdatwms['1512_1'].title)\n",
    "print('Bounding Box:')\n",
    "print(sdatwms['1512_1'].boundingBox)\n",
    "print('Bounding Box (WGS84):')\n",
    "print(sdatwms['1512_1'].boundingBoxWGS84)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ektG_GHaaMNW"
   },
   "source": [
    "Then check what methods and format options are available for the WMS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "285cNKgXaGAT"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['image/png',\n",
       " 'image/gif',\n",
       " 'application/x-pdf',\n",
       " 'image/jpeg',\n",
       " 'image/png; mode=8bit',\n",
       " 'image/svg+xml',\n",
       " 'image/tiff',\n",
       " 'application/vnd.google-earth.kml+xml',\n",
       " 'application/vnd.google-earth.kmz']"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[op.name for op in sdatwms.operations]\n",
    "sdatwms.getOperationByName('GetMap').formatOptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uIUh5vgkaURz"
   },
   "source": [
    "We will use the WMS \"GetMap\" operation to download the images in PNG format for years 2003-2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yBq-X2_paSAC"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2003\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABJ0lEQVR4nO3RsQ2DQBREQVqhXddGTzi3RIb03/kmePlq5zg/5/3v3de1TMf0WVCApAMSC0gsILGAxAISaxuQVVCAxAISC0gsILGAxAISayuQFVCAxAISC0gsILGAxAISazuQOgqQWEBiAYkFJBaQWEBiAYm1JUgZBUgsILGAxAISC0gsILG2BamiAIkFJBaQWEBiAYkFJNbWIEUUIAEEIEC6TQMAAdJuGgAIkH7TCECAtJtGAAKk3TQCECDtphGABECetgB5AeXNHUB+QKZ3AIkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWEBiAYkFJBaQWF8dBw3fr4+2ZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2004\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAAA/0lEQVR4nO3RMQ0AMAzAsFIp3ZHeYCyHD/+RMnv20jG/AzAkzZAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpCYBwMiDcsRaP5SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2005\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABMUlEQVR4nO3bsRHCQBAEQaWidIlNOYkAKBlQfN0g2hj/ftv+bX/spzpt0wcISDogsYDEAhILSCwgsYDEAhILSCwgsYDEAjLYeRwvAVk88LsBWTgukNi4fwEyPditQaYfX+xrINMPuUuXINOH/Wvb9AECkg5ILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCEgtILCCxgMQCMtTV72cgC8f9JCALxwXywxBAgLSbhgACpN00BJAoChAg3aYxgABpN40BBEi7aQwgQRQgQNoBiQUkFpBYQGIBCQYkFpBYQGIBiQUkFpBYQGIBCQYkFpBYQGIBiQUkFpBYQIIBiQUkFpBYQGJNgDwBdFcrK0TZikUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2006\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABKElEQVR4nO3coRHDQBAEQaWidJ2PkYOTExARcN28v8Hwq21+x/V5X//e+TqX6ZgeCwqQdEBiAYkFJBaQWEBibQOyCgqQWEBiAYkFJBaQwe5uBPLjgZ+2FcgTlKn7tgOpByQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkKF8AxoY13um2LhAFoYAAqTdNAQQIO2mIYBEUYAA6TaNAQRIu2kMIEDaTWMACaIAAdIOSCwgsYDEAhILSDAgsYDEAhILSCwgsYDEAhILSDAgsYDEAhILSCwgsYDEAhIMSCwgsYDEAhJrAuQLIi7fHbzci3cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2007\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABLElEQVR4nO3bsRGDMBQFQfXfD5FboCdcgIfAHjQ6wwaXf72NNY59P9RprD5AQNIBiQUkFpBYQGIBiQUkFpBYQGIBiQUkFpCVvbaPgEwe+NuATBwXSGzcZ4AERrsvSGCAWteBBB5zh85BAsc9sbH6AAFJByQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgGJBSQWkFhAYgFZ1cnvZyATx/0lIBPHBfLHEECAxAtgAAESLoABJIgCBEg4ILGAxAISC0gsIMGAxAISC0gsILGAxAISC0gwILGAxAISC0gsILGAxAISC0gwILGAxAISC0gsILGAxAISDEgsILGAxAISawHIG3hS09LgdqI0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2008\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABW0lEQVR4nO3ZMXKDMABFQc6fG6RJnwOkyhV8J7t2DAUZM3qgLV6P/jYasdw/P+6X7/fnNC3Dx4ICpByQWEBiAYkFJBaQWPOAnAQFSCwgsYDEAhILyMi+v14CcvDAe5sLZA/KG8YFstWgcecACYx2XZDAALXeBxI4zBXaBgl83IwtI28UApIPSCwgsYDEAhILSKxl9GOagKQDEgtILCCxgMQCEmsp/EcWkGxAYgGJBSQWkFhAYj2DQBkekFhAYgGJBSQWkFhAYr2CQAEiINmAxAISC0gsILGAxFoHgQJkum631YAcOO5/AnLguED+tnW2tQZDnBNkz8B7C2CMBTly3EuC7EEZPeY0IDMVwAASRAECJByQWEBiAYkFJBaQYEBiAYkFJBaQWEBiAYkFJBiQWEBiAYkFJBaQWEBiAYkFJBiQWEBiAYkFJBaQWEBiAQkGJBaQWEBiAYk1AOQBayfMJvybTq4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2009\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABaElEQVR4nO3WwUmDURBG0RRrD4It2IN7C3BlZyLEAlyYHxLnJnMWt4HvwLx3+np5Oj9658+Pu+k0PRYUIOmAxAISC0gsILGAxFoDci8oQGIBiQUkFpBYQCZ7f/sVkBsPfLRVIIdQrjAukEsQhoZ+WJDCWXlYkMpDWupqIOM/lkcHGR9uK8j4QECWg8RRgMQCEgtILCCxgMQCEgtIrJ0gYRQgsYDEAhILSCwgsYDE2gsSRQESC0gsILGAxAISC0is3SBBFCABBCBAwgUQgAAJF0AAAiReAAIIkHABCCBAwgUggAAJF4AAAiReAAMIkP/t+/X54qYh7hLkyMBHm4YYB7nluOtBpsdcAzI9EpDFIBUUIEC6TWMAAdJuGuNPkG0o0xhAgPQDEgtILCCxgMTKg2xDARILSDAgsYDEAhILSCwgsfIg21CAxAISDEgsILGAxAISKw+yDQVILCDBgMQCEgtILCCxJkB+AGc+0wFR6yagAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAAA/0lEQVR4nO3RMQ0AMAzAsFIp3ZHeYCyHD/+RMnv20jG/AzAkzZAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpAYQ2IMiTEkxpCYBwMiDcsRaP5SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2011\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkCAYAAABw4pVUAAABS0lEQVR4nO3awW0CURQEQVIhXcdGTjgASzZIRtPw69D3t1OX1Wov16/rXZ0u6wMEJB2QWEBiAYkFJBaQWEBiAYkFJBaQWEBiARl2v91+BOTFAz8bkBeOCyQ8NBAg79F6cCBA2q0HBwKk33p0IEDarUcHAqTdenQgQNqtRwcCpN96eCBA2q2HBwKk3Xp4IEDarYcH8iYoQGIBiQUkFpBYQGIBiXU0SBEFSAABCJBuawAgQNqtAYAA6bdGAAKk3RoBCJB2awQgQNqtEYAA6beGAALk3JE/FmQ91keDrB+82r+ArB/ik/oVZH3ciQGJBSQWkFhAYgGJ9edr7/rA0wISC0gsILGAxAISC0ish772ro88KSCxgMQCEgtILCCxgMQCEuvhv07Wh54SkFhAYgGJBSQWkFhAYj31s/X62BMCEgtILCCxgMQCEgtIrG9rUUAspvTOCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# specify x/y spatial resolution in decimal degrees\n",
    "resx = 0.010 \n",
    "resy = 0.010\n",
    "\n",
    "# The date range is 2003-2011\n",
    "mortDates = [2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011] \n",
    "\n",
    "# odd numbered array starting from 1 to 17\n",
    "odd_i = range(1, 18, 2)\n",
    "for i, j in enumerate(odd_i):\n",
    "    layername = '1512_' + str(j)\n",
    "    img = sdatwms.getmap(layers=[layername],\n",
    "                     srs='EPSG:4326', # WGS 84\n",
    "                     bbox=(longitude-resx,latitude-resy,longitude+resx,latitude+resy), \n",
    "                     size=(100, 100), \n",
    "                     format='image/png', \n",
    "                     transparent=True)\n",
    "    print(mortDates[i])\n",
    "    display(Image(img.read()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4jUUPYjsaac_"
   },
   "source": [
    "In the above images, higher tree mortality is displayed in darker red and no mortality in green. Mortality peaked between 2008 and 2009.\n",
    "\n",
    "We need a legend for the maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "OCvh1dTBD047"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKEAAADaCAYAAADOtCK0AAAgAElEQVR4nO2dCZQkRbWG7wFENlF2VJBmlU1WQTZhAAdwRFGQYVMZFBVQ1CebjNs8FxxAcAFEQLBxgYe2ooIi8tBGBR4Iok+U1ecg6mMRRXgi4Iz98iMiTt2OjsiqyqzJru66/zlxTmVGZWRk5V+x3PvHDRGDwWAwGAwGg8Fg6EcsXLhw9qJFi0Z6kSir5FabFmnrIm1VpHWKtEQjD+iwcpE2TpynLst1cP0mRTqhxv03FPfs8TOv788vVaNs6r9vkfjtNyr53jFFelmN+5Th+UV6XZH2LtKKXV8NecZ6BMoqudVdRXq4SL8v0v8V6bdF+mCRntVhVYeL9IquH9DhkCL9MnH+KXEkaIfdinSV/7x8kX4tndcbXFuksSLtoM4tWaQ/+vMrdVGWxluL9EiRbi/SD4v0aJG+kvnud4v0xor3KcNbivR0kW4p0r1F+pO4BqdzNEzCV/vPvMBXFun+Ip3eYVWvE/dvq4K6JNR4jjjidEtCiHKmOrdLkX4l1Um4q7j67yOtFnaFIp2c+f7iIuH2RXqB/8wf65IifbOrEiaJhAH7FemvRVraH58q7p9ESzkqrS702CL9vUgLxJHpbeJapCvF/fv/XKQvFem5mXt3QsLzfbk3FumxIn25SM/2eTtKq4X5hjji/Lcvkzqu6uuYIyYk/JC4XmBJf+6sIn1AxpPwheJa3CeK9ANxhDoxUybP+4VMXgqahNyP3+5Qf/xTn/c/RfpLkY4u0pZF+oW49/PvXdzncHG/TeeYZBIynqEp39kfv17c+ILWZl6RbvDnGWdAjjerfNIBPo9z/KinZO7dCQkhCt0J46Z1i/QbcT8omFmkm/znIXHEWdvflz/QS/y5ZTP3p+wj/PPsJG4M+AdxBNYk5Pf7apFWETe+erxIn8mU+WCRDs7kpRBISNl0nRArtKB/K9I1/tloWfldvu+fi9+HRqHTLpY/6wVd1GvSSQgejM7TmvHw2xVpYZGW8edz3fGa4iYYtGI3JPJBpyTUkw/Gq2f5z5qEqe54OV/OkpJGIOG7i/TpIu0ubgy3krRIyHP+q0hrqev+Q/IkXOTL6RSQ8D1Fuk1cC6zxt6gs/oxHqmMIeVgH9+A90pKu3UW9Jp2EtCJPihvfgE+Km7z8pEhXi3sp4YFiEjL+4eXSxf2oSP8lbrKTAiS8PTpHKwDJNQkPUvnvKtLF/nM7ErZDICHdLePg84p0lIwnIa3vUzJ+Bv1xyZOQl/36LuoACenmv5XIg4Trq+Ofi2uJA74mrovmj3KTSsur72zt6zSzizo5TDIJ9xTX1NONbVakB6Q1xedc6PZATMJ3FOkKab00frQcCWf5sjUCAdb1xxDlQJXPGK/XJAQ/LtI/irSajCfh8/xn/WI/J3kSMjY9M5OXAiSkdedP8KYoDxKuq45vlfFkukwcCWnpd1cpmJaGxI3lj5AqmEQSYnujdTrPH/NQjMMCqWi9NAkh3Bx1PWOaz6ljysmRkBdMq7eTOodp4T51v05JyA9PV7iq+i5DiNdJ3t6nSbirvzfQJARMBI5SeZhwciR8rbg/8GbqHCTJdZthTMjY7gF/fUCnJEwBG+wdRfpwJr89GjRWQ0JIR0twp7gZ6Oel9c+nheFH/6Y/z/iOrjqQkJfD+JEJCN3mS8XNmM/11zBxyZEQfETcQJ+BP60IEyJtsuiUhOLrwLN8p0jrSecTkxgxCekZmI0yHLm7SN8TN0RJgT8PE7EnfX0u8tdcm/m+nh3z2z3k7wfqkPBj/hlullY3nery+wJ7iGsJse7TIj078R3sTbx8BsV0y3Sj2qPBv/hV0vIMMA45rkj7i/tH7inlgCzHixvvDUV5O0rL3gXWk9Z4cfUivVzl0RXv7J+HVva5/nOuJaTsFyXOL+2vW1qdo5veq0hriCPWv5U8D2CcydAEU87eJXVg1q8nDPyWYTKCzVYPA2itV1fH28rE30uX8+ooVXUqGPoADFE2EEck/rS0UGWuOIOh56AFYeLwT3FmkmoDfYPBYDAYDAaDwWAwdImh+UOzizTSo1RmrMa2hkAAIQHGYjM9GBwgz7qnrTvWi0RZJbd6nzhFMirlT4gzRSxd8n3DoKAhEuJiwqU2Sx3jYtqviWc09DkaIiEuIPyLWvmMKnj+Yn9AQ/+jIRLis0UwoLVyEPCLi/0BDf2PhkiI8vnJ6ByCzS8v3qczTAk0RELUKXTHWhGDBCsnUzIMEhoiIaoQtIBaDoVmbXEsQTRMNTRooqHVQ3DJwiTk5SzVXKGJZzT0ORo0ViNSZUbM2BBV8m6NPKDBYDAYDAaDwWAwGAxtsfCmm2Yv+tnPRnqRKKvkVoQQIwIAISUOLPleJ8AHTfQBbI2EADkkymfVGmtfWUzOksk6QSgNixuQZ+yWW8Z6kSir5FaEMoMYLIL/YM1q43MmnBnGb8KbEWMlmHxYgE6EASI07OHv956a9zMsTjRIwgBW89clIXH0Zqnjr4uL/wdojXV8PBbGE6aiyfDEhm4wRUlIAMnj/GcU26i1Q8BHuulz1XeJGYPf+nk172lYXJiiJBwSF2ePLhmx7KekFRuQWDMfUd/lPCR8cc17GhYX+pCEtFwzfdox8x2CPNLlHuE/E5ZsC59HYMl56rshDqGOXmXoJ/QhCZlFX+/TxYl8Jh6EydhSnSP66Zf858/6FMBMmpZQh3Iz9BP6kITtEEilI0xRXogYT+y/61UeLSrxoW1i0q9okISMySAEEVcv9p9fUPL9MmD/u1Bc60aLuEBcuDfABIQ4hAf78omH+LGK9zE0gQaN1cxmfxmlfStWG2M02kQilf6vuBmxDt9LdCvMMsjGdCBOg8FgMBgMBoPBYDAYElh41bdnL7r6ipFeJMoquRWhQJipslsQ+6i9KspnWwN2cUKexR52dex67Dz0Xl8WOz29TZ1nls6WXrj92FTx+dG12CE/Km72jnmn6s6ihk4Becau+e5YLxJlldyK3SvZS4PtI94ubgut4EqDoOxrguSK7QzYz6SbzQNjoKphPw10i2y99VZ/HgM3e4RAcrY64M/AnyL4nTHzQFr2cnuNOFGE3mrMsDjQIAnR/G2jjiFDaKEwNP9A5bFL+Y8qPhJ7j2CsTqlmIJtuYdnxEu9L2NeNVu930t2WYYa6aJCECA7CllcYj3Gl7eGP2ahZb0u1uThDc5Uu+SRxW7ayyzq7ar6s5LtsXAMJh/xxCE0y019roocm0CAJIRbbieG2Y2fOuTJ+X7lj1HeJ0hCHkusUuATv9mlY3M6TOV81+8Zdo46pG0psxosIImi9480IDb1GgyQ8W9x+dbxU9mSDiEM+b1RaGwuC0EKtkiiHMWWQeq2cyB8WtxPSyur7jD9XjL7HJOohGb9dFvVAnxi2PGNM+MfsExl6g4ZIyP5rkGpNdQ7d32n+M4uf5qo8xA7ItVLdMTuLB6nXNol84h6OqmPKeDL6LmtT2MhQ7/oZ6jSsjkOLnCK7oVdoiISQig20l1HnzpHWVrHvFxdMPYBoXbdWfCRm1XqXd4Iu6X2N2RDwz+LWnsQgrva31TFCCbaVtdjaixMNkZCXyESEVopWBcU0RAh2RQiCImZv/5lJzDsrPtJyvuw54qRejPuCGYZ7s8cxJpyZKoWtXtfy9cBMhNmI5+nbbVOnDRo0VqOYZkzIYB9CMlnQ3e0B4nSB7GFMN13HTMLMmJaULh3D9Cb+/IYyUU5G0l01RnQmJ4Q3hoBVNY8Gg8FgMBgMBoPBYJgyWHjJhbMXXXrhSC8SZbW5HXKty8XFrMZAvZo/j8lmWNxM9T/FKV9i7ClO5ID5Btdc1dkz98JPjZkm3taM+qC4uahIF0R5G4izZcZpa5+/ry/zV0W6ski7VKzf4AHyjF120VgvEmWV3ArvBHItTDPYAxEahLXDuPTYdBHzCKHccLPto65FkvWItKRehPl4tlQDXhHWKN8rbtdRDe5/oziDdeyuW6NIx6p0hjhD9uo+H3POyb5+uCXx0mxcsY6DhQZJiGB1biZvyegYsaluibDbHSq9BUSLSRiwu7T3GfNnulodx8+AIPY4MbRHQyQkdActw7ZFOrJIhxfpOSXVYquJEFlrSJzLb50iHS3Oy7JM+rKuUIeEGNl/I84rkwPBmt5dqWaDhoZIiMcC/y0vju4QQStek5RKBoEprrPt/THdG904MQmRVyGzIgJDXX9uHRIS9YEhQy7cHNdT53jpgCGFhkjI2AgShu6YlgTtXqzzQ27F5ESHdtvHXxu6Y8aCyK1SW5Lh+w3+4K3aPHodEjJ+vTyT9yJxUSF6PXyYvmiIhLQYEGkHdW6euBlyAC0bazsulfE+5a1kogxsWFoyMA3WhQSZV7u9lKuSkLqxBCBlCUAIwe7289rc26DREAl5cQvEda0BqFvOU/mYRWgd4/EeUixEDZuocxDo5JqPXpWEmHgYLsTxbaj3qLh42hYBrBs0ODum62WwTpdJi4jsfm+fB6FYnIQpJnSnW6tricCFDQ5zCET+h1Rf/7GJLx9TzAX+c2hll/fHx4szCfF5++h6/jxfTZSL7RKZ2Cz1DBYdthM0aKympSCsL60aY6bjVR52wlhedbbKRwc4LG4ywASlzrqPuYl77eXz1knk6UCdtHBXqe9rXJm41kw0BoPBYDAYDAaDwTAlsPDc+bMXnTt/pBeJskpuVbbBIu47fMqYYoajPGasKQlVMJ3gUWFHACReBFyaUfW3kHpSLjwlqGu+XKSPJ8reTpxtEpcjrsuNatRzegHyjJ132lgvEmWV3Kpsg0UCnRNFC4nVI1Eeyza1hApPCB6UtXw+9kM21MF+yKaKrJLbQqqhjpSLSGOj4vSQ10fX4klhsf08cbFxUAnheoyVN4OJBkkYULaPCd6ImIQx8HJcp47jF4nw9f0d1KMMVaRcAewyFZNwV3EhR4InBaM4f6QX1ajj9MEUJCEL2Y8qyaclOqmDepShjpQrRULIxnrrsIaZrpl4OcvVquV0wRQjIW4w9qlbLZO/szi/bt0Wpo6UK0VCwMbgyNdGxQWD2rNeFacRphgJWRtyVSaPVoaX/OaS64mHGPy6Za1QHSlXioQIeJmQINZ9g7+WYYONCcEUIiFdIPENU35jdIi85HbyrTBpIJWF96gq5QIpEvLHuE1aY0Jkaw9LK0joYGMKkRBTCEsE4sCZvFBaFaK99qplqSLlCkiREJn/tdE5ZsevqVrBaYUGSVi2weKy0grR+7j/vEN0/akyPnxcAJpECLOvtLraTRLf6wR1pFxhn2bCDf/af36JzyPgEqYjSMc4kljdTFRM/g8aNFaXbbC4ZiLv0uj6b0t6Q8ZvJq7NreprhzpSrh0T135C5dMls8yU1py1yfEWGgaDwWAwGAwGg8FgmDQ8Nf+42U+fesJILxJlldyqTMrFRjqxROpclY9RF3kUJqDdomuxDTLzxjbHvnQoaqpGZyCSFmofZGFssqOjdqGEwWPDisFbxKloYpsls/zPipv9cv3LE/eg/jzfdhXrOP0AeRaeftJYLxJlldyqTMqFH1XLpJB1Xabyr/HHrNKLIy8gocJdh8kD0wlG4BOr/RrPSKw+Jm7ZJmodbHu7+zwIxYo6IiuwBx5CCu26g5B4UrAhIivDJLO3jAd/Djw7hAh5dcU6Tj80SMKAMo8JoGWDVKktXiFnTEI8JTpsCK3i9zqoR+7eGvxxzsrkYYBGTxgW6xPSjkhcZQvf2f72dHF/FCNhQB+SEO0dMqdU5K0UCSEgXXx4+aizP9tBPToBu0d9KpNHwE9ayhAnESM2u5Xy52EPlqHo+3TV94hrMY2EGn1IQsaCw5m8FAlxteE1YRxHF4kfORctqxvg+mOh/paJPAjP0EIvASD6AkMNIo4h38czogUKw9Lad9lIqNFnJCQEMLsxxWOpgBQJGQvSfdMKEb8QMuQ27O5UygWJEa3m3H+MF9lFdCV1boE4v3hokbk2CBlw6fEHWcofGwk1+oyELFpC4pSLR50i4c3iJicBEPDuzPWdSLnoXtkFKhfYKOwOumF0npn5PHWMAOMx/5nZMC3nB31ivQkTmHimP5joMxIiFvhcJg+kSHiHjDf50Mr9qYN6pMDkg4kOs/GUmQfSPCquZYvxeXEKmwBmyPf5z5igzlEJcjJ5MhEDaJCEZVIugJyLl5OyrW3pv4+o9VQZvykiL58989iYkUDstHaXVPw5WK0HyQ6QVrcdxoTr+7yzVd5MX2+AHRT5Fy0gKwGxJ346cx/rjjUaNFaXSbkAhuJRSQtTz0hcGzZFJGIXs1gIAElYArpGxZ/j/MR9zvB5uyXySDp4J7G4GZ+iFRyWiRt9B1wp1hUbDAaDwWAwGAwGQx/hiROPmP2P971lpBeJsjK3webGvnSYT1h/i393LZXfqw0WkXvF0bS6QZmUCxBJiwVYKGFYeKXlWO2icjFj/4IvmxV5ZTtaDRYgz9NzjxzrRaKszG0I20G0K8K/YZYhhBs7fQYi9WKDRSJ74W1pF8ioDGVSLkxH2PewS77M1wHPR7BXlkXlEp+H/ZI/HGuQh2vUc3qhIRLyArULDPsZ+9VtpvI1ut1gMej0UK/UIWGZlIuWjkhawefM8+C+2zW6JrX4HUM9goawYH49f7ySGBojYQwEApAwt89Htxss0nrincCDUYeEMbSUC/IhMdvWHzOcwCgdB19KkZDQJTdG5zBqz+hVRac0JomEuN6QPKUEAt1usIjHgm4Sz0kvSZiScjEuxR886u+Tap1TJKRr/050jmc4qBcVnfKYBBIy2EeutUEir8oGi6hdjvGf25GwjpQLGRYTIyT9RNb6gq9rPMFIkfB94lQ0GrcW6bCSOgwOGiYhA3oG8ym/aZUNFplxMolBeIAyh1nz4/5zavZZR8oFgR+WVgtMHq1ZHIouF5Xrx9E51qO8IlOHwUKDJER/x0D+kERe1Q0W6YK1RIou7wn/OV4J1wnKpFy04L+NzjHLfVd0LkVCxpEMKYI1gD8PwT4tIBJoiIR0s7xAyKNlUKv6/F5tsFh3TFgm5XqBvy8RteiuibD1tKpnWVQu/mSYpD4qbkbMH+7KGvWcXmjIWJ2KdkUK4tBebbBIeddU/jHKpVwADSBrijGvEGFrTnTvsqhcm4qbjEHcK2S8sd5gMBgMBoPBYDAYDJOKPpFyAWa82AmJeEW8Fm1sLovKBTCHMPPEP1snIj4zXKI5MAPGuIwNUIsazvb3wEg9LONDfWBHjCOLaVslZhpspKzCYyG8njkPNvpEykXEBYy5kBhVys9k/Prjsqhcm4uzMb7JX4+db6eKPwc2yHf5+lBX3ItHq/wzi7S/fwZsnixBDV4V/MMYqXV0Me0aZL89vCR4TzD1nFyxjtMPfSLlOkWcOCEAY/FvEuWkFr9DjPPVMe68L3b4+O2A8Tq1bQXA6Iw7cW1/DAnPy3wXCRe2xW0y+YONPpFy0ToQRiPEasH3G6tOQIqELHw/XB3TSt3ZYT3KQF0QHZySyUdR8xdpLX6HhAw3GHbQNWu3Hy0rxm1cl7SQiGZtS7GAPpFy8cLx9xI6DVIxLkx5FFIkxB04Sx0TAeGvHdYjBVxtxLeh6+ePEKttcLkxZqSr1tEiiD3IMxGmDq0gw40gWkXveL84GRitNF6fXAs7eOgTKRex/ogfwziJVpGu+ITEtSkSLpDxSwFY9/F45t6dSLkgDqRGpQ0R3xDlM97Ex82YkHFuGNcupb6D6AL/8XH+mNaRrjt0x/jA/y4WMtihT6RctB56tkhrhlhgqeh7uahcWlwaur4UOt1gMYDJzq8yeXS3qIJyW8YiVgjhhIkU9miUz/CjbEfSwUGfSLmIcnpC9F1ajmWj76VIyFoULTJgPNmrro4JUizfCmAoQRebC2xEvS70n5GhIVzQMjXWzaRCIg8e+kTKhfKY8SDmFroq1nbcoK4vi8rFMk26d5YDbCFuPJbaA68ToAWkFX6huPUfDAtO83mMFVFCD4lb+olki9ZtBZ+PKYfuljUntHBP+3oGXO/L4plZmUePsHLFek4v9ImUixaCVW0PiFtfgoxejxnLonKBo8S1ssxW2aGzLHh5GYj6v0Cc4JRJBKQJK+RozRgDIpplwRPqaz2mw66J8ppg6owH3yrjQeg6YhJCTozdpqo2GAwGg8FgMBgMhr7BowfuOvvRg3Yb6Ukqyiq5Fd4QpFoYf7Hj6dntdv4ceUidcIFpQ/U5Ml4iFXtTMOsgZCDiFcsw95BqaCfl2ljcRjnkEx4kjiXzSl//n/rPAZh8mGljWGeGjS0ztWPVYALyPHbI7mO9SJRVcit+eDwHvBzizGAnC3ZC7Gm4tnC/YWN7UJy/NQAbI/bBIJHSvmJMKOxb8jV/fo5U35qhTMqF4XyBOLLv7O+noyqEtcWH+vSYtOLWYF/E9knIO34DfMe57coGDw2SUAM7HkSbmcmn1fyuOoaEW2e+C0luluq2wTJoKRcG8N9Lq2VECfRPaUm5sBPqPfX4HII6xYoZNIn397qyUxaTREKCHmH0XS+RxwtDwHqmOgcJ2WuYVmrT6Pv4Z+meUVQf68vuBWIpFy3jDdF3MI4f4D+jltYBjmj1fpEpm7zf9aaa0wANk/BD4uT9dFWxDxltIS0aERSISbOCysNnDNnwMePN+IDKu8snCEIID8SjejzWLXJSLoJl/kUdD4nzb7/DH1NvPRad4cuIgfABqdr7a9RxeqFhEuLbxWlPtNMF0hoTAtZjMKYjwhYv9C0qT09SGK/h+lrFH+NPRpESujwmNbdm7l9HykX5EH20SCeJmwQxZgxyLVyOM1Q5aA0fjspmyMAki4lLHGl2cDFJ3THghb4jk8cs+56Sa3m5YSfQ0SLNV3lE8SJcSGqMWFfKBXHp8gnICaEROAQlDKTcX30X3eTtUXn4tZGZrdbBvQcHk0hCFi+lhKuA2W0usBGzYUgWxAMQQq/toNt8sIt6lKFMyoXcjIlJ2Fs5DpjO56+qY4Yfqd1BDQ2RkH8+3RsvgMnIu8WN7cILocXB7IE9jdV2tCphpolsao44FQoLo5it0kqG7heZF+YeJF20bqhbPl/x5yiTcgFIOSRONkZLrlcEooqhS97GpwekpZShTIYQjGXDcMBUNAENGasZZ9HyMSFB1s6STm3Lw0bIbBFiYgZhZhzWHSPpYiCP0poZNZrEjaPyscvd57/zFakekLxMygX4kzHxYScBosnGBmdm73/2ScctpBWMpWg/r1hHg8FgMBgMBoPBYDD0HA/N2nL2w7O2GulFoqySW5VJuQChPG70Kd5gERceBmlMN0Qw0FHvUd5c5vMQPcyo9ktMAIblr2fyWPTEM8zK5DOrPic6h12TmT3+ZORgG/WgjtMDkOeRV20z1otEWSW3KpNyvVacWQOPA0IEbGw6qgJ2P/R9RD8YFuemC4AoH/DlonrBHrdFvV/lGdsjdViUyMMT8wNf32MT+Zie8OhoIzcmI553nrgAAOf7fItHAxokoUYs5WI7BS1KoBUKAZGw07EMNKhnlvHHW/nj+EVCkDriAOpGi80y0hQJ+aNAfMSzMQmpG646jPGahBjgH5KWK5FnQvwQ74s3mJgkEsZSLozRelsIuvXg/Mc7gXpF+4LxAR+VKZu8kzqsRwoYlukuUdPEJMR3jGAC702KhAw5GDYQlUGTELLxvMFfTdf8Nynf2mxw0DAJc1KuWD9IXBpaCmRPvNC7o3IgyYcT5eP6o5Ws2sLgpcElyGL9FAnnScs/HJOQKGKIE1aUiSQEPC/RIUbFeYVyMWwGDw2TMCflYpwEeYi0wKAe1TEkpNtivBjHG0QY8NHoHK0ML7ksyFA7KRfkD5s7xiSk9YNYwZ0Yk/BSde+YhFzDhIRNGZGGoY1k2GBjQjBJ3TGIpVxMBth1nckFs2MmBnTBTEYeiq5F3PpOdUzrw0ueL+Uok3KhWXzKl8EkivjUY/4zs3H+IKP+mMQ6EWbjSLZoefFbf9jn8Sf5q//M2hTIeZu0hhS08Aw3qi7Iml6YRBKWSbloDb/iPxM0CMlU6GJpPXSASl4orcolUq9lCYE6Q6LLH/OfWUdyUJSP1IygmIxlV4/yiDnzmP9Ma85E5drofrSUdYK8Tx80RMJ2Ui6UMugAkVDRajBof6m6nnJpXSAkM1/GbaFVQUsIIViIFLpavSNoVaTGhBqpiUlA3B1jE8V0BOnQH7JJIxMV2+UTNGSsbiflYiUdLw1iIs2PtXa8LMZRvEiM2doOyDrhWCY1V+qDP0iZ3ApbXyrWIiAIaLyLJ38uJi5IwTDYV12WajAYDAaDwWAwGAyGnuMPOwzN/uMO6430IlFWB7fEzEIwoTimcwDmG2RS8fYRYAOft586h22QBegYoFHXIIRYeuKlXSMl5SqLykV9WSjFjBqj+MslDWbOPIPtYRIAeR7Yaf2xXiTK6uCWBEjHWxDr7QD2PTbVwUj8rET+Ff5abeQmEj7uOkwee4kz9ZzY+S+QRErK1S4qFzsQEI4E89IccbbOF8t48OfAs4OpKrf1xOChYRJi30Mxg+84JmHQ6UGqFAkRkEJCXGWahHhKPqKOaRW/19WPMLEeKSlXu6hc2AC1GBc/sV6zDAhRwja6/FGMhAENkpCXOypO+RzcWxoIG2hZwkY6moTo9NiWATVyTEII+DVpeVDYwEaHaOsWOSlXu6hcCDC2V3kohn6ijtcU5+kh5o6RUKNBEqKqDrH+YhIGnR5jvhQJ6erCrk0xCfHN4jVB3s94jNb0eVINZVKudlG58NS8TX3/chkfT2dYWuNgI6FGQyRkjQjdVRCxxiREfRK6rpiE6PRQrARixSRkLMiYkIgHBMwkcsLBmXrUkXK1i8pFC8oOT/PEbSmGcjwEwiQMMX+QMNkyEmo0REJaALrTIIO62aegv6MrO93nEeAIEs4TN95CRPpTdS0C16ul1Q1SztvVvQ6WiSLYgDpSLlAWlQvMEBfWmD8JoU2u8+fpAb6lngGZFwKomfgAAAZ0SURBVIKMeKPJwURDJJwh46VOt/sUWrSzVB4vh5fPYihazv2ja2npaJFC/MI7xG39GkAr96cKP0U7KVeMOCpXDNbIhHUzx0RlMztm8mQiBjAJJhqQmpgEpMaEGnF3zEwbUqJ8hiy0dpd0WI8ypKRcZVG5XurPUwdUPARNykm1rDvWmARjNZgreblV2IwxR0K6yTnqGOM3AlheOmMyDMxrdFiPMqSkXGVRuTCgE9mV1hHZ2laSBzIv64oNBoPBYDAYDAaDoW9w3/orzb5vw1VGepKKsjq4ZU7KVRaVC3ce9kMkVNgUV5E0MGxf0OmzJ8Aa529EaQOVf0GU955MOamoXMzYETXgafmktBbRGyDP/RutOtaLRFkd3DIl5SqLysWmM/f465BQYY75fqLcV/hyc1tPdAKM3viwj1VJm3so++Mqb++4AElH5QKj4uyXuPBYLjpco57TCw2TMCflKovKBbnwwwa/K35obHV6D7ug02Ntc10SXl2ST9lla5pzUbnQFVLnsBPAev646i4D0wsNkrBMylUWlYvYLbHRGLfcYeqYbWnx5+Kyq0tCyIwYAh/08lE+ZbP4nlYwFqyCXFQunu3G6LuILmbUqOv0QYMkLJNylUXlovvC1xrGUAgP/iUt9Qo6PV44Y81ekJCF+fiN7/FlDal84skgG6O1xKV3vMori8qFm/E7Mh6Q/SAxNEbCdlKusqhctKA/EqekQUIFEWglQyBMwgcf4z+3I2E7KZde18IfAPfbWZl83HQEQVrRH5dF5UIPqZcCACJNHCaGxkjYTsoFclG5AOs7EI/S5bKOhO4Z/R6Thqf8NZR7cZEe959Ts89uN1hEqXNbJm8Jfy98wJ1E5fpxdD07WNnWYqAhEs6QcilXDB2VKwYCB0QCRMJaOSqXLu8J//m5Xf8YE0FrG0fTCqB8hgUvkfZRubb1x0GUwZ+HuDsWEAlMgokGxN1xu6hcCEfX9ecQil6UKbfumJAJBLZC5FjItmjNjvR5xMx+o7gxIsSD8PyRUlvaxt3xEv67BPZcydc/Dpg0uJgEYzWIpVztonLRvdHdQYpTJB/rmUnMNZ0+ewKYiSAxLRxLCpjthtV1kJBhAN0/kcWYoGyQKAOkonJxPQufiCzGqsG1atTTYDAYDAaDwWAwGAw9xZ0is+8SGelFutP5fHOYJ+NlUKervO1kooRqU5mIVFQugGmGmSf+2ToR8QnwjmEdE0pKElYm5dpPJj5DsFUulcg7VF2LLRFXJS5DjNpHyiAB8twjMtaLdJdbjZYDRl/0dEEGpaMkYFe7V8ZLqF6YKCMVlWtzcZ4LbHz8CVhxt1O1X+MZGRZEZjfOlL2xTMpFna6PniGYkjBS44Z8r8rbUV1LmDnMPzPEeYuwkw7O9hINk/B1mTxIeF0mLyAXlQvhw/nqGG/LF9uU1Q4YznMkzEm5qNN5mbxAwmUz+RBQt4y0uGdmvjv90DAJPyOuC2ORuPY0QEJ8y0ioeBkrRNeWReViEfrh6pgQbvE2ZN2ijIQ5KRd1wjfN89E160CdgYRH+xRHdEBpHrZJwzjO0OIYGRQ0SEKiVOFFoKvDv6vHXJCQ1gAJFQRCvLChyi+LyoWnRW9+vYM4z0od5EhYJuWiq4U8EAqtIG66IFqFhPf7PFp8PEN6OIIcjSizDAUItMmwpRfRZqcGGiShlkGhsEaEsHkijxcGWQNJ20XlWiDjN+hmkvN4pg7tpFwBORKWSbl0Hi05LXfQPC4R5aO81gE3aVkhHxMSoswSyUGPGac3GiRhDGaCb8/kobO7w3/uJCqXHk8xWbg3U26nUq4cCTW0lCsFutfLM3ks1KJ7XseX81hUzjyZKIKdvpgkEtICEF3rwEx+CIQO2kXlosU8Q10LUb8h9dAJCbWUKwXqdWEmbzN/Ld01Lf/CqBxaxh92Wtkpj4ZIiO6PH3YjnyATZoggPKXl20Vc1xtW3r1zYjHPIO6Od/HfZ7JDN894bN+KPwf2OrpqxnqP+M8hBHA7KRcTDjZSROCKHO1pf32oI5aBtX15jP20ymZU3BIHdIbMvol7+KGKzzD10JCxmvEcXSoSKF4O3aIe8zAOgjy0DgvEvYDUFhIgjsoFCHLOnsiE84VAKY1fJwgRwXS62Oe1k3IRJg4bJhMWxoN6XfUMccRiHPyoL1MvJWXZA3s40y0/6MvKmXMMBoPBYDAYDAaDwVAH/w9Fg4WALTLaJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "legend_url = sdatwms['1512_1'].styles['default']['legend']\n",
    "urllib.urlretrieve(legend_url, 'sample_data/avhrr_legend.png')\n",
    "display(Image('sample_data/avhrr_legend.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8WhcgbHEEUah"
   },
   "source": [
    "To retrieve the tree mortality rates, we use SDAT's WCS. Then we plots the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ynPVNxUnaXI6"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done!\n",
      "      mortality\n",
      "2003      10.25\n",
      "2004       0.00\n",
      "2005       2.25\n",
      "2006      97.25\n",
      "2007     246.00\n",
      "2008     453.75\n",
      "2009     984.75\n",
      "2010       0.00\n",
      "2011       4.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Tree Mortality (Mg C km-2)')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAD8CAYAAABKKbKtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9//HXZyYr2YZ9SQYBRWRNWKug1qWKWxWtVqwK2nrt7bW3tz97baW3621tvbW3v9af11bq2lYBrYpUrV7FWpcqECGAgJE9G0tYEkLIPp/fH3MCA2SZhJk5M8nn+XjkMTPfOXPOm5CZz5zzPef7FVXFGGOMiTSP2wGMMcb0TFZgjDHGRIUVGGOMMVFhBcYYY0xUWIExxhgTFVZgjDHGRIUVGGOMMVFhBcYYY0xUWIExxhgTFUluB4iWAQMG6IgRI9yOYYwxCeWjjz7ap6oDI7GuHltgRowYQWFhodsxjDEmoYjIzkityw6RGWOMiQorMMYYY6IiagVGRB4Xkb0i8nFIWz8ReUNENju3fUOeWyAiW0SkWERmh7RPFZH1znMPiohEK7MxxpjIiWYfzJPAQ8AfQtruBZar6v0icq/z+DsiMg6YC4wHhgFvisiZqtoC/Ba4E/gQeBW4DPhrdwI1NTVRVlZGfX19N/9JvVtaWhp5eXkkJye7HcUYkwCiVmBU9R0RGXFC8zXABc79p4C3ge847YtVtQHYLiJbgBkisgPIVtUPAETkD8AcullgysrKyMrKYsSIEdiOUNeoKvv376esrIyRI0e6HccYkwBifRbZYFXdBaCqu0RkkNOeS3APpVWZ09bk3D+xvU0icifBvR2GDx9+0vP19fVWXLpJROjfvz+VlZVuRzEGgKVrynng9WIqquoY5kvnntljmDO53Y8H44J46eRv6xNfO2hvk6ouVNVpqjpt4MC2T+O24tJ99rsz8WLpmnIWvLCe8qo6FCivqmPBC+tZuqbc7WgmRKwLzB4RGQrg3O512ssAf8hyeUCF057XRrsxphd74PVi6ppajmura2rhgdeLXUpk2hLrArMMmO/cnw+8FNI+V0RSRWQkMBpY6RxOqxGRs52zx+aFvKbXKyoq4tVXX+10ubfffpurrroKgGXLlnH//fcDsHTpUjZu3BjVjMZEQ0VVXZfajTuieZryIuADYIyIlInIV4D7gUtEZDNwifMYVd0APAtsBF4D7nLOIAP4GvAosAXYSjc7+Ltj6ZpyZt3/FiPvfYVZ978VV7vfzc3NYReYUFdffTX33nsvYAXGJK5hvvQutRt3RK3AqOpNqjpUVZNVNU9VH1PV/ap6saqOdm4PhCx/n6qerqpjVPWvIe2FqjrBee7rqtpuH0wkResY744dOzjrrLO44447mDBhAjfffDNvvvkms2bNYvTo0axcuZIDBw4wZ84cJk2axNlnn826desA+NGPfsSdd97JpZdeyrx58/jBD37AkiVLKCgoYMmSJaxcuZKZM2cyefJkZs6cSXHxyYcLnnzySb7+9a/zj3/8g2XLlnHPPfdQUFDA1q1bmTJlytHlNm/ezNSpU0/p32pMtHzrkjNPaktP9nLP7DEupDHt6bFjkXXmx3/ZwMaKQ+0+v6akisaWwHFtdU0tfPvP61i0sqTN14wbls0PPz++021v2bKF5557joULFzJ9+nSeeeYZ3nvvPZYtW8bPfvYz/H4/kydPZunSpbz11lvMmzePoqIiAD766CPee+890tPTefLJJyksLOShhx4C4NChQ7zzzjskJSXx5ptv8t3vfpfnn3++zQwzZ87k6quv5qqrruL6668HICcnh6KiIgoKCnjiiSe47bbbOv23GOOGvhkpAPRJ8XKksYUBmSl878pxdhZZnOm1BaYzJxaXztq7YuTIkUycOBGA8ePHc/HFFyMiTJw4kR07drBz586jheGiiy5i//79VFdXA8FDXOnpbR8GqK6uZv78+WzevBkRoampqUu57rjjDp544gl+9atfHd0jMiYeLV5VQv+MFP76zfOYcd9yvnzuSCsucajXFpjO9jRm3f8W5W10GOb60lny1XNOadupqalH73s8nqOPPR4Pzc3NJCWd/N/SeopwRkZGu+v9/ve/z4UXXsiLL77Ijh07uOCCC7qU6wtf+AI//vGPueiii5g6dSr9+/fv0uuNiYW9NfUs37SXL587kkFZaYzo34e1pVVuxzJtiJfrYOLOPbPHkJ7sPa4tVsd4zz//fJ5++mkgeAbYgAEDyM7OPmm5rKwsampqjj6urq4mNzf4Le7JJ5/sdDsnvj4tLY3Zs2fzta99jdtvv/0U/xXGRMfzH5XTHFBunB68sqHA72NtabXLqUxbrMC0Y87kXH5+3URyfekIwT2Xn183MSa74T/60Y8oLCxk0qRJ3HvvvTz11FNtLnfhhReycePGo5383/72t1mwYAGzZs2ipaWlzdeEmjt3Lg888ACTJ09m69atANx8882ICJdeemlE/03GRIKqsmRVCTNG9OP0gZkA5Pt97D5Uz+5qG2Mw3kiMTsqKuWnTpumJE45t2rSJsWPHupQoMfzyl7+kurqan/zkJ20+b79D46YPtu7npt9/yK++mM91U4LXYK8uOch1D/+DR26dyuzxQ1xOmPhE5CNVnRaJdfXaPhhzsmuvvZatW7fy1ltvuR3FmDYtXlVCVloSl08YerRt3NBskjxCUWmVFZg4YwXGHPXiiy+6HcGYdlUdaeSvH+9m7nQ/6SnH+kfTkr2MHZptHf1xqNf1wfTUQ4KxYL8746YX15TT2Bw42rkfqsDvY11ZNYGA/Y3Gk15VYNLS0ti/f799UHZD63wwaWlpbkcxvZCqsnhlKZPychg/LOek5/P9Pg43NLNt32EX0pn29KpDZHl5eZSVldmcJt3UOqOlMbFWVFpF8Z4a7rt2QpvPF/hznOWqOWNQViyjmQ70qgKTnJxsszEak4CWrColPdnL1fnD2nx+1IBMMlOTKCo9yPVT7UtQvOhVh8iMMYnncEMzy9ZWcNWkoWSlJbe5jMcjTMrLsQsu44wVGGNMXHt5bQVHGluYO+PkadBDFfh9bNp1iPqmzi8yNrFhBcYYE9cWrSpl9KBMpgz3dbhcvt9Hc0DZuKv9UdJNbFmBMcbErU27DrG2tIq5M4YfHfC1PQX+YAGy62HihxUYY0zcWrKqlBSvh2vDGANwcHYaQ7LTKLICEzeswBhj4lJ9UwsvrC5j9oQh9HMmGOtMvj/H9mDiiBUYY0xceu3j3Ryqb2ZuG1fut6fA35cd+49QdaQxislMuKzAGGPi0uJVJQzv14dzRoU/8V2+c8Hl2jI7XTkeWIExxsSd7ftq+XDbAW6c7sfj6bhzP9TE3BxErKM/XliBMcbEnSWrSvF6pMtX5WelJXPGwEzr6I8TVmCMMXGlqSXAnz8q48Ixgxic3fXBVfP9PtaWVtmgtnHACowxJq4s37SXfYcbuGlG+J37oQr8PvbXNlJ2sC7CyUxXWYExxsSVxatKGJKdxmfPHNit1x+94LLMDpO5zQqMMSZuVFTV8fdPK7lhWh5J3u59PI0ZkkVKksc6+uOAFRhjTNx4trAUgC9O697hMYBkr4cJw7JtZOU4YAXGGBMXWgLKs6tKOfeMAfj79TmldeX7fawvr6a5JRChdKY7rMAYY+LCu5srqaiuZ+70joflD0eB30ddUwuf7rEplN1kBcYYExcWryylX0YKnxs36JTXZR398aHTAiNBE0VktoicLyLhj9tgjDFhqKxp4M1Ne/jClFxSk7ynvL7h/frg65NsHf0ua7fAiMgIEXkY2Ar8GrgduBt4R0TeF5FbpbMJGtpf9/8RkQ0i8rGILBKRNBHpJyJviMhm57ZvyPILRGSLiBSLyOzubNMYE7+eX11Gc0C5sQsDW3ZERMjP89kV/S7raA/mF8BzwBmqerGqzlXVOao6HrgeGAzM7+oGRSQX+AYwTVUnAF5gLnAvsFxVRwPLnceIyDjn+fHAZcDDInLqX3GMMXFBVVmyqpTpI/pyxqCsiK033+/j0z011DY0R2ydpmvaLTCq+kVV/ZuqnnQahqruUtVfquqT3dxuEpAuIklAH6ACuAZ4ynn+KWCOc/8aYLGqNqjqdmALMKOb2zXGxJkV2w+wfV8tN0agcz/UZL+PgMLH5Xa6sls67IMRkQwROa2N9vHd3aCqlgO/BEqAXUC1qv4vMFhVdznL7AJae/pygdKQVZQ5bW3lvVNECkWksLKysrsRjTExtGRVKVmpSVw5cWhE1zspr3XofjtM5paO+mC+QHBv4RURWS8iU0Ke/mN3N+j0rVwDjASGARkicktHL2mjrc1R7FR1oapOU9VpAwd2b5gJY0zsVB9p4tX1u7hm8jDSUyJ75Lt/Zir+ful2waWLOtqD+T7H+km+CiwSkaud57rVue/4HLBdVStVtQl4AZgJ7BGRoQDO7V5n+TIgtOcvj+AhNWNMgltaVE5DcyAi1760xTr63dVRgfE4h7NQ1X8AFwE/FpG7aGcPIkwlwNki0sc5C+1iYBOwjGMnDcwHXnLuLwPmikiqiIwERgMrT2H7xpg4oKosWlnChNxsJuTmRGUbBX4f5VV1VNY0RGX9pmMdFZha5wMdONp3cgFwAzCuuxtU1RXAn4HVwHonw0LgfuASEdkMXOI8RlU3AM8CG4HXgLtUtaW72zfGxId1ZdV8srsmansvEHLBpe3FuCKpg+fuOvF5Va0WkUuBm05lo6r6Q+CHJzQ3ENybaWv5+4D7TmWbxpj4snhVCenJXq4uGBa1bYwfloPXI6wtq+Jz4wZHbTumbR2dprxaVTeHtonIZaraqKpPtfc6Y4zpTG1DM8uKKrhy0lCy05Kjtp30FC9jBmdZP4xLujoW2c+iksIY06u8vK6C2sYW5kboyv2O2BTK7ulqgTmVs8eMMQaAxatKOWNQJlNP69v5wqeowJ/Dofpmtu+rjfq2zPG6WmD+JSopjDG9RvHuGtaUVDF3up9uDmfYJQX+YBGzCy5jr6NOfgBExENwDLARQJKITAdQ1QejG80Y0xMtXlVCsle4bkpeTLZ3xqBM+qR4WVtazbWTY7NNE9RpgSF4PYoSPKXYpoczxnRbfVMLL64p59LxQ+iXkRKTbXo9wsTcHOvod0E4BWaEqk6MehJjTI/3+obdVB1p4qYoXvvSlgK/jyfe30Fjc4CUJJtnMVbC+U2/LiIXRT2JMabHW7yyFH+/dGaeHtt5C/P9PhpbAmzadSim2+3twikw7wJ/EZHDInJARA6KyIFoBzPG9Cw79tXywbb93DjNj8cT2xNS820KZVeEU2D+L3Ae0BcYCAxwbo0xJmxLCkvxCFw/NfrXvpxoWE4aA7NSrR8mxsLpg9kMrFG7SskY001NLQH+/FEZF501iCE5aTHffusUyjYmWWyFU2AqgLdE5FWC44UBdpqyMSZ8b32yl8qahojPWtkVBf4c3ty0h0P1TVEdnsYcE84hsjLgPSCb4KGx1sNkxhgTliWrShmUlcqFY9w7ut7aD7POJiCLmXD2YH6vqiWhDSfMbmmMMe3aVV3H28V7+doFp5Pkde8U4Ul5xzr6zx1t35FjIZz/7RdaZ5oEEJFZwB+iF8kY05M8V1hGQOHGae4dHgPISU9m1MAM6+iPoXAKzF3ASyIySERmAw8DV0Y3ljGmJwgElCWrSpl1Rn+G9+/jdhwKnCmU7Zyl2Oi0wDgzUN4NvAH8FLhUVXdGO5gxJvG9t2Uf5VV1UZ21sivy/T4qaxrYfaje7Si9Qrt9MCLyIsExyFr1AaqA34oIqnpdtMMZYxLb4lUl9O2TzKXj42M2ydaO/qKSKoZOTHc5Tc/XUSf/QzFLYYzpcfYdbuCNjXuYd84IUpO8bscBYOzQLFK8HorKqrh84tDOX2BOSbsFRlWXxzKIMaZneWF1GU0tGpNZK8OVmuRl7LBsu+AyRmxYUWNMxKkqi1eVMvW0vowenOV2nOMU5OWwvqyaloB19EebFRhjTMSt2nGQbZW13BhHey+t8v0+ahtb2Fp52O0oPV67BUZEBojImDbazxKR2I61bYxJKItXlZCZmsRVk+Kvn+NoR78dJou6jvZgHgSGtdE+EvhNdOIYYxJddV0Tr67fxdUFw+iTEs5gIbE1sn8GWWlJVmBioKMCk6+qfzuxUVX/ChREL5IxJpEtKyqnvikQ81krw+XxCAV+G1k5FjoqMB199bChSI0xJ1FVFq0sZfywbCbm5bgdp135eT4+2V1DfVOL21F6tI4KzFZnaJjjiMilwPboRTLGJKqPyw+xcdehuDo1uS35fh8tAWVDhY2sHE0d7aXcTXCq5L8DHzlt04Dzgc9HO5gxJvEsWlVCWrKHqwty3Y7SoXxn76qotJqpp/VzOU3P1e4ejKp+AkwEVgBnOT8rgEnOc8YYc9SRxmaWFVVwxcSh5KTH91H0QdlpDMtJs47+KOvwFA9VrQd+H6MsxpgE9vK6XRxuaOamGfHZuX+iguHW0R9tdqGlMSYiFq8s4fSBGUw7ra/bUcKSn+ej5MARDtQ2uh2lx3KlwIiIT0T+LCKfiMgmETlHRPqJyBsistm57Ruy/AIR2SIixW2deGCMcdene2pYXVLF3OnDERG344Sl9YLLtWW2FxMtbu3B/AZ4TVXPAvKBTcC9wHJVHQ0sdx4jIuOAucB44DLgYRGJj6FZjTEALF5ZSrJXuG5KfHfuh5qYm4NHsMNkUdTpZbYisobj54UBqAYKgZ+r6oGubFBEsgmeiXYbgKo2Ao0icg1wgbPYU8DbwHeAa4DFqtoAbBeRLcAM4IOubNcYEx0NzS28sKaMS8cNoX9mqttxwpaRmsToQVnW0R9F4ezBvEFwj+Irzs8bBD/cDwJPdmObo4BK4AkRWSMij4pIBjBYVXcBOLeDnOVzgdKQ15c5bScRkTtFpFBECisrK7sRzRjTVa9v2EPVkaa4HNiyM61X9NsUytERToGZqar3qOoa5+c7wPmq+jOC45J1VRIwBfitqk4GanEOh7WjrQO6bf41qOpCVZ2mqtMGDhzYjWjGmK5asqqEXF86554xwO0oXZbv93HwSBOlB+rcjtIjhVNgskRkausDEZkCZDsPm7uxzTKgTFVXOI//TLDg7BGRoc42hgJ7Q5YP/WqUB1R0Y7vGmAjbub+W97fs58bpfjyexOjcD5Xvdy64tI7+qAinwHwV+KNzdtcW4E/AV53DWr/o6gZVdTdQGjIVwMXARmAZMN9pmw+85NxfBswVkVQRGQmMBlZ2dbvGmMh7trAUj8AN0/LcjtItZw7OIi3ZYx39UdJpJ7+qfgiMc+aAEVXdF/L0om5u91+Bp0UkBdgG3E6w2D0rIl8BSoAbnO1vEJFnCRahZuAuVbUR6oxxWXNLgOcKy7hgzCCG5qS7Hadbkr0eJgzLsY7+KGm3wIjIN9ppB0BVH+zuRlW1iOC4Zie6uJ3l7wPu6+72jDGR97fiSvbWNMT9wJadKfD7+OOHO2lqCZDstWvPI6mj3+avCZ5KPBQYAAw84ccY04stXlnCwKxULjxrUOcLx7F8v4+G5gDFu2vcjtLjdHSI7DPAjcAVwIfAIlV9OxahjDHxbXd1PX8r3ss/f/b0hP/WXxByRf+E3PidwyYRdTSa8ipV/XdVzQeeBr4oIhtE5KrYxTPGxKPnCksJKAl57cuJ8vqm0y8jxTr6o6DTrx4i0g8YC4wB9gBdunLfGNOzBALKksJSZp7en9P6Z7gd55SJCPl51tEfDe0WGBG5VUReIXi6cDpwi6pepKr/iFk6Y0zceX/rPsoO1vWIvZdWBf6+bN57mMMN3bm0z7Snoz6Yp4D1BE8j/ixwfugoqap6XXSjGWPi0eJVpfj6JDN7/BC3o0RMvj8HVVhfVs05p/d3O06P0VGBuSRmKYwxCWH/4Qb+d8Nubjn7NNKSe86g5vl5xzr6rcBETrsFRlWXxzKIMSb+vbimnKYWZe70xJi1Mlx9M1I4rX8f6+iPsI76YJaKyOUiclIREpHTROQHIvLl6MYzxsQLVWXRyhImD/cxZkiW23EiLj/PZx39EdbRWWR3ETxM9qmIfCAiy0Tkf53xyJ4ANqjq4zFJaYxx3Uc7D7K1spabetjeS6sCv49d1fXsOVTvdpQeo6NDZOXA3cDdInIGwSv664BiVbVLXo3pZRatLCUjxcuVk4a6HSUqjk6hXFrFpT3oBAY3hXUJrqpuUdV3VbXQiosxvc+h+iZeWV/B1QW5ZKR2OkZuQho/LJskj7DWhu6PmMQe48EYExMvFVVQ3xRI+IEtO5KW7OWsoVmsLa12O0qPYQXGGNOpJatKGDs0m0l5PXusrvy84BTKgYBNoRwJ4QwVc5mEXmFpjOlVPi6v5uPyQ8yd7qenfxQU+H3UNDSzbV+t21F6hHD2YG4DNovIz0RkdJTzGGPizOJVJaQmeZhTkOt2lKgrCOnoN6eu0wKjqnMJTg5WDiwSkXdF5MvOlMnGmB7sSGMzL62p4IqJQ8npk+x2nKgbNTCTzNQk6+iPkHDPIqsCngGeBIYDNwFrReRfohfNGOO2V9fvpqahuUd37ofyeoSJuTm2BxMh4fTBXC4izwHvAlnA2ap6CZAPfCfK+YwxLlq8soRRAzKYMbKf21FiJt/vY+OuQ9Q3tbgdJeGFswdzK/BbVZ2gqj9X1V0AqloL/FNU0xljXLNlbw2FOw9yYy/o3A9V4PfR1KJs2nXI7SgJL5wCs0NV3wptEJGfAajq/0YllTHGdYtXlpLkEa6bkud2lJiyjv7ICafAXNZG25WRDmKMiQ9L15Qz8+fLefS97SR5hfe37HM7UkwNyUljcHYqa8vsgstT1e6YDyLyVeCfgTEisjrkqSygMNrBjDGxt3RNOQteWE+d0/9Q3xRgwQvrAZgzueefptyq9YJLc2o62oN5FrgBeMW5bf2Zpao3xSCbMSbGHni9+GhxaVXX1MIDrxe7lMgd+X4f2/bVUn2kye0oCa2jAtOoqluArwCVIT/1IpIdi3DGmNiqqKrrUntPNdl/bIZL030dFZg/O7cbgI+d2w0hj40xPciRxmaSk9r+SBjmS49xGndNyMtBxDr6T1VH88Fc7tz2jiusjOnFjjQ2c/sTq2hsDpDsFZpajg32mJ7s5Z7ZY1xMF3vZacmcPjDT9mBOUUed/JM6eqGqrot8HGNMrNU2NHP7k6so3HGA38wtQDXYF1NRVccwXzr3zB7Tqzr4W+Xn+fj7p5Woaq+6DiiSOpo56H86eE6B8yOcxRgTY7UNwT2Xwp0H+PXcyVydPwzoXWeMtafAn8Pzq8uoqK4nt5cdIoyUjg6RnRfLIMaY2GotLh+VHOQ3cyfzeae4mKACf18AikqqrMB0U1hzn4rIWcA4IK21TVWfiVYoY0x01TY0c9sTK1ldUsVv5hZw1SQrLicaMySLlCQPa8uquHLSULfjJKRwBrv8HrAQ+B1wOfBr4PpT3bCIeEVkjYi87DzuJyJviMhm57ZvyLILRGSLiBSLyOxT3bYxvdlhKy5hSUnyMH5YNkV2Jlm3hTNUzI3AhcAuVb2V4CjKYe35dOLfgE0hj+8FlqvqaGC58xgRGQfMBcYTHLbmYRHxRmD7xvQ6hxuaue3xYHF5cO5kKy6dyM/zsb6smuaWgNtRElI4BaZOVVuAZhHJAnYDo05loyKSR3A8s0dDmq8BnnLuPwXMCWlfrKoNqrod2ALMOJXtG9MbtRaXNaXB4mKHfTpX4PdR19TClsrDbkdJSOEUmDUi4gMeJzgG2Upgdccv6dSvgW8DoV8LBodMBbALGOS05wKlIcuVOW3GmDDV1Dcx3yku/+8mKy7hah1ZuajEDpN1RzhTJn9VVatU9X8I7nV8VVXndXeDInIVsFdVPwr3JW3Famfdd4pIoYgUVlZWdjeiMT1Ka3FZW1rFQzdN5oqJVlzCdVr/PuSkJ9sFl90UTif/0TlfVHWLqq4ObeuGWcDVIrIDWAxcJCJ/AvaIyFBnm0OBvc7yZUDoaAJ5QEVbK1bVhao6TVWnDRw48BQiGtMztBaXdWXVPPSlyVxuxaVLRIR8v4+iUhu6vzvaLTAikuIMajlYRLJEJNv5yQOGd3eDqrpAVfNUdQTBzvu3VPUWYBkw31lsPvCSc38ZMFdEUkVkJDCa4GE6Y0wHDtU3Me9ocZnCZROsuHRHQV4On+6p4Uhjs9tREk5HZ4PdBdxNsC9kA8cOVR0ieMpypN0PPCsiXwFKCE4NgKpuEJFngY1AM3CXc9KBMaYdh+qbmPfYSj4uby0uQ9yOlLDy/T5aAsqGikNMH9HP7TgJpaMr+f+viDwI3KOq90dj46r6NvC2c38/cHE7y90H3BeNDMb0NKHF5X9unsLs8VZcTkV+SEe/FZiu6bAPxtlTmNPRMsaY+FFd18Stj61kQ0U1D1txiYgBmank9U2nyDr6uyyc05TfEJFrop7EGHNKquuCfS4bK6p5+OapXGrFJWLy/TaFcneEU2C+DrwoInUickBEDorIgWgHM8aEr7quiXmPrThaXC4ZN9jtSD1KQZ6PsoN17Dvc4HaUhBJOgRkAJAOZwEDnsZ0DbEycOFpcdh3it1ZcoqK1H2adHSbrknAutGwBZhPsZL8PuMTO4jImPlQfaeLWkOLyOSsuUTEhNxuvR+yK/i4K50LL+wgO67LN+fm2iPw02sGMMR2rPtLELY+t4JNdNfzuFisu0dQnJYkzB2dRVGYXXHZFOKMifx6Y3LrXIiKPExyL7HvRDGaMaV9rcSneXcPvbp3CRWdZcYm2An8Or67fbVMod0E4fTAA2SH3s6IRxBgTnqojjdz82IcU767hkVunWnGJkfw8H9V1Tezcf8TtKAkjnD2YXwCrRWQ5wav5LwB+EM1Qxpi2VR1p5OZHV7B5z2EeuXUqF541qPMXmYho7ehfW1bFiAEZLqdJDOF08v8JOBd41fk5X1WfjnYwY8zxjhaXvYd5ZJ4Vl1g7c3AWfVK8rLGO/rC1uwcjIpNOaNri3PYXkf6qui56sYwxoQ7WBovLlsrDLLx1KheMseISa16PMCE3x4bu74KODpEVERxgsvWiytBeLQXOj1YoY8wxocXl9/Om8dkz7TI0txT4fTz5jx00Nge3xtPUAAAVhklEQVRISQq3C7v36qjAfAe4FqgCFgEvqar1bhkTQwec4rLViktcyM/z0dgcoHh3DRPzctyOE/faLcGq+oCqzgS+RXAOlr+LyDNtHDozxkRBa3HZVnmYR624xIV8f7Co2MCX4Qmnk38zsAR4mWBn/7hohzKmtztQ28iXfv8h25w9l/OtuMSFXF86AzJT7Yr+MHXUyT+c4IyT1wK7CRaZs+wwmTHR1Vpctu+r5dH50zhvtBWXeCEiFPitoz9cHfXB7ADWAS8S7IcZBNzRegWrqj4Y7XDG9Db7Dzdw86Mr2L6vlsfmT+fc0QPcjmROkJ/nY/knezlU30R2WrLbceJaRwXmZwTPFksiOIKyMSaKrLgkhny/D1X4uKyamWfY/1FHOpoy2cYaMyZG9h1u4Obfr2DngVoev206s+yDK27l5zlTKJdVWYHpRDhDxRhjomjf4Qa+9PsPKTlwhMfnT7cPrTiX0yeZUQMyrKM/DHalkDEusuKSmPL9PuvoD4MVGGNcUlnTwE0LneJymxWXRJKfl8OeQw3srq53O0pcC2fCsYEi8oiIvOw8Hicit0U9mTE9WGVNcM+l7GAdT9w2g5mnW3FJJK0jKxeV2l5MR8Lpg3kSeJrg0DEArRdePhmdSMb0TEvXlPPA68VUVNXh9QgC/PGOz3D2qP5uRzNdNHZoNsleYW1ZFZdNGOJ2nLgVziGyQar6DBAAUNUmoCWqqYzpYZauKWfBC+spr6pDgeaAIh6xQywJKi3Zy7ih2ay1PZgOhVNgakWkH8FrYhCR6UBNVFMZ08M88HoxdU3Hfy9rbA7wwOvFLiUypyrf72NdWTUtAXU7StwKp8D8O/AXYJSI/J3gyMr/GtVUxvQgjc0Byqvq2nyuop12E//y83wcbmhmW+Vht6PErU77YFS1UEQuBMYSnBNmo6o2Rj2ZMT3A1srD/NviNe0+P8yXHsM0JpJCO/pHD85yOU18CucssnTgbuBrqloEDBeRy6OezJgEpqosWlnCVQ++R/nBOr48awTpyd7jlklP9nLP7DEuJTSnatSADLLSkux6mA6EcxbZ48B6gkP1A1QAzwF/jVYoYxLZwdpG7n1hHa9v2MO5Zwzgv7+Yz+DsNCbl+Y6eRTbMl849s8cwZ3Ku23FNN3k8Qn6ej7Wl1W5HiVvhFJjRqnqTiNwAoKpHpHVIZWPMcd7fso+7ny3iQG0j/3HFWL5y7kg8nuDbZc7kXCsoPUy+P4dH/r6N+qYW0k7YQzXhdfI3ikgax84iGwl0uw9GRPwi8jcR2SQiG0Tk35z2fiLyhohsdm77hrxmgYhsEZFiEZnd3W0bEy2NzQF+/uombnlsBZmpSbz4L7P4p/NHHS0upmfKz/PRHFA2VBxyO0pcCqfA/CfwGpAnIk8BfwMWnMI2m4FvqepY4GzgLhEZB9wLLFfV0cBy5zHOc3OB8cBlwMMiYl8VTNzYsvcw1z78Po+8s40vzRjOy/96HhNybb723qDA6ei362Ha1uEhMudQ2FrgBmAmwbPI7lHVvd3doKruAnY592tEZBOQC1wDXOAs9hTwNsHRA64BFqtqA7BdRLYAM4APupvBmEhQVZ5ZWcJPXt5IerKXhbdO5dLxdlV3bzIoO42hOWnW0d+ODguMqqqIvKyqU4GXIr1xERkBTAZWAIOd4oOq7hKRQc5iucCHIS8rc9qMcc2B2ka+8/w63ti4h/NGD+C/b8hnUHaa27GMCwr8PtuDaUc4h8hWisiUSG9YRDKB54FvqmpHBzDbOojd5qWzInKniBSKSGFlZWUkYhpzkvc27+OyX7/D34sr+d6VY3nq9hlWXHqxfL+PHfuPcLDWLg88UbsFRkRa927OJVhkikVktYisEZHVp7JREUkmWFyeVtUXnOY9IjLUeX4o0HoYrgzwh7w8j+Cp0idR1YWqOk1Vpw0cOPBUIhpzkobmFu57ZSO3PLaC7PRkXrxrJnecZx35vV3rDJd2mOxkHR0iWwlMAeZEcoNOv85jwCZV/VXIU8uA+cD9zu1LIe3PiMivgGHAaCebMTGzZW8N31hUxMZdh7jl7OH8xxXjSE+xc00MTMzLQQTWllZzwZhBnb+gF+mowAiAqm6N8DZnAbcC60WkyGn7LsHC8qyIfAUoIXhiAaq6QUSeBTYSPAPtLlW10ZxNTKgqT68o4aevbKRPShKPzpvG58YNdjuWiSOZqUmMHpRpezBt6KjADBSRu9t78oS9j7Cp6nu03a8CcHE7r7kPuK872zOmu6wj34SrwO9j+aa9qCp2HfoxHRUYL5BJ+8XAmB7r3c2V3P3sWqqPNPH9q8Zx+8wR1tdi2pXv9/FsYRllB+vw9+vjdpy40VGB2aWq/xmzJMbEgYbmFh54rZhH39vO6EGZPHX7DMYNy3Y7lolzrR39RaVVVmBCdNoHY0xvsXlPDd9YXMSmXYeYd85pfPeKsTa+lAnLmCFZpCZ5WFtaxefzh7kdJ250VGDa7A8xpqdRVf60ooSfvryRzNQkHps/jYvHWke+CV+y18OE3Bzr6D9BuwVGVQ/EMogxbth/uIHvPL+ONzft5bNnDuSBGyYxKMs68k3XFfh9PL1iJ80tAZK84VzD3vPZb8H0Wu98Wsllv3mXdz7dxw+uGscTt0234mK6Ld/vo74pQPGeGrejxI1w5oMxpkepb2rhF68V8/j72zlzcCZ/+PIMxg61jnxzagpar+gvrWb8MBtNG6zAmF7m0z01fGPRGj7ZXcP8c05jgXXkmwjx90unb59k1pZW8aXPDHc7TlywAmN6BVXljx/u5L5XNpGVlsQTt03nwrNsWA8TOSJCvt9nHf0hrMCYHm/f4Qa+8+d1LP9kLxeMGcgD1+czMCvV7VimByrw+3jw083UNjSTkWofr/YbMD3a28V7+ffn1nGovokffn4ct80cYUN5mKjJ9/sIKKwvr+bsUf3djuM6KzCmR6pvauG/XvuEJ97fwZjBWfzpjhmcNcQ68k10HR26v7TKCgxWYEwPsHRNOQ+8XkxFVR3DfOnccvZwXiqq4JPdNdw2cwT3Xn6WdeSbmOiXkcLwfn2sH8ZhBcYktKVrylnwwnrqmoIzOJRX1fFfrxWTmerlidunc6HNz2FiLN/vY/XOg27HiAt2oaVJaA+8Xny0uITKTE224mJcUeD3UV5Vx96aerejuM4KjEloFVV1bbbvOWRvbuOOAn/wIsu1pdUuJ3GfFRiTcFSVdzdXcutjK9B2lhnmS49pJmNajR+Wg9cjrC21fhjrgzEJo6klwCvrdrHwnW1s3HWIgVmpXDlxCMs/2Ut9U+DocunJXu6ZPcbFpKY3S0v2ctaQLOvoxwqMSQA19U0sWVXK4+9tp6K6njMGZfKLL0zimsnDSE3ynnQW2T2zxzBncq7bsU0vlu/38fLaCgIB7dUzoVqBMXFrd3U9T/xjO8+sKKGmvpnPjOzHT6+dwAVnDjruTTtncq4VFBNXCvw+nllRwo79tYwamOl2HNdYgTFxp3h3DQvf2cayteW0BJTLJw7lzvNGke/3uR3NmLAUOH+ra8uqrMAY4zZV5YOt+3nknW38/dNK0pO93PyZ0/jyrJEM729znJvEcvrATDJSvBSVVHHt5Dy347jGCoxxVXNLgFc/3s3Cd7bycfkhBmSm8K1LzuSWs0+jb0aK2/GM6RavR5iYl0NRWe8+VdkKjHFFbUMzS1aV8th72ymvqmPUgAx+ft1Erp2ca8O6mB4h3+/jifd20NDcQmpS7/ybtgJjYmrvoXqe/McO/vThTg7VNzN9RF9++PlxfG7s4F59to3peSb7fTzSEuCTXTW9tv/QCoyJiS17gx33S9dU0BQIcNn4IfzT+aOYMryv29GMiYr8kI5+KzDGRJiqsmL7AX7/zjaWf7KX1CQPX5yexx3njmLEgAy34xkTVUOy0xiUlUpRSRXzznE7jTuswJiIa24J8NqG3fz+nW2sLaumX0YK3/zcaG49+zT6Z9pMkqZ3aJ1CuagXX9FvBcZEzJHGZp4rLOPR97ZReqCOEf378JM5E7h+Sh7pKb2zk9P0bgV+H29s3EN1XRM56clux4k5KzDmlFXWNPCHD3bwxw93UnWkicnDffzHFWO5ZNwQvNZxb3qx1hku15dVc+7oAS6nib0eW2DWl1cz6/63bFyqKNpaeZhH393O86vLaGoJ8Lmxg/nq+aOYNqKf29GMiQuTWofuL6uyAhPPROQy4DeAF3hUVe/v7DXlVXUseGE9gBWZbmhrEMlrCoZRuPMgC9/Zxpub9pDs9fCFKXnccd5ITu/FQ2IY05bstGROH5jBmpL474dpfb+nDDljaqTWKartzagRP0TEC3wKXAKUAauAm1R1Y3uvSR06WofO/3XwfpKHz40bTGZKEn1SvWSmJtEnJYnMVC8ZR+8fey4jNYmMlOBzyd7oTJkT7yMAnzgVMUCKVxiak8bOA3X4+iQz7+zTuPWcEQzMso57Y9pz/W/fZ3VJFarE5Xsdjn+/73rqmzTs2hyRY9uJsgczA9iiqtsARGQxcA3QboEJ1dAcYNOuQ9Q2NFPb0EJtYzPh1tUUr4cMpxBlpCSdcP/Y42DROn65zDaKV2qSh5eKKk6aR747e1qBgNKiSktACagSUIL3nfZAwGlz7h9bTmkJcMLjY+0BVX76ysaTpiJubFHKqur58dXjuWFaHn1SEuXPxxh3LF1TztqyagLO501H73VVpalFaQ4EgrctAZoDSlNLgObj2oP3O3qu9fVNAaXl6HqOtZ247udXl7U59fipSpRPiFygNORxGfCZsF/sS+etb11w9LGqUtfUEiw2Dc0cbmjmSGPwfm1js9PWwpGGZg43NnPkhOVq6pvZc6j+aLGqbWimqSW8iuX1CIGAnjQTY11TC996di33//WT4z70W5wicXwR0KN/sLEWCCjzZ45wZ+PGJJgHXi8+6bOhrqmFu58t4ofLNhz3gR/L93SSR0jyCskeD0le4Uhj5IsLJE6BaWt37aT/DhG5E7gTIGXIGUDbsxuKCH1SgnsXkTq809DcwpGGlqNF6HBDsPAcaXSKVWPz0bb/+dvWNtfRosr5Zw7A6xFEBK8IXo/gEcEjweLk8QTbPR6nzbnvdR57nNcce53g9Rxr98ixdXg9HL+d1nULfGPxGvYdbjwpo01FbEz4Kqrq2mwPKMwpGEaS13P0g97rEZK9EmzzCMkhzyU57cmeY88neYUkj6eT1zjrDm1zPl9Czbr/LcrbyXoqEqXAlAH+kMd5QMWJC6nqQmAhBPtgcmN4vDM1yUtqkjesEYCXrqlo8z8z15fOL67Pj0a8LvveleNO6oOxqYiN6ZphvvR23+s/vmaCC4nads/sMSe93yMhOj3YkbcKGC0iI0UkBZgLLOvoBRNzc3j/3ovirjMNgv+Z6SeMGBxvH95zJufy8+smkutLRwi+IX5+3cS4/H0aE68S4b0Ox7/fIykhziIDEJErgF8TPE35cVW9r6Plp02bpoWFhTHJ1h3xfhaZMSYyEu29LiIfqeq0iKwrUQpMV8V7gTHGmHgUyQKTKIfIjDHGJBgrMMYYY6LCCowxxpiosAJjjDEmKqzAGGOMiYoeexaZiNQAxW7nCMMAYJ/bITqRCBnBckaa5YysRMk5RlWzIrGiRLmSvzuKI3WqXTSJSGG850yEjGA5I81yRlYi5YzUuuwQmTHGmKiwAmOMMSYqenKBWeh2gDAlQs5EyAiWM9IsZ2T1upw9tpPfGGOMu3ryHowxxhgXJUyBERG/iPxNRDaJyAYR+TenvZ+IvCEim53bviGvWSAiW0SkWERmh7S/JiJrnfX8TkS8bW3T7Zwhzy8TkY/jMaOIvO20FTk/g+I0Z4qILBSRT0XkExH5QrzlFJGskN9jkYjsE5Ffx1tOp/0mEVkvIuuc99OAOM15o5Nxg4j8IlIZu5NTRPo7yx8WkYdOWNdU5/e5RUQeFJG2Jll0O+N9IlIqIofDDqCqCfEDDAWmOPezgE+BccAvgHud9nuB/3LujwPWAqnASGAr4HWey3ZuBXgemBuPOZ3nrwOeAT6Ox4zA28C0BPg//zHwU+e+BxgQjzlPWO9HwPnxlpPg5Q17W3+Hzut/FIc5+wMlwEBnuaeAi13MmQGcC/wz8NAJ61oJnEPwM+mvwOVxmPFsZ32Hw95+pH7Zsf4BXgIuIXgx5dCQX2axc38BsCBk+deBc05YRzLwF+DGeMwJZALvOX8QESswEc74NlEqMBHOWQpkxHvOkLbRTmaJt5zO+6YSOI3gB+LvgDvjMOd04M2Q9luBh93KGbLcbYR8eDvLfBLy+CbgkXjKeMJzYReYhDlEFkpERgCTgRXAYFXdBeDcth6iySX4Bm1V5rS1ruN1gt/CaoA/x2nOnwD/DRyJRr4IZQR4wjmk8/1I7dpHMqeI+JzHPxGR1SLynIgMjrecJ6zqJmCJOu/oeMqpqk3A14D1BKcuHwc8Fm85gS3AWSIyQkSSgDkcP/V6rHO2J9fJ3Kqtvwe3M3ZLwhUYEckkeFjrm6p6qKNF22g7+mZV1dkEK3cqcFFEQ3LqOUWkADhDVV+MdLajG47M7/JmVZ0InOf83BrZlBHJmQTkAe+r6hTgA+CXcZgz1FxgUaSyHbfxU//bTCZYYCYDw4B1BPci4iqnqh50ci4B3gV2AM0u5mx3FW20RfSLRQQydktCFRjnD/t54GlVfcFp3iMiQ53nhxLcK4Hgt4DQbyt5BL9tHaWq9cAy4Jo4zHkOMFVEdhA8THamiLwdZxlR1XLntoZgX9GMSGWMYM79BPcCW4v1c8CUOMzZuq58IElVP4pkxgjmLABQ1a3OHtazwMw4zImq/kVVP6Oq5xA8LLTZxZztKXMyn5Q/jjJ2S8IUGOfQy2PAJlX9VchTy4D5zv35BI8xtrbPFZFUERlJ8Jj2ShHJDPnFJgFXAJ/EW05V/a2qDlPVEQQ73T5V1QviKaOIJIlz9pDzR3wVEMmz3SL1u1SCfW0XOMtdDGyMt5whr7uJKOy9RDBnOTBORAY6y10CbIrDnIhzVqNzltS/AI+6mLNNziGqGhE521nnvM5eE+uM3RaNjqRo/BD8kFWCu+NFzs8VBM8UWU7wm8lyoF/Ia/6D4BklxThnZQCDgVXOejYA/4/gt8W4ynnCOkcQ2bPIIvW7zCB4plPr7/I3tHE2lNs5nfbTgHecdS0HhsdjTue5bcBZ8foectr/mWBRWUewePeP05yLCH6Z2EgEzxY9hZw7gAPAYYJ7LuOc9mkEv5xtBR4iQid3RDjjL5zHAef2R51t367kN8YYExUJc4jMGGNMYrECY4wxJiqswBhjjIkKKzDGGGOiwgqMMcaYqLACY4wxJiqswBhjjIkKKzDGGGOi4v8DfVnExAGooV4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sdatwcs = WebCoverageService('https://webmap.ornl.gov/ogcbroker/wcs')\n",
    "mort = [] # array to store mortality values (Mg C km-2)\n",
    "\n",
    "for i,j in enumerate(odd_i):\n",
    "    layername = '1512_' + str(j)\n",
    "    mortality = sdatwcs.getCoverage(identifier=layername,\n",
    "                           bbox=(longitude-resx,latitude-resy,longitude+resx,latitude+resy),\n",
    "                           crs='EPSG:4326',\n",
    "                           format='XYZ_FLOAT32',\n",
    "                           interpolation='NEAREST',\n",
    "                           resx=resx,\n",
    "                           resy=resy)\n",
    "    \n",
    "    # output from above is in bytes; lets convert to string\n",
    "    mortOut = mortality.read().decode(\"utf-8\")\n",
    "    mortValue = []\n",
    "    \n",
    "    # loop through the output and store mortality values\n",
    "    for k, lines in enumerate(mortOut.splitlines()):\n",
    "        if k > 0: # skipping the header line\n",
    "            mortValue.append(int(lines.split(',')[2]))\n",
    "    \n",
    "    # mean mortality (Mg C km-2)\n",
    "    mort.append(np.mean(mortValue))\n",
    "\n",
    "print('Done!')\n",
    "\n",
    "mortdf = pd.DataFrame(data=mort, index=mortDates, columns=['mortality'])\n",
    "\n",
    "# print the tree mortality due to bark beetle\n",
    "print(mortdf)\n",
    "\n",
    "# plot it\n",
    "mortdf.plot(marker='o')\n",
    "plt.ylabel('Tree Mortality (Mg C km-2)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "q7gB7ooRaitk"
   },
   "source": [
    "### 4. Summary\n",
    "The bark beetle infestation is often triggered by drought, with beetle activity becoming detectable within a year or two. As soon as the tree recovers from drought, the beetle [populations collapse and tree mortality diminishes](https://www.fs.usda.gov/detail/boise/news-events/?cid=STELPRD3841444). In this tutorial, we demonstrated basic methods to access data from three different Web services hosted at the ORNL DAAC. Now, can you put together the three datasets we just downloaded (vegetation index from MODIS, precipitation from Daymet, and tree mortality from an archived dataset) and draw conclusions from them.\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "ornl-daac-webservices.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}