Calc_Sandwich_slopes_and_roughness.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.ndimage\n",
    "import time\n",
    "%matplotlib inline\n",
    "\n",
    "def slope_filter(z,dx,dy,idegrees=False):\n",
    "    \"\"\"\n",
    "    Filter kernal to calculate ArcGIS slope\n",
    "    Call from scipy.ndimage.filters.generic_filter\n",
    "        \n",
    "    Iput:\n",
    "        z - Raster of elevation values\n",
    "        dx, dy - Spacing in x and y directions\n",
    "        idegrees [optional] - when True, converts returned slope to degrees\n",
    "        \n",
    "    Returns:\n",
    "        s - Raster of slopes in fraction (rise/run)\n",
    "        \n",
    "    Use the mode option in generic_filter to treat edges\n",
    "    Suggest mode='nearest'\n",
    "        \n",
    "    See Garrard \"Geoprocessing with Python\", p. 256.\n",
    "    Horn (1981), cited in Burrough, Macdonald, and Lloyd (2015), equations 10.9 and 10.10\n",
    "    \"\"\"\n",
    "    dzdx=((z[2]+2.*z[5]+z[8])-(z[0]+2.*z[3]+z[6]))/(8.*dx)\n",
    "    dzdy=((z[6]+2.*z[7]+z[8])-(z[0]+2.*z[1]+z[2]))/(8.*dy)\n",
    "    # slope in fraction (rise over run)\n",
    "    s = (dzdx**2.+dzdy**2.)**(0.5)\n",
    "    if(idegrees):\n",
    "        s = np.degrees(np.arctan(s))\n",
    "    return s\n",
    "\n",
    "def slope_dir_filter(z,dx,dy):\n",
    "    \"\"\"\n",
    "    Filter kernal to calculate ArcGIS slope direction\n",
    "    Call from scipy.ndimage.filters.generic_filter with size=3\n",
    "    \n",
    "    Iput:\n",
    "        z - Raster of elevation values\n",
    "        dx, dy - Spacing in x and y directions\n",
    "        \n",
    "    Returns:\n",
    "        az - Raster of slope directions in 0-360\n",
    "    \n",
    "    Horn (1981), cited in Burrough, Macdonald, and Lloyd (2015), equations 10.9 and 10.10\n",
    "    \"\"\"\n",
    "    eps = 1e-6\n",
    "    dzdx=((z[2]+2.*z[5]+z[8])-(z[0]+2.*z[3]+z[6]))/(8.*dx)\n",
    "    dzdy=((z[6]+2.*z[7]+z[8])-(z[0]+2.*z[1]+z[2]))/(8.*dy)\n",
    "    # azimuth in geographic notation\n",
    "    # this will return a number even with zero slopes\n",
    "    az=np.degrees( np.arctan2(dzdx,dzdy))\n",
    "    az = (az+360.)%360.\n",
    "    return az\n",
    "\n",
    "def tri_filter(z):\n",
    "    \"\"\"\n",
    "    Filter kernal to calculate topographic roughness index\n",
    "    Call from scipy.ndimage.filters.generic_filter with size=odd_number (typically 3)\n",
    "    https://gis.stackexchange.com/questions/6056/calculating-topographic-ruggedness-index-in-arcgis-desktop\n",
    "    Although this cites Rily et al. (1999), it produces numbers that match the post by whuber and does not\n",
    "    agree with the numbers presented by Riley et al.\n",
    "    \n",
    "    Iput:\n",
    "        z - Raster of elevation values\n",
    "        \n",
    "    Returns:\n",
    "        r - Raster of roughness values\n",
    "    \n",
    "    \"\"\"\n",
    "    i=len(z)//2\n",
    "    r = np.sqrt(np.sum((z-z[i])**2))\n",
    "    return r\n",
    "\n",
    "def mad_filter(z):\n",
    "    \"\"\"\n",
    "    Filter kernal to mean absolute deviation topographic roughness index\n",
    "    Call from scipy.ndimage.filters.generic_filter\n",
    "    \n",
    "    Iput:\n",
    "        z - Raster of elevation values\n",
    "        \n",
    "    Returns:\n",
    "        r - Raster of roughness values\n",
    "    \n",
    "    \"\"\"\n",
    "    i=len(z)//2\n",
    "    r = np.mean(np.abs(z-z[i]))\n",
    "    return r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 13 files.\n"
     ]
    }
   ],
   "source": [
    "# list of files to read\n",
    "fdir_in = \"F:/2017_Karen_Sandwich_maps/\"\n",
    "fdir_out = \"F:/2017_Karen_Sandwich_maps/slopes_and_roughness/\"\n",
    "bnames = (\\\n",
    "          \"2016-01-22_SandwichTNB_PT\",\\\n",
    "          \"2016-01-25_SandwichTNB\",\\\n",
    "          \"2016-02-11_SandwichTNB\",\\\n",
    "          \"2016-03-30_SandwichTNB_AS\",\\\n",
    "          \"2016-03-30_SandwichTNB_PT\",\\\n",
    "          \"2016-09-21_SandwichTNB\",\\\n",
    "          \"2017-01-09_SandwichTNB\",\\\n",
    "          \"2017-01-25_SandwichTNB\",\\\n",
    "          \"2017-02-14_SandwichTNB\",\\\n",
    "          \"2017-03-16_SandwichTNB\",\\\n",
    "          \"2017-04-28_SandwichTNB\",\\\n",
    "          \"2017-05-04_SandwichTNB\",\\\n",
    "          \"2017-09-18_SandwichTNB\",\\\n",
    "          )\n",
    "nf=len(bnames)\n",
    "print(\"Processing {} files.\".format(nf))\n",
    "fnames=[]\n",
    "for i in range(nf):\n",
    "    fnames.append(bnames[i]+\"_DEM_10cm_trimmed.tif\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10226, 'width': 11987, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.09999999999999806, 0.0, 376286.64299908397,\n",
      "       0.0, -0.1, 4625440.8344103005)}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\csherwood\\AppData\\Local\\Continuum\\miniconda3\\envs\\IOOS3\\lib\\site-packages\\ipykernel_launcher.py:55: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time: 2293.811865091324\n",
      "Elapsed time: 2294.911279439926\n",
      "2016-01-25_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10156, 'width': 11420, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.10000000000000205, 0.0, 376287.14850415837,\n",
      "       0.0, -0.1, 4625434.75106331)}\n",
      "Elapsed time: 2205.3239262104034\n",
      "Elapsed time: 2206.4275889396667\n",
      "2016-02-11_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10204, 'width': 11673, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.10000000000000099, 0.0, 376296.13149943855,\n",
      "       0.0, -0.09999999999998174, 4625448.832959766)}\n",
      "Elapsed time: 2242.999901533127\n",
      "Elapsed time: 2244.2254796028137\n",
      "2016-03-30_SandwichTNB_AS_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 9858, 'width': 11504, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.09999999999999899, 0.0, 376292.3872081661,\n",
      "       0.0, -0.1000000000000189, 4625440.284709241)}\n",
      "Elapsed time: 2104.4556589126587\n",
      "Elapsed time: 2105.5313110351562\n",
      "2016-03-30_SandwichTNB_PT_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 9858, 'width': 11504, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.09999999999999899, 0.0, 376292.3872081661,\n",
      "       0.0, -0.1000000000000189, 4625440.284709241)}\n",
      "Elapsed time: 2174.043838739395\n",
      "Elapsed time: 2175.1224870681763\n",
      "2016-09-21_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 12379, 'width': 12221, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.1, 0.0, 376232.5822683347,\n",
      "       0.0, -0.09999999999998495, 4625453.651816779)}\n",
      "Elapsed time: 2831.643162727356\n",
      "Elapsed time: 2833.0549829006195\n",
      "2017-01-09_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10491, 'width': 11301, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.1, 0.0, 376323.5035454615,\n",
      "       0.0, -0.1, 4625436.139639027)}\n",
      "Elapsed time: 2183.004531145096\n",
      "Elapsed time: 2184.284298181534\n",
      "2017-01-25_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10175, 'width': 11191, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.1, 0.0, 376317.60028574534,\n",
      "       0.0, -0.10000000000003662, 4625436.443681666)}\n",
      "Elapsed time: 2154.5384035110474\n",
      "Elapsed time: 2155.6910960674286\n",
      "2017-02-14_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10544, 'width': 10865, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.10000000000000214, 0.0, 376300.70337594853,\n",
      "       0.0, -0.09999999999998234, 4625434.457048389)}\n",
      "Elapsed time: 2306.416239976883\n",
      "Elapsed time: 2307.538923025131\n",
      "2017-03-16_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10578, 'width': 11474, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.09999999999999898, 0.0, 376307.0163862964,\n",
      "       0.0, -0.10000000000001762, 4625455.607236579)}\n",
      "Elapsed time: 2233.764681816101\n",
      "Elapsed time: 2235.049449443817\n",
      "2017-04-28_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10975, 'width': 11993, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.10000000000000098, 0.0, 376255.2562659597,\n",
      "       0.0, -0.10000000000003395, 4625499.890746)}\n",
      "Elapsed time: 2480.957988023758\n",
      "Elapsed time: 2482.245751619339\n",
      "2017-05-04_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 23630, 'width': 24283, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.04999999999999904, 0.0, 376239.69198716746,\n",
      "       0.0, -0.050000000000007885, 4625495.902980097)}\n",
      "Elapsed time: 11636.93161725998\n",
      "Elapsed time: 11643.742728233337\n",
      "2017-09-18_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "{'driver': 'GTiff', 'height': 10236, 'width': 11762, 'count': 1, 'dtype': 'float32', 'crs': CRS({'init': 'epsg:6348'}), 'transform': Affine(0.09999999999999802, 0.0, 376276.8919304822,\n",
      "       0.0, -0.1, 4625437.520928205)}\n",
      "Elapsed time: 2402.8985652923584\n",
      "Elapsed time: 2403.993199110031\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAFpCAYAAACGbcLwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XuYzfX+9/Hne8wwGGMcZ9Q4TEYh\nZ4pq2kapKKEotNmpKO1UKuWYNVOKZJeiEEWZQlRC0aYa0WZyGoeQkRksmcnMYM7m4H3/sdi/fvt3\n/+79u+8b67vWej+uq6uNuXbf1nW9Pr3W5/QVVcUYY4wxxhcFefsBjDHGGGP+X1mRMcYYY4zPsiJj\njDHGGJ9lRcYYY4wxPsuKjDHGGGN8lhUZY4wxxvgsKzLGGGOM8VlWZIwxxhjjs6zIGGOMMcZnWZEx\nxhhjjM/y6SLT+9ALWnPiUW0X95G9Z8EYh5ObUbkBlYZYXo1xOHkRlXGev7z9LP9OsLcf4P/VNx8/\np89VC2ZVizKOF17n7ccxxvwfyGcoqUABuBBw/thoTMC6ml3K/LYQCq4MgSnOzqvPzsicq1DqNKhO\n83ZvUzUnx9uPY4z5P9B+CAWfwTBIjHb2oGhMoDtIWyFzEyyHxMHOz6vPFpmef3lDTmeX8MLnRRx7\ncjH1n9qrYZN/14YPbXb+p25MQOrHt/NuhC5Qn6OaIKJ9WWx5NcaJYuPQ1wR6AStFE0Q0kwhH5tVn\nl5YA4nq9KloxTvW6KbzQYBJr8/LYs/NJbz+WMeZ/QxUR6a86RxDZReIWheeB1739ZMaYf6VpiMiX\n6lomSA+lnh4lcWkjdIC3n+y/ElVHFqz/Kxu/HKc/rTtKdLP2PNMtg5Pfv07Z01XF289ljPmvRFT3\nawwtJB34CcI6o/lYXo1xIGmP6lBBJipEgMstJKg6Kq8+u7T0RzVvi2T4W3HcO6I6r66qQ71uzxM+\nye37Dc0Yv7SBJXKEE1qLOM3DVSCI7LO8GuNEqdtJGAU6UDiTURkAWeqs3fp+MSMz6bcEPV1ezsj6\n9Tm74xRfzjvAyhei2Jr6KDqohaOaozGBThqgZAKs5DndzSGaUkEwqzfch3a1mRljnERuRjkAZLv5\nSXtSQSWKqMZIZrGPDo7Iq18UmQuePfaijo6MZOO7+0FiKPxzKcOOHOFq110cWN3bER+4McZDXkGZ\nuBQmD4CJbmgejU4Q3EMg2mFT18YEOtmAEp8CGR1gRAhEgWuhJ6beXmryi6WlC95o+LJMzcyk/Yhr\naHNjMcdHp5LTsAXhr39P0Lvl2qnTe/7T2ozxcToBodMA4iasg73R0ATiB6+h4QwlkTGaIGJ5NcYh\ntCvCiM781LgDrdZshR7QTa/nE01lk4hX8+pXRQbg7UaTZWpmJulXh/Lge3/iy3mnmHm0OnfcNo7t\nfztF1ddOafNeK22ANMYBdCuy52xrCCuDbIiXnjAHEqQ7iRQioirdnLUeb0yg0tnI19xJFpGwAr6X\nn0h7vy038zOJUYoMQeXXy7/M41dLS3/01NGJWoYypHYdGmdVkFL7HKnFxbx8pBgOjaTpl/s5tPx2\nm742xgFkOIobWJsEswbDyJUwsjf6jSBphyE6Bj1m+2eMcYLObNCf/tYVRu8jrKAeBWEFsC0GVych\nsZfCadCNly+vfltkLuiwb5QOqF2bP+08S5sbGzCrOJdv8vL4Pu0GQn65g/o/78P9/vU2QBrjAIdp\noG0Ld1PQtx4J68aQ0PA1cBdBXDXYBkSAnrBCY4wj3CX64lfjmPzhqxQMrETYwArYBAwDlgNNQNdd\n+rz6fZEBqPxGgVbJz6fgxtX81uwMRfllLK1byraiIr7YewtVf+7IFTt22AyNMQ4gb6FXPf0zAJFk\nAbBZgoDGwI/AYOgFusoKjTHeJkfK9KfGHdhHS2qQTwSnuVU+BqKBHyEqHvp6lqUu2TMEQpEBeME9\nSbPLy/lk2Vh+vXcuodWCqR1ZjQW5OYx2uzm18xWC80JpkpxM2hc9bYA0xoukJ0owsDoPlodD/30w\nsSVM3gSpcTAQoven4ZYUVAdbXo3xInkcpQkwH4I2FXIuSiC5GsTn0UF3s+O1OB4f8waz5UFU61z0\nvPrdZt//zrTol6RucDBVzpwmetvtzJd8ilFuSK9gdqNGLOkzl/J2b5PTrBnhk9x61f3rA6PhGeNA\nugahP7AkHMYCREIc0CMO2q2DA3m4pRldtS5J9Fe7UM8Y79HZCO2A1XBucHWgGm27boEZ4exoHwcz\nYHatZxmp75NJhIpc3A38ATMjc0Hd0Qe17oEDZHSNp3L3MeTvHY/7xiUsDSsh4bffyN83mrrbSjn5\npl2kZ4y3dWaD9pR4EsnkUx3G/T+vglZAKFCSBsQCn6I6wPJqjJdN40ktklkk1lX0QSHilROcaRcF\no4FhChECp9NQbXZR8xpwRQZAPsxSIlKhqBENvi+kclERtcd/zs6f/0yVjGsoeT7cBkVjHKJSVqGe\nS6jOotn9GSLjYGxHiAIOAEtAT9leGWOc4i0e01Hvz0XvFmSCEjI1D4CyzHAYcWlOMwVkkQF49USi\nNq1ShZTCQvYWF5M+5lbO1qjB0Y/ibFA0xmFE1igjekJzqPf0UU4OacSqRbfSi28tr8Y4jNyMMgro\nUkKLK/ex/8cO6BUCMZfmBuCALTIAZ4vf1aAgIaiSsKGogFvCR9ugaIxT7RWVuxUygB7n99EYYxyp\nLEK08hCFLXhOGbrs1NIl9UP+3/RPNZ6zQdEYHyCSphd7jd0Yc2lICKpll/ZLhxUZY4wxxvisgDl+\nbYwxxhj/Y0XGGGOMMT7LiowxxhhjfJZPF5mnjk60DT7G+IgGHLa8GuND4lnrE5n15SKj/baXefsZ\njDH/AwkiminfePsxjDH/QwkiukGyvf0Y/yM+fWrplx2T9fnwPFbGTrOjmMY4nNyA0gl0pt3/Yowv\nkOdRBoO2dXZmfXlGBoCbpp+k/q4nfLeNGRMgkjd3xjVLkNsu7gvjjDGXhkYJrnaC/OLszPr0jAzA\nptXjNa7Xq+wumkHfX3/lcOuZjm6OxgQykRx1UReAxFBFi539Tc+YQCdVUVeJJ6aJcXpJ3pX0/8vn\niwx4lpjSdudQfrYpz974i5UZYxxMBCUacL8CvIBqiOXVGAeTZijxwPyV0KU3utlZZcbnl5YAKlrU\n4IfuNbjhYWHahqu4YvdI329nxvirvkAcwIOAmxqFJy2vxjjZZKAXQE84AImMcVRm/aLIJOXmkn+u\nggpVOt9exts/NOX6/c846oM2xpy3FhgFEALEUDCnHgkilldjnGoO9OqzDAiGvpBwZLKjMusXRWZb\nYSFVJYiXM09Q3CiUeldmMjY5hhfck1S++dAxH7YxBs+3u9MA+2EY0AkSmYzIbuUl5wyOxhiP6O/T\nyKEOsB/6wu2NvyaRqUg91AmFxi+KzN+vni7Feo4QhOlZWTTv1Yhb71PqT/mN31vlcNWeJ73+QRtj\nPPQ5hAJgaDzMB5YDMybgoi30BrnO2SckjAk0x2gmh4iFOS1hFPz9195E6X24soVYQJ7xbmb9YrPv\nHz11dKIGixAXFkaX6tX54LEfaRsXybt/CmHtoQ5oz0GO2qRkTCCTqihdgOQcrtJMDt93LbiBLQpR\ngp5w1qZCYwKdtEfpD0yFd/OH8te5CyEVmANEgx67/Jn1uyID8OyxFxUgs6yM6dHR1Cg8R371IKZm\nZrLk1ClOHngUHdDWBkhjHEDuQVkLTAS2AAeATsCSV4DWQDGqAyyvxjiEPI+yF2qvPk5uwpWQgWcz\n8MA1QGPgd1TjL1tm/bLI/FHPtOe1RqVKvJIXTvaJQvZfV52k3Fy+P1kV9k1Ch9vRT2OcQh5HmZME\n9GSvXk+rWr/C6ZUQ3BvKQdVmaIxxEvkWpXs6RMSgvQQ5oLBtO7TrCCWg+y99Zv2+yABUfe2UApRE\nlVI68HPWJv1CWq86rD1zhnVZNQnZ+gylz4bZAGmMAyzjbv1ZVpPIYegeA9lA6j6gGYSFQMEaVHta\nXo1xiodEExZCYrDCQiAZmJ8HweGeqxb2gp68dIUmIIoMwF/SxynAsk9foPCFT9lbXEzIz/l831iY\nn51NavJ0KpWVUfZ0VRsgjfEykRSFzlAXz5LTKDef6aP0k/cgNBpKkiE2Hk2zGRpjnEAkS2keCdHA\nCKA/lJ4WKkcoRACn06BXM3TVxc+sX5xa+p/4KGaKRIWEUP3kSSp9MJiFOTk071ifFj/kMatRI0b2\nSqAibhQAsmxrYLQ7YxxKtbNcpT+jPwuMWgiE0E9GQP9oiAdoDIc8P7uDlpZXY7xMNVJ67v8c13qB\n/skAVI7IhBnAVIBmsPf8D4+5uEe2A6bIAJScO0dZaCih2aW8u+R5gr4N4hbpQ2T6Wa6sXJn4GjWQ\nxftV77vOvuUZ42WH5VoSIoHBQ2FbJIzqTc2kTFjrhi4xQBFSFe3APsurMQ6wpv29nv8xP55WuhVW\nR/Lc05NhNDAScIO0QHlNL2pmA6rIvN1oshTVq0d5aChno4/D6XbUSonkrnHX8nNxMZsKCtBBLWxQ\nNMYBVJGbAMqhbcctJLw5hjObonDREAqAJtXspZPGOIjuPJ/HaGjNHk7cFcHcwsdwFYjnWoXul2bz\nb8DskflXjx6ZoPE1arCxIJ/CinN8/1Jvji24wQZFYxxIZKXSo7fn9QYshOShaFcrMcY4lchuZWIb\nz7JS+Rr0lTth/MWdifnnPytQiwzAPwre0FZVq1JdhUrBI2xQNMbB5BmU03D7gi/5hj6WV2McTuai\nFMCK5+6gD9/YqaVL5e3fJ+tT9SfaoGiMD5B66KU8xmmMubjkPlSXXdrMBnyRMcYYY4zvCqjNvsYY\nY4zxL1ZkjDHGGOOzrMgYY4wxxmf5dJEJnlWqsmqVNnxos230McbhOrNBE0Q0QS7urZ7GmEtjKLN9\nIrM+vdn33l/H6JrlYykeU8tOMRjjcCIobEI1zvJqjA+Qt1BG7UO1paMz67MzMl/MeUqzy8spqQ/1\nntnvu23MmAAgN6N0gpolsZ67JYwxjibDUbbA9XqSeNY6OrM+PSOT8s1E7VLaGX7rhT4mjm6MxgQ6\nqYqSCjS3WRljfIHEoIPSP2DxkIfRRc69v8lnZ2QAMvafoqLXb9B8GrJ4v9Yal+67rcwYP+cqEUgC\nKEbErSKzLa/GOJgrQ1i8/WFIykKqoiJZjsysT8/IAKxe+Izu/Ucu746swbH8WoSkjKH02TDHNkdj\nAplIkTKnGiwH1icBzVDtbHk1xqFE0Cg9TOa4q2DqdiAa1UhHZdanZ2QAbrqzMTfddQUvry7llauC\nKWs/D5uZMcaZBuvHMBWCkgph4WAgA5HtlldjHOpVfYbMQVfRYcomWqgAaThtZsbnZ2QAis7N0SW5\nueg7h6n1VDP6pVam7soenHyzhaNaozEGRNYptOF7oug2TGE+MBR0gXPX4I0JZBKCUg4uhMRtCiOA\niaB9nJFZn5+RARhx5AgrTp+m1ehrKZx7mG87VpB9/0LCJv/u+y3NGH/T/DYe1M/pxvdcP28D9Icl\nCwWRFMurMU40HTZqRxLZxaiOU2AUuPoKIjmOyKxfFJmPYqYIwLd5efx5TAc2PpFCbtMGXNF3MrJ0\nl9Z/aq8jPmxjDOh+5ENpwi1aQMqP8SxcNoD3dBXs7YwIKvfY8WxjnESfRm5uvx23dqOmjEcbCD20\nLdGai7RGZbt3M+sXS0sXPH50gtYLDqF7jRqkz/yFBk3iSO76K1MyzhG2/q/kuxo4YhrMGAOTGa1v\nVjxDbvBnsOUp6LIbpraBWYDbjWq05dUYJ9kr+lirGbwn3eiph1gj3alZUsSZ+ChIBS32zlKTXxUZ\ngL+kj9PQoCDeiYqmpKiM9xP385fprWmyZw/5uxPQB52129qYQCbtUboDBef/SgWmAz0UyEA1xvJq\njIPIOJTmQF0gGMiEGx78js0NbwG3onr573TzuyID8OiRCQrQqmpVhlerxe4fT1Cra31GHj3KupxK\nkDEUHWQbgY1xAolBiQJKgNQ0oCqQC13awJZ1uLidBFXLqzEOIa94Nv8SBcwABgIJRTCqGszYjYu2\nlzWzfllkAF5wT1J3aSnbioqY37gxeZ+5+dOQqxl59Chbi4o4cLQLNTbfTN5L0VJz4lE9M7mRDZTG\neIlImTI0BCKAGUnAYGgFZAAF+4BwoDaq1SRBRK3YGONdUgtlOYS0y6Osbg60ivmPzG4pA0IgDDSf\nS55Zvy0yF9QdfVBDCgupCAmh1ciFTEkP5YqYcF6ucobVZ85w4lg3aqR0Ju+l/1iPl2/fUb31CRso\njbnMJjNaX6w1HSYCscAKCJt1koKwX4D9QE8Ijf7Pa/FpojSzYmOMV/wuKgMUFpbRoXEKO56PI+H1\nMSTIn4EsIB46haBb/yOzF7vY+H2RAc/sTEZpKcuORFLj+/tIuXsBxdfWILZKFUYdO8bavDxOpLxJ\n9axTfPTEAu6tNcYGRWO8RCRNmdMMXSzIhp95QZcxTa6HET0hG1ieDsQA6ei8q2CYlRhjvEkEZQks\nGSgMpJD39DEelbmwqRosAWYptBNI9dxFc7FnZwKiyAC8kfWyHigpYf3z3UgfVIl3bvqRW91QO7Ia\n66qWsvzUKbYUFpK3chIFE+vbwGiMF4mowqdAd+hbB4bBe3cN4dGPF8HgMiAXgiOhHFSdcSmXMYFM\nnvQsNZEJzIK4J9ax8aXbaTjpIO6qzTx74PoCm0BPXtzMBkyRuSBq5G4NLi7h92tb0uUuF5llZSw6\nHs7CFsovJWf5fvPzVPm9KucqBdk7m4zxMhmCMhFoB5S4ISLas7lwIuBOAm46/5N1UA23vBrjZfIl\nWhovVO6unpOIwcBkYCGwYjsEd4RyoO7FKzR+cSHe/43MWW0kv0EUddLS+G3CHVQb0YkHZ15PSkEh\nBRUVVMqvBWAlxhgH0EXIyGumQSdYqM/BYCBOwZ0CxAMbgBgrMcY4hPZBDta8CoaBviXoSOH+Ph/C\nijygI5QnQ/TFnZUJuBmZP3r1RKKmFhfzxSoXbRct4kzDhhRERXFidnsbFI1xGJFkhaqQ3BndK0iy\nErdsHRu5zfJqjAOJuBXC6aNf8Zg8wJ3Tldeee5IXmGlLSxdT/qm3dV/lCn4sKGBhTg67W86wQdEY\nh5IYlAJgFDDVc7TT289kjPnvyT0oJXguupwI+sXFz2zAF5kLzha/q1Wq/tUGRWN8gCSi6rISY4yv\nkG9QvePSZNaKjDHGGGN8VsBt9jXGGGOM/7AiY4wxxhifZUXGGGOMMT7Lp4vMS78laPCsUtvkY4wP\nkF2otMfyaowPqJRVqLLUN/Lqs0Um5e9rNFiEishdRI3c7RMftjGBKkFEB7X9AMpB3vKNwdGYQJUg\noo9FzoUkkM+cn1efPrW0a2OCtqMSHH0A/XNTO4ppjINJCEq5G/gQ1QmWV2McTLqhJOcBaah2dHRe\nfXZG5oKPfnRTo+V0ZNUqrTUu3XdbmTH+bhbADlg4gSeZpiLbLa/GONXCMiCDKK3FMu5WkX2OzatP\nz8icpzNGPc7uUbVYkF6PsB8Gku9q4Oj2aEygShDRRBRC8bwNl32otrS8GuNACSKa2EM9b7QOBvaC\nFjvvIkqfn5EBpGO3CK6Zkk7OzVUouGU+tV845PPtzBh/lKAqAK4SYaKOh/ktkeucvwZvTCBKUBUO\ngCtV0BlCn+LFnjfSO4w/FBm69GhE62kd+PrN3fzcKZyCgdOtzBjjUBlan0QWUEmmQDKwbR3yivMG\nR2MMaHchkeUkxMGXGwZB0m5kg7Py6hdFZsLJEyzJzaXNE8359ZN0jrVuTeSDM63MGONATRr+DmFD\nSeTvsBdw3wzR0I4tlldjHEYOKbTqRyLHIAkS9GPqdT3KeCY5Jq9+UWSmRb8kJ8vLSS4ooOvQa5jz\ncDIrztajx4j5yIdZjvmwjTGgxxAK0oE9kAocCIW6sOueLshXzvqmZ0yg0+8R9ipwJayF+QyjMylM\n+fglHuADR+TVHzb7/tMDh8dqn4gI2lerxvZ393PyeCUYE8Vot5uyE7ehd/Vz3CYlYwKViFshEkiB\n9XHQPf38nzRhoQ7kQZZaXo1xCBFVogXcWdyuW/i79AGKoEk1dJTA0+q1vPpVkQFPmQmrVImXr7iC\n4OxSvph7kkcmNabx3j0cS3kTHR5ig6MxDiGyVKE7kAXzW+IaJqzSjeyQSCAY1RjLqzEOIbJJIQ6a\nAJtKcEVXJVif48Vm06Eu6GbvnGjyuyID8Jf0cRoaFER8jRrcU7kGqxfsp+3wqxntdrP6zBmCN/yN\n0mfDbIA0xgFEkpRWg2FvGfQKgQLwXMSVAhRDcG+0zHlHPo0JRCJuZWA0LAGS8CwPTwfY7fmBTm3Q\nrZc3r35ZZMAzM1MzuBK/lpzl5SuvpHT979S8LZKk3FzW5+Wxc89D6IC2Njga4wAinylx/WBTGfQN\ngQw8AyTJQAtoEommW5kxxglEUpT+nWEFMPb8byYBGekQFgO9QBdfvrz6xWbf/51PrpoqZ8orKAe6\nfDqUJ5oVU3DuHDd9fYYVsbE81HkpjXaPVFm6y148aYyXqfYT+gIDQ6AXMB9u0O+gSTyMioRgkBdR\nkX2aIGJ5NcaLVDsLk4GhwGDo8PImJqaPhxExMAdwg3yJipRdlrz67YzMBb0PvaDBIgSLsHrJaEoi\nIvih0SSuvb4+a4JLSC0q4vVfIqmyYwglz4fbNz5jvMjz2oKfYMbjsARu3/wl/yi8kYKIelA+DygG\nboK6HdGTNkNjjDdJDErGOki9DVbAdldLvuZOXrxuOmzbBOQDN0GvcHTVpcur3xeZCy5sAp73wwgI\nLuDLnqvoHh5OvruQb6qXsq2wkJm7OsO5ymifO2yANMaLRHIUsvhExzOOqRyRXXhKTBZEjYHBwF5w\nrZV/3hZsjPEOERTcHNQbGMFsvpMGQFXgR+g+3JPXcnANuzR5DZgiA57ZmbrBwXy0OoGKpp8xtX0a\no2rW5fPZe0gbGEl2eTkzd94MmbdT/eRJCibWtwHSGC8RUWWhsOrBW7l77rcwFcjIg4hwOJ0G7Iap\n/c6v0e9GtY3l1RgvkQYoK0C/FRqOP4h7UDPPhuB4YBNQvgnWxsFoYG8aqs0uWl4Dqshc0KnTe//8\nlz5bowY/33cf8++fxvq8PJacOkXl7163ZSZjHEIeQlk4mxbamf3fdPC8vG4gkJ0OwTFQPg/V4ZZX\nYxxA3vJcarnw6QG8yTO0ZB+Lb3sY1iu0Eti7DtXbLmpeA7LIAMiHWXrd9VPZerAH0SvCidyzh+0T\nmkF2nN01Y4zDyFsoScC27Xgu0SsDQoDjno2HxhjHkG9RlgNzcoA6ni8f5QAXdybmn/+8QC0yFxws\neVtn/f47AMsn9+W39zrZoGiMQ8lrKGO3Q0RHiAbdYxt+jXEq+RElrgi6V4NhoAMuTV4DvsgAHD47\nUyNLhOo1R9qgaIzDSQuUA+l2668xPkCeRJmlqIqdWjLGGGOM+Vd+eyGeMcYYY/yfFRljjDHG+Cwr\nMsYYY4zxWT5dZB44PFZrJJ6wTT7G+AD5BpVuWF6N8QF1Ko6rvOMbefXZIvP1h/O1bnAwBVeVE/3I\nTz7xYRsTqBJElAggFWSQbwyOxgSqBBGNrXQIZoA87vy8+vSppf1bX9aWebXB3R99MNKOYhrjYCK7\nlbA2kATax+5/McbJpAZKFESnpXGMi3+J3cXkszMyAKUlFSzasYeQ1i5k6S4Nn+T23VZmjN+LhII1\nMPD8DM2rYnk1xqkKgEPbcT/ejAQRTRDn5tWnZ2TO0zeeeoy9z9ZhQXo9qm4eQtH4uo5uj8YEqgQR\nTUShC7AlDYhE1d5rZowTJYhoYpPzeV0CdALd6rzZVH8oMmz8cpzu25pFX1dHorbkUmvVIHKnxTru\nwzbGgEiWuojijL7KWaqQQRO+5l7LqzEOJIK6EBIOQmazmvyDG7mXrx2VV78oMnpurq7KO0PevMO0\neaI5nQ4coP47D+N+/3pHfdjGGMgkQhvUOI2rQEicpTASXAgJqpZXY5xmhKhsUFwHhBklJzgzNArX\nEmfl1S+KzAvuSXq6ooIR9erx26dH6PBAU25NS+P40ufJeynaMR+2MQZE3EqXaNgLRAGDgekQkpFH\naR1bZjLGSaQ1ymRgMlAXQpfkUjK0Nu99MYThLHJEXn16s+8F06JfktMVFWwrLOTGwc14/9EfWFkW\nSY+hs5BFx3y/qRnjR1SjhS3boaAMDuXxvusBFuX3pywqnDf4q+XVGAfRPQgJ53+xCX6sGYc2FR5t\nsAiOOWMDsF/MyFzw6JEJ2qtmTWKrVGHX7ANkHYWg8Vcy2u2m7FRr9NYnHNEejTHQgMOaKT8AsTA/\nDoatBHJppdey48z1hNR0ztS1MYGuI5t0h8QB6USpkim/AfmM0828+tDLsMB7efWLGZkL3mv8iqw+\nc4b52dncNPJahoy5msozTlPSvj0x9Q8g88r8p7UZ4+NOcJWQNBRoQ4tHdkBCb1w8xC8511B5oHqW\noIwxjrCdOCGjDHrFcA8rYFscLu6kPTuR/oqI92ZF/GpG5oJHj0zQSgI3h9WgX9VwVs7fR9vhV/P0\nsWN8k5cHKe9w7q/B9m3PGAcQSVNaNYO9aRDbzPOb0UByOvSPgVTQNOcd+TQmEMkglNPA2izoG+k5\nmn0IWA/MAKaDbry8efXLIgPwcMZ4DQsK4kBJCRMbNCBoQw7ht9QnKTeX9Xl57NzzEDqgrQ2OxjiA\nyHYltCOUpEBYZygogrBqUOAGdgA9UQ2xvBrjADIcJRtYkQPxdSAVaAJkApn7IKoleuLylRm/Wlr6\now+avCol50ta1+UPMzymkNMVFdz8TR4rYmN5qPNSGu0eqbJ0l8rHv/pnmzPGR6h2FCKA4M4wDHAH\nUS//KMRHw4jeEBGCPIlKDRx9w6gxgUDnIcQC3evAWLj/1Ic8vvMNGAhsagnxIF+i0uLy5NVvZ2Qu\n+Ev6OC1XpRz4+pNnKYmIILlJAtdeX581wSWkFhXxelotqmx7hJLn7einMd4k21E6AeyGEW34fXYN\nbuE79sp1539iJUT1hjBbbjLG25pwQI9c1xy24XmH2kFhdeIt3N36WzgNuHdDrzYQDPrFpcur3xeZ\nC/6SPk6rVwpizoYREFTKl3eskCxlAAAgAElEQVR9RvfwcApPFLGqSglbiwqZs/NPcK4yek83GyCN\n8SL5G8omSPhiDAkpr0GXJKADEA4cgi7xkAquEmddzGVMIJLjJcrUUH6a2ZpneYNNAnDb+T9Nh4me\nvW6u1ZcmrwFTZAAeODxWawZXYu7qRLTpx7zSJoPna9dn+Tu72TOgPgUVFczcfity4mbCMjPtMj1j\nvEhiUIIhJ60qdW4o9mwAXp4MdeMhHli+G3q0gQzgwFJUB1hejfESeQglDPRLQR5TaAX0TYP4ZjAU\nGOqGhGhwA/PTUY25aHkNqCJzQadO7ylAeWgo5aGhHOjbl/cGTmdtXh7LcvKp/o2Lgon1bVA0xgHk\nNZRZQDC8n/4Aj6R8Al3S8TSbECAL1UjLqzEOsIy79f6eq6AdFCRU4qYqP7JLWkJYONQFMsou+sb9\ngCwyALJsq97Wdinrfu1I1MrmRKWmkvpCE+TEzXY02xiHEclSOkXCtqUQMQDa4TmezXpUh1tejXEQ\neQglCpiaDN3jPa8iWQtkp6Da2ZaWLraic3N0RlYWJarMe7kPJ2a3t0HRGIf6K2/obGkN7W7j5Z2j\nmch0y6sxDvUBD+gjMhNm1EHvEGh+afazBXyRATinczXnt0LqXfmsDYrGOJxIkhI/GP3eTi0Z43Qi\nKKNBX7dTS8YYY4wx/4XfXohnjDHGGP9nRcYYY4wxPsuKjDHGGGN8lhUZY4wxxvgsvyky9Xc9ofWe\n2W87l43xAfIkKvKZ5dUYHyBLUZFkx+bVX4qMNqlcmew2Ud5+DmPMv1F4tpKyHljYz9uPYoz5d1qI\nsglqljT39pP8t/zm+HVpSbFWSXkbjj6ADmlo90sY42C7uFrb3XcQSkBX2X0wxjjaNlE5ohAFepPz\n8uovMzJUDq0qb36RwaCb3kG++F4b/WWTfzQ0Y/xQWw6Ka7nAahBJUxEsr8Y4VScVV3+BkZ4L7qSZ\ns/LqN0UGIDikMvPqjefLbls5dnNnbz+OMebfcCEQ2gxCvf0kxph/x5UqnvecOYzfLC1d8IJ7kk6+\n4gp2FhfTZfUtNF94mP1r7nHcVJgxBqQF2mH/JnZIPvTtCSt2o9rG8mqMA8l96GfL7qSfvAxJHWFw\nDqp1vJ5Xv5qRGX/cpSXnzpGemk3mxxnk9PuRjFe+odk9a/yrrRnjB+Q+lAMruVtuxsWd0AsIbYOT\nT0cYE6jkRZTlKeyWNbjoROM/H4DmdRBxez2vfjEj897JV/TQ2bNklpVRokowMKdxY3avOcaPq9LZ\nNaYxi398wjYBG+MA9TmqJ6UAmrRkUPoHLJa6sKI3N/T5js3SCeLDoRx0o/M2FRoTaOJZqxvkDhgl\n5LxZlTqyhtDTbfiqZi9uHf4PaA5kg07xXl59fkbm4YzxurO4iAMlJZSoEiqez3LMcTeHb6jB0Flx\n3Dg3k7fu+BD5/k1tMvgH329uxvgoeQg9KY2gS0uIgMUfPgy0hoGwWSpD3XBe+/5J2AbyN1RknuXV\nGC+Rt9ANLXpAf4G60LjwKLSKp6RHbW6VjyECCkZWgrVQKavQa3dD+XyR+aDJqwIQFRLCsg1jOFle\nTjlQWHGOzYUF1AsOpvfLnei8s4Qlbc9Se/znyFefaYueX9gAacxlpgsQSIISKE0XtInQTPNgFMBx\n6AVj7p4JJcDo3dB8OLIBlfucdUrCmECgTyMcWAmh4JooPBdWn9Q9V8McgBiiX08jbE4FAOeiKmBY\nP+pUHFf5+PLm1S+Wlh4/OkFDEArOneOTZWO5v/+U//TnwSJcExpKl82FxN0dw8jfjjHnRBGhyZMI\ny8zk5JstbArbmMtEaqGUACVZ0C4SUpcCTQA39OpH0PxCzq2vDgnAYCAYmJgCfTuzfIXQT9Xyasxl\nIt1QDp3/RS9gTjJwE5ABM5rR5+nFZBDDrkFdIAHqXXOUkxIBCeG4EoSEy5BXn5+RAdh6b2Oyy8sp\nV+X+/lMIDQqieqUgQoOCWLT+JbLLy0nOz6fn/ucJ/qIjLd46w2dXR/Hkn98k+8YC6j2zX5vds0ZD\n3ir2/VZnjMO5TguMBYiEAmDqAIjvDE36QRi8GBXGjD8/BodWQjZkThTgSmgH/VHPPRay2+6eMeYy\ncCUL1x/bAFFAGLAlHhaGQK9mEAHt5AFSU26AJSmE1M3jr9IYCIdhJSR2UqQWKlKk8tWly6tfzMgA\nPHV0omaWlxMMLF43hYfumMDyD/5K/4ffpUqQMGfjwwyP+4ByVVacPs3x2GspLSlne+UyXs/KYu1v\ndYlc1pOsQcvQm162b3zGXEJSDyUW2JIHceHQCUgFksugbwjFSULVMCUkO4+yLeFsv6slHe/ZB4eA\nAtC7BHnnS/TOPvCVzdAYcynJEDynaBaWwegQgkYXcm5Jdc+S8BzQrwTZqgw68QGptGPf6I5IhUIG\nEAauJCGRLbjocklmaPymyAAcLHlbs8vLWZSbw9zViQy9c9J/+ZlyVX568gZ+ua8zr9w+h1ZVq7Kt\nqIglubmk/TKAsL1Xke9qYAOjMZeYSIrSpTNsmcdenUor2QBNoiFjDVD1/E/FAushdigcgkH6AYvl\nLmAHsB94FlU73WTMpSaSroyKgRnJ6GvdkAkKQ4H5aUAkUBWCQ6B8DQzrCVsgZ09V6kgxkA4UQ1hL\nNP/i59UvlpYuuDr0Kbkx7FnJKa9AY5YBUCVIWJAWQ7AIyfn5fDf5HrJatwbgpaRnmNQriqScHNKK\nK6Ag1kqMMZeJamfRzQithrNMDkOPaBgI43QzUBsX3WBiNJAHC4F4uJF/QHAkhPUErrcSY8xlohoj\n+ibC4HgSxuDZwzYU9POrgWBcVKZmQSZQRM95n8MIqL2sxHM8u1MMEHJJSgz42YzMH9376xhdfeYM\n8WFhRAQHExYUxE/D2nH0xhvp/uBMYqtUYVT9+kzPyiIpN5ff275jA6IxXiK1UO0vkAOyU8GNZz3+\nNLAET5FZ+zZQFUYOR2dagTHGWyQG3ZghrAcShyqkel5fkBimsBZYDUxNAdywtx967aXNq98WGfCc\nZloxpT9l1atz5xMLWPvmn8lu3pyBPSbSJyKCjNJSnqlfn8pBI2xQNMbL5BuUCKDLPuAQdOnt2ROT\nvQ9atfSsx68HXWwlxhhvu5PPtTvreU7eBdJhZAxsA7YopApk43lb9iUuMeDnRQY8ZWbOr3UZ1DiL\nrz74K3lNa9Bg86/0mvAFYUFBvNHQNvYa4xSek0g5wHEIbQN7FWJzgfWQcS/aOMTyaoxDiKCEAiXp\nEBHDo6fe4r2Qp6E8hWYaykHaXpa8+n2RueCBw2N1cZYS/V5fsps3B6B4TC0bFI1xIGmBHjsgNOQg\ncByi49FjNhNjjBPJM6hrhpCIAvugS0t08+XLa8AUGYCQt4r1ob6TWTW1Pydmt7dB0RgHky9R+i6F\n7gPQdVZijHGyiLOZeib0OMzviD5yefMaUEUG4Oq9T+nBVm/boGiMD0gQ0ctxM6gx5v+ft/IacEXG\nGGOMMf7Dr+6RMcYYY0xgsSJjjDHGGJ9lRcYYY4wxPsuni8zJ42+oLNuqUSN320YfY5xumqj8jMoz\n9tZqYxzvIdGm/Kwy1/l59ekiU+/KZyEilYKoKG8/ijHm37jzhc8gVCHT209ijPl3PlgwiN6sBB/4\nz6tPFxlA3mqbRWH7jdQal+741mhMIPuae4XY43Dac0zT289jjPnvPcwnMmPIOFjr/Lz6epGh4ee/\n82PXI5zu8Smy4JSjP2xjAt6caJgKicGKfOX8KWtjAlnjRQe4f/aHJK5XnmSaY/PqF/fIfPjqo3rs\n0Bk+H92AtGVjyHc1sAu0jHEoeR/lEGiUICj6tN3aa4xTNeGAHklpjquLkLhF0c7Oy6vPz8gA3P9U\nW+56sDkD3sqi/cDpNPrLJt9vZ8b4KdcwQfsJ0rcUUi+8KNIY40RDpQWHuwjb9VPYAnKP8/LqFzMy\nAI8emaATGzTg02dS2fd8fVa/8QC/v93Kcc3RGAMiyUqreAiFiVvHM1m6oNrb8mqMA4moMlogE35d\n1ICm951AlzlnZsYvZmQASs6dY1thIbcNjObGhTBw9BIqv1HgHy3NGH/TJB66wMit05g87lUgEmno\nvG96xhhgqEBfeH/RAzTdfgKW5yE9nZNXvykyH8VMkeT8fGLb1OXPo6+kzQfnOPDEAiTlGcd82MYY\nD01HyIZZx5+CqftI1qfAvY+GpFlejXEYXYCwFw4RC51Ah9aEtTn0ZbEj8uo3S0t/NP64S0dHRrI8\n8Qhnco5xbOxVzFv8HMVjajlmKswY4yGSrhACC6MhAmgCtEtDtZnl1RiHubCnrXb5cdpV2kldcvhU\n7kE13Gt59csiA54yM7BWLTa/upuW10Vyqmst+h8+TNXPx5PX/hja7yYbJI1xCBGUeOAA0B1OLIqg\nQevTkA2lxULIabW8GuMQLdmh+x/qAG5gLJy5S6j5ikIquJKEBL28efXbIgPw+NEJWlhxjpk1PFcT\nLtUCdhYXMWdXF2ptbEHutFgbHI1xCHkGZcY6GHsbxALDANYAzYEYVJ2zudCYQCcpKF2AGRAyOI+y\n5uGQ7YbQaAC0+PLl1a+LzB8tzZ2qV+0oJDq2JjND8lmfl8fWLZPgXGX0keo2QBrjJHNFpZNCGNB8\nH9ASyEET68Ikm50xxkkSRHSzruAsVdggd0CEQAG4yi/P7IzfbPb9d5IL8um25UlmhuRz3yGY2agR\nU3vOhFYTkbmqtcala82JRwOj1RnjcNJFYS9oZYE5LYE8aFIHWaTIfahIkTrp1IQxgSxSH+QYDRkj\nPWGLwGmFGZDYSZG5qNRC5W+XLq8BMyMDnrtmPp/5IG+P/5yupyoBIA2rMj87G9fO1nC6HdXTa1DU\n80XOdZxr3/qM8SKRMmV6CCQDwcAmPFPX7AG6AoegeRtcBy7/mrwx5j+T1uhVe37m8DvXelaDRwF7\nYaKOZ3K9VyE7D4aF45p/8fMaMDMyAJUEKkJCeOL1AUQfjiC0WjC1SoXY0FCW3JAGYYcovOljKzHG\nOEIGjN4Nq/EsMc2HmiXBMKonUI0Ommclxhin2JvHYbkakjy/7LpnLQt1AJN3vQphMErfuSQlBgJs\nRib6kZ+01uFfKYiKok5aGgBFtWszMukA7rIyphyMgMwe6KAWNjAa42Uia1RvuhP58e+8px9xXDwj\nZCJfA7WBzkAKqp0tr8Z4WUPS9BG5msQuij4rJNzv+f1EtgHREBwJ5VmoRlqRuRiCtj+mg2vX5kBJ\nCXu/mEDLzz+nICqK002aUDkvj6MfxdnAaIxDyJMos5IgdjAcWoOLOwFI5BwgdprJGAeRX1WJ3QHD\nOsL8HFzUBSAxwrN5X49d/LwGZJG5oPTcHN1UUEBGaSnzs7M59tp9HFtwgw2KxjiQPISy8A2gNRBJ\nB81jO/alwxgnkqUoA9OAYqAJi/RhBrP8kuQ1oIsMwMrT0/SusHDycks4WF3pXP0ZGxiNcSiRdMXd\nAKJ/BvJRjbe8GuNQIugN+h2b5Rbg0t3WHfBF5oLdRTO0TbVRNiga4wNGM1mnM9HyaowPWExfHcSK\nS5ZXKzLGGGOM8VkBdfzaGGOMMf7FiowxxhhjfJYVGWOMMcb4LCsyxhhjjPFZPl1kZMNrtlPZGB8i\n39iLHo3xFVXP5PpEXn22yHzz8XNKaCbyfqFPfNDGBLoPGaBMBzleYpk1xum+Fi0ZXZv+JDk+rz5b\nZKrXrMKGltEQO8vbj2KM+R9oLJ/iWi9QUsXbj2KM+TcS7gLXfCGC095+lH/Lp++RObhzu5Zcs5G2\nO3Pg8KPokIZ2QZYxTvaQqCw8B/zECb2DKE5bZo1xqAQRTRymMD8NbXY1HHTmm+Z9dkYG4Or2HeXo\np0dImHOA6zq9gWx82XdbmTGBYIEKwQIDO9NA5pIgYpk1xqESVIVNwMJmSNo2x+bVp2dk/kBf/+sD\nLB5Znx0tZziyMRpj/kOCiN6kcdz+2kZ0jL292hgnSxDRhMFQaXoBFZHVHZdXn56RueCBw2NpO7UD\nz35VwqTfEtQ2ABvjbImtlNtlJCyHLBGVz+w0kzFOlThSkaSNnBtWnQQRvZtljsqrXxSZT66aKvOz\ns7nxqWtpuvA3/tbrDeSL7x31QRtj/oPuQYgdgN4ifKyPQzJIQyszxjiRzkQYGodrtZCwE77OuhO5\nxzl59YsiAxAMRIWE0H1AM3jFzdk+v1B/1xOO+aCNMf/i0HYSpkGezIYwYARIT+cMjsaYP1ju+VtC\neziXWZ2wpJM4ZSbVX/bI0GHfKI2vUYPQoCD6RESQOmUPhXkhnB7bgJePn0Sve8dx63rGBDIRdJxO\n4h/cyAapBdQGmkF3mLHuMZ5mrmXWGIeQBujfT9zMNjoxXobDsJbQDigBbSFwp/dONPlNkQFPmWke\nGkqwCBMbNKDOqQoWTP6ZNlOaM/nECTZuegl9xHkblYwJVPIkyqx06BQD5RCanEtJdi1C656ipHtt\n2JaOaoxl1hgHkA9RRgI9gGiY8eZj7KQd13CQ8U++CbMUVbnsefWrInNB70MvaN3gYKJCQhhRrx5n\ntuWQ36YG87OzSc7Px71uCqXPhtngaIxDiGxXaAJxdWA+0B2IBZLnQehwKActs9NNxjiBCEookACD\nx8wj6Z7hMF0hdgd06gjRoF9cvrz6zR6ZP1oZO02yy8vZUljIHWlpXHl9PU5+coS5DRoy+coradVj\nIrJqldYdfdD/WpwxPki1o9CuDgwDmpeBezbMgk91JZQUQflS5AZUJNkya4yXqSL0gBvGfEdSzHBY\n8Rm3N12JbukE24AV+5DnUZGcy5JXv5yRuaDqa6e0cn4++VdcQe3Dh6k9dBbv7KnMrQObMdrtJik3\nl5NbZnHlT1txv3+9fdszxsvkNZQIYEQeZ6rU5KOSR3jynvmQcf4HUt8GngLOD6bGGK+RDSixJZAc\nimuwkDAHpJ16NgafBuavgdE9IQn0xKXLq18XmQv+kj5OQ4OCmJf8FIQdIvO6LeTlnqVp2zr0P3yY\nL34LI/ybB6mSl8fvb7eywdEYLxPZp4xtiXYTZI/CaAg9nUvJyNqQBLQCDgARoCet0BjjTSJlyvIQ\nFvQXnit3k3vgSj65ti8PPL4CkoEueL6MuEHTLn5e/XJp6V99FDNFokNCeKjrW9TdVJsXNZf9TYI5\nlVnE8weD+blTODc+OJOTt6UT/chP/t/sjHE41ZbCapAnFEavJE7XUTK2NoVfCDfodzALKF8J7S7s\nrzHGeItqiLAcHhqs5AbnsOLaO3jgxxW45gjv7R9C1ILDkJwC3UEk/aLnNfhi/x861aQrEmT8cZfW\nSUtj3qErWPr93dS8fxrHjs7k6LHXGVSnNmvrf8fxP42lweM79cTs9vYtzxgv0j3ILq7Wg1xDZUop\nnV2Z6p0UJB0mAhTD+nVAbWQDql1tZsYYb9HFCK1FeRpYDIsG9WeIuxgkFLYAHIU524FomnBAM2h+\n0fIaEEtL/0rmqlL3H4QcaUvs2rWcatqUm59fyrLt/aCoEfpgpA2IxjiIyNvKsKegHOIWrGPTh7fB\n0E0wPQ5G56Eabpk1xiFEVirTe0M2JE/pzG056yirC1GaTaZcddH3twVkkQF49MgEbVK5MsfKSskq\nK+fogKaUVqvG7h8G24BojAOJJOn12pCfGnSFTHAhfKt/ZyO3WWaNcRiRNfp3fZXb79kIJeBaK3TT\n6+lKykXPa8AWGYD8U2+rhodQroq7tJQ21UbZgGiMgyXRX4fUWw4FeG4UtZNLxjjX16IyW8ENZICe\nujR5Degic8GOohnawUqMMT5DBLUSY4xvkBao7rfj18YYY4wx/0VAHL82xhhjjH+yImOMMcYYn2VF\nxhjzv9q7+zAf67z/48/3mGEwM+6GGTVCbqIkRZuKK5VaSiKpFJu2+1KpVJSamahou1GpdGMbm02K\nkmVpVTstrigkucuIwZTBDIPBmBnz/v3xZdtrr926rt9lfM/vd16P43AcGR35HOdxvD69zs/5+Zyn\niEjEiugis/Szx7z+8ju0yUckAmSYudmryqtIhAhldnrgMxvRRaYwvyUpcXHhHoaI/E8Nvz3cIxCR\n/4WGfma4h/CLIv7U0hdzHvHzUgopPWOCjmKKBJzZfqddLfxbHZ0WiQTWHKcX+EvBzWxEr8gANEpL\n4OUZZegRk0gk+BLahd4DE+6RiMj/QO5+6AHWJLiZjfgi0+LUEXbCSXV5fUErfpv7cGAvtIiAezdL\nf9egF5jlKa8iAedey9J7GbwLVjOYZSbiiwxAy/bJtDptN6e9kMctmx4J5IUWkZCMH8CbGuk0wZ4I\n5sQoIj/JGA/pXYz0EsM2lQUus1FRZE48pT4TmpbRefTp3Ls9ng92jXWb90bgLraIgD3v2D4nc4LD\ncrBM3GyF8ioSUNWuLCYz21nnE2FGHIn7dniQVmeioshcvymXgxXOy9u3U6N9XbY9vZbdF5RjCx8N\nzIUWkcOeAbLmwXhgGvgq4zTfH9hla5GqruK22jABpnT9LYyD+xMa8emBc7CLgpHZiD+1dMQ5a+/z\nTrVqkRwby8AGDfj6zXWsW76D9+9NZdnJ4wK721qkKqpeuMfL+iTBlcDQpTC3I6N+PYxH7QHcU5RX\nkYB5inv9NW5l06o20K6MVr6aN7iFbomL8b3hPdEUNUXmiGs3DHeA8SecwN71e3n3+eU0eeo0rluz\nDz93lCZIkQAxW+xceRasB1oCcyG92MjkI9x7K68iAWO2zXk3hbOv/owtNCHvolakf2Jk8g3u7cOS\n2agrMgC3bHrESyoq6J6UxHk7Ypjx2iouG3Mm12zYwJdtn9fkKBIgobf9Xgdswz9oTYe+X/DNoM4w\neUXYJkYR+ffMcDoBzSB9mtHcr2LwiKkwBtyP/epMVBYZgKs2DPeEmBjK3RnYoAFxf9lOi94nMG3X\nLsZt28bm9uM1QYoEhDXESQOWl0GnOCgHlmfB5MEwGLwsuC/jEqmKrCdOG2AcMBxYCcxaTZrHkde8\nFb7x2GU2Kjb7/ivvnTjGyt1ZUFxMrzdv44Ljjey9e2nxUQEL2rThgnXD3BbfG50tTiTC+A6MunCi\nr4NhwBCg+2AYDDQDs6luNlZ5FQkIn4NRDCP8MXo+9QF1puXD4JPJO7MV9AGzbDfLPiaZjdoVmSOu\n2jDcy935cHUX4ld3oax2baZXf4BTr2/J8v37Gbd9O/PXXoxf/mvd8YmEmT0aOgXhTY2a/QspGVof\n5gL5U6HT1bBkMXBWWJavReS/sz/iJMArfYxp/ic++7gXZAHvTocr+8E0h+GGP1V5mY3aFZkj3jtx\njBWUlzPglIVc1n8s8UVF9DnpOGrHxHD6j87b3ogxZy2m0Td3es1lt0V3qxMJOB+FMXo69rlTUrcM\nslZAN/iVp0IyEH8WAGarnRRTXkXCzK/D6LOfO2Y4n9ml0KMQOsP9vgbSgC4Gk0PfWXvDKiezUV9k\nAP520rP2zuFHTfuaHiR5Rm9SZ/yKN+oe4NEaRXT7ppTRxx1Hw9hYYl4p1+QoEkbu/czfxmAxWf4E\ndbLy+fKi8yAZfnXgc0gAiMO2OxYXjPdYiFRl7rXML8cgl+KSRoy6ZxjPjhgJHaDX/PehO8ABbmnp\nWOOjn9mof7T0zzp1et3zT+tA/O4iKmJi2Hv88RQOHMdTxx9PVmEhG/88ipIHkrRsLRIEP5pzM1T7\nfTEVRbWZf1JHumxZhl3isLIMiNNjJpEguc084zXI2AC3Nh/Hax8OZc8AqHOZwywg4fCeuKOoyhUZ\nAHtrl/c/byypsbG8UVBA3IyRtP7zn/nhzDPZ+urpmhRFAiT0pew8YCeMbw8tIb1HKKYZ7sqrSMBY\nTZwSgG3ULy8H4K7YNKByMlsliwyE3jUzuEEDTi2LZR4HyDl4kIdSH9OkKBJAVg+naAWw8/BPzsBd\nK6ciQWX9caaVAblALNAMd6uUzFaJPTL/yutNn7AZRUW8crCIZjVqqMSIBJjvwujVHsjB7zlfJUYk\n4Px9jAlxwDLSObHSSgxU4RWZI57bNsrvS3lUk6JIBKi5e6cfqFNfeRWJEDcx3t9kSKVmtsoXGRER\nEYlcVfbRkoiIiEQ+FRkRERGJWCoyIiIiErEivsi06TXT7a1d2ugjEgk6m2dU0mvKReToyrDIyGtE\nF5kpz93u3z30DQ2//iHcQxGRX2AjcJtaSiaLwj0UEfkFdQ/m+9u+kkx2h3sovyiii8z+vWUQn09J\n3brhHoqI/IL6o3+AcXHQ5qxwD0VEfsFdNV5kw8enQJ+kcA/lF0X88evnto3y+1ek4hfdrHdLiASc\nWY6H3vAZp7yKBJwZThr4lmB/zyyiV2QAev14CYNabsBmT4/sRiZSBfj41sBUrKe+Wi0SdOkY5Dn2\nbLDzGvFFpvXpHa3vx/tYlrCCmKW3un08KdAXXKRKu9ON+IGkzzWsJhGxkVCkqspwNzoZ6cMMOz24\neY34R0tHbM19xqe/upKnr0ugaNpw9jyeFuilMJGqrHrhHj9UXo2K1NrACtzbK68iAfUcd3hrvuOy\nxp9Cvlfqd5P+f0T8iswRX9dzznzsNO6asJO6V46h0d0ro6OhiUShsmZJVKTW4lm/A1LbY7ZNeRUJ\nqPsHvMJlTT6lYpXBMMMsWCsgUVFk7t480qft2kVWYQF3PHUOd4zfS9ptb1JvxMZAXWwRAbONTmeA\nZdw/6BV6bX0fuqVglq28igSM1cTpAORBzBznnd/1gTGG2f7A5DUqisyiffvIKy2lWfUafBdXzoD7\nTua2KTW4665J2KRt3qrvnMBccJEqr3NzWASM6QiT5zFren/IzgPOxc7E+zBFeRUJimHA8sO/BuZw\n7XczIAOgJvYAvooWYc9r1OyRAbhqw3BPiImhe1ISXRISyH52J70eSKbbunVseuc+do8+IVDP9USq\nMrNCh09gydXQGVjrjG1xNw8NeAlSwZ8P9pFPkarEDIcV3OJ/5fWu9/DK/MHc/tkkbI5DO/Drw5fX\nqCoyAL3XP+jJsbF0rjvAX/sAABlcSURBVF2bdjVr8umQRTQ/pT6rr03hqQVX0HBhNerm5rJuZi9N\nkiJhZrbCIRG6NcdHG7beufz6KXzUdgCkAQsgvcRCpydEJKzM8hySYFwSux8yWpZs4is60ezs7Ydv\nRiB97rHPa9QVGYCL1w3z5NhYiisq6JaYyD0NGmLVjMk7dzKjqIgPvxxMkzmFbP5DF02OImFmTXDy\nsoA9nO3t+KLVBbAeSAXqAmsBNuLeXHkVCTMbhDM5C2jJFz6Es19YDlmEshoPrATywP3YrdBExR6Z\nf/aX1s9YiTuxZizZt48z161lV/5+Et/fyjNpadxw9ttsGbCFWk8W+HG3LIm+JicSQXwLxrTBwN18\nYXVg/VKu8knQicMlZilkNMcScesa7BdziUQ7fxtr5adDuy6cXXM5DHU+/foc6EXoBiRvDywB64vb\nwmOT16hckTmi1pMFHlNRQXxREQcTE6mxdy+fPf5nVry8lorfNmXJvn289PmtpP5nAXW2bGHtrN66\n4xMJE4vD6QEUw1vZxg28Ahm3wwQgfw/EJ0EPIBsoOrZ3fCLyX9kTOO2ARZA+xsjkXVh5NVxDaFWm\nC6Fy8wmQC55TeXmN6iIDoQ3A7504xo7882ev/JZ1Yz7Bdxxk/94y5tUvZ/TWrWz57GnKO7xIy3Ht\nyfmwpyZIkTAwy3b3bhb658UOyxjnyxk66DWYnANkw7CbQ+UmDVhbpu82iYSJ2TZ3Tzmc1zKHqfjs\nQdhWh5sApv5Ubq4E1odWdI76OKK9yPyjazcM9y8e+TU7W7Sg429eZPZxzfl02npOuKop2cXF3DPn\nZhK/L9NbgUUCwBriFAC3AROWMsqncC1/pIU9BfQEUgCtzIgEgQ3CmQW8CVy5h9KiOsR1AVv5V6Ab\ndKfSVmaico/MvxNrRp0tW2g1dy5bHrqQk28+mWvyRpBfXk7uwYOQOpe9LeLCPUwRgdCqywzIePUh\nzvbdTONKWizdCt0GH/4X9gOrwzc+EfnJmBLq5Ofjc42RPoaX6txOi29XQkY3wOGTMli/p1L+6iq1\nIgMQ80q5183NpeUNL/PV5g7UXHopJ3/wAUvv7wp72uC3BusbEiJV2U2M94kDhsC7rxK6pUsinVQy\nmoBtWYd7K+VVJChuMLfODrdNBc6F2DTSy0MRzWQf7rUqJa9VrsgA/GbjCB+aksLp8TV5fsd2Sioq\neGRtEn7+vZoURQLGHsDrjM5nd59UmPtG6IdZN4f1BVwi8q/ZbPyWS1/g9QfugWemArHEFVxEaYMk\nbfatDC9uH+1dEhI4qSyW2nWGaFIUCTA7HWckUBf8QpUYkSCzF/C0e3LYezCRohqplZrXKl1kAB7/\nMcMfOy5Dk6JIBDB70d3vVl5FIoBZjh+Lx79VvsiIiIhI5KpSp5ZEREQkuqjIiIiISMRSkREREZGI\nFdFF5vbNjzhAvREbtdFHJAKYrXaAVxmszIoEnFnoo48L6BjovEZskflmfoZf+tl+7Ks7fddTzXWK\nQSTgrCfOuydzNZP8drKUWZEAs4l4U18LfzbvwtJA5zVii0xRwQEO7C2jXpwTO7400G1RROD+OaPh\nmtW8d+b1oS9di0hgrbuxCZuatMG+cqxJsPMa6cevfeq4O7nm9BPw8x4KdGMUEeBa84wpkDnU8ef1\nUjuRIMswc4DMyY5fF9y8RuyKzGHWpmNDXv90HdWX3ea1niyI6FYmEvXecctoB+883wezpW63B/tO\nT6Qqy3A3AP/aMHOfxYWBzGukr8gAMG/KMF/08WbmP3wC8xY+gt9QL7DNUUTA6uG8CYwO/d6/Du7d\nnkhVZ33xuDf3ULYoCT4hcKupUVFkAH7/+EhPbdqMl875jrnZT+I3xwXqQovIf/X3ZWtKYXgc/lSw\nJkcR+cnf81rXYXkZ3jQ4/4+NmiJzy6ZH/Jp69ch5ejWf/7YR75w4JjAXWUT+O7OlDimM81EMtYtx\n76fMigSU2TanUwr+sGFXbMA9OKeFo6LIXLVhuMdb6JqOPv54/njnAlKeas+NCzvil/86MBdbREKs\nOU4bIA14sxBvlYx1dIgHf0srMyJBYufj9ABigeGwpdxoMtkhGfzX4c9rpG/2BeC9E8fYguJikmNj\n+aK4mH4vncO+Md/xhzUf0Hv9g54wenvktzWRKOIbMW9hkAf0aoDlfMMrUwbjdxo2Ebd62gQsEhT+\nVyx9uEEBkA1NcLZfl8i0HoZ9g9v54c1rVKzIHHHF9w85QLwZk5s3Z1f+fsY/+J80fq4Dty44B+97\nftibo4j8xGyb0zIF1m+jqe9i08dtKO1sVK+bD7EpeFn47/ZEJMQMZxwwtIwB/jZbaMLs+Iupc9Ch\nB/ic8OQ1qooMhMpMQkwMJe50S0jg9uSGLPp4MzXOa8gbBTuYMO9hSFqN9z9TE6RIAJjtca5Mgmmr\ngTSYlgQTgE9++nfSsb8fBRWR8LHTcfoAGUACxKzfR8WC2qETiHWB9ZCed2zzGhWPlv7RBy3GWn5Z\nGRmNG5NdXMxZ360l5j+SWf3yWu5tlMLEy1+gyUmTqPVkgbfpNTO6WpxIBHJPMqatALYD82Ao0O3w\nL14EppL5jGOGm+UosyJh5F9jLAIog+IyKobUDv3BcCB7MeRtJPMTxxrjZnuOSV6jbkXmiNjxpV4R\nE0Pijz9SrbSUstq1uWjgOAbNO8CFN7fl3Z07uXXJiSQtPJdGq1ZRZ8sWliy5RXd8ImHSiM3+NaeT\nZq/Am1dDLtCd0MrMDKATsBZYBCP9YUa3fRJfo0dPIuHwKoO9u02iNYs4z3fx7aFTOanad3wx/QJY\nQiivecAY8FcMe87x+ZWT16gtMgA9cx7wmjExxJtRDuQePMjKDx9h9/3vU1RwgF0NYskqLOT5yUOp\nvWMH1ffs4cfXO2liFAkTewFnFrAcGAID099gst3ELf4ijW0omWyA1OaQXwi8BJyqY9siYWKf4+QS\nutkYAkvPOpmOcauZXnYJK2wOmTgMAbKA4nlATdy7HPW8RnWR+Ue/zX3Yiysq+Ntz17LtwkJmX7CE\nnolJrFi4lUMd6zJu2zbemZ3JcV9+yeY/HP0LLSL/OzYR56YVQFsYHBf6YR7wyXTgCkg2KNgPTMD9\nPmVWJIzsG5wODtcYDAYSCJ1y6pMDHVpBB2AyUL4n9Dj5KIq6PTL/zu+bPWklFRXEFxXR8YkcBj7R\nlwYrvqHDugeoW60aHWrV4lCTv5LXuXO4hyoigN+I0aw9ZMRB1lLIWgFvlkHdfrzt/aFgJrAGut0X\n7qGKVHl+GsY4g5FAj0LoUka7y7+Cga3wTIOsHCjfD8OSjvrfXWVWZI6oN2Kjl9Wuzb5Tv4aiDjT5\nNIdqZWUcHDGdrYuf16cNRALG+uMsAFIJ7ZEZCHQG4oGBL+J+tzIrEhD2As4nENrPAQyF+Gt2Eht7\niOKE5bhfpEdLR0Pv9Q96cmwsw1JSeHfXLpbv38+3j/Zk4xS9Z0YkiMxedbrczo3zxzPR+kObFOos\nz6eoRqoyKxIwZoudYWdROtKonhraJ/Pk7+5lBM9rs+/R9MGusV7mTtv4eJofjCGxnu7qRILMMvHv\n0xvTIm4rlE/XJl+RAKu2bZ8f6peA7XZYWYZ75T3tqLJFBuDuzSM9NS6Ohxuna0IUiQBmGx2a465j\n1yJBZxb6plplvyahShcZERERiWxV5tSSiIiIRB8VGREREYlYKjIiIiISsSK6yHwx5xHfWf6K28R9\nnvRYnjb7iATYkY8+2ljcGqK8igScNcHNyrwFq9xODW5mI7rInN3zCftu3g807PQg1ffvD/dwRORn\nhE4aLYRFQMHScA9HRH6Bb8FgMRtePgVW5oV7OP9WRBcZgLN7PkG/enUpPL0abXrNDGxjFBEg9zqY\nBTAXs0LlVSTghvh/hr6RBFjNYK7KRHyRAaxJ5iYatnuODRdcEO6xiMjP8KZxll5uQFPo1iDcwxGR\nX/ASD1r6IgOWwdBwj+Zfi4r3yGxaO9aznljCRyOO5+uNXWk2ub4+NyASUGbTHVKA9UBboD7urZRX\nkYAyy3GIBZYBrYBE3JsHJrNRUWQAtm1+1vM37WVp2+rcOP9svK+KjEhQZZh5JgVAA2Ay7gOVV5EA\nC2XWIQ3Im4N7z8BkNhoeLQFQVHCA2S2N62rXZdJ/LMI++jg6GppIFMpMdkIrMquBqw+v0ohIUGX2\ncKAQ8sqgTU/MFgQms1FTZF5OLmFtSQmrvJQzcsrY1DMHmz8qMBdaRP5BQRZwPLAdUuMgtR9m+5VX\nkaCauxQSGgA74SagTxeCciQ7KorMbzaO8NzSUt5efyLP5Oezp10CvmEfGa+v4NoNwz152LpAXGwR\nATN36E5ojToH8l+FfIBamM10szeUV5EAsZo4CR0hHmAbDFsAM4BiMMt2s6lhzWxU7JFpvfJur1ut\nGgANYmOpV60a3RITGdigAXNfK6X4qn1c/1l3mn1QrE3AImFmL+OsBCZMB/YCp0KPjjDXoaXB+sWE\nNgHXxD1OeRUJs5q7d3rJovrQYwGwDWgL15wMyUAxkJUHCWlQTFi+TB8VRQbgjNVDvVmNGuQePEha\n9eokxMQQHxPDxBOaMnfydxx/5Qk8s20b78zOpPHSpWx562xNkCJhYjbTIQ569IS5q4EPgXOhRzeY\nOx04wOFdhZA6EN967CdHEfmJ2WKHlpDdALoBPAf0gMknw8BCYCFwKrAW2vXEvz12mY2aIgNQ/bli\n73X5KACKDx2ibmwssUCXhARazNtF2/7NmVFUxD3fOYnzbmDP42nW5IYvXKVG5NizTJwMB9Yc/nUA\nSKKdN2ZlqzNh/ZFC0xZYEzrZlGXOYFdeRcLApuPcBJQAJTnACqAt6/zXtO65BeauAHKB3sBM3Hvb\najM/2Ss3s1GxR+aI0vsS7MMl1xJrRnJsLO9//hCzsobwyPPXcdWae2j69feckr2bH86qx+ghWdiU\nNV5Wu3a4hy1SJXk65o/FANshrR8sHwjEsdLOhPVjoUs/4AxCJQfMtrm9HD03XiKRxvth6UUGJauh\nZSsu9lgAWtfbAnNfheHtOVJiIAUz/JThlZ/ZqFqROeKqDcN976FD5JaW0jY+ngMVFWTv3ctldesy\n691h7L7/fYoKDrCjfjUytm5l5vSHidu/j73pjXWnJ3KMma1wWAO9roZZoXehpzOITN4GBsIw4Jl5\nwA/AuUAZ8CXug5VXkTCwejhFM2Fwb8iaCdQknYvJZBqk9oN3gW57gDmQfDWUA0VTcb+6UjIblUUG\n4NoNw72c0LsIp2xsQt9mm0mIiWHWqzey69QafHr5VDqWxlF68BDzapZy3deJJCy4ggY5OeRO/g9N\nkCLHkMXhlGcBB/DMO7D0VYf/JInQ46UUIIfXfRy32IWQPBgKngCOV6ERCQNLxCmeDpTzCtdwB+uA\n+hDbABI4/PhpNXnelTR7Hdr1g5XTAXDvd1QzG7VFBsAWPuqXpRykzJ3GcXG89cUgOj47n9Jatf7+\ntexLZv5IcmwsD+XlUfL1s/iNtTUpioSBvYwzJJvQ46QcuKkjzAWugTqj89m9NhU6FALZQD/gVdxv\nV15FwsS+welQCKkNID8PstP+/mcjznuM+XRlgZ0LrIcO7WF5VqXceER1kTni2g3DPaFaNb66piXf\nXXYZB7q+So/UPczNO47GU7py3LJlbLjgAnY+3VKTokiYma1w2rUnfaWRyV8hvhuUTIZmAyF3MU29\nDpusRqC+9SJSlZltc7qnkP6JkUk+dEqBJdPhmX4wDIb404y363FP0aOl/6uE0ds9prSUPje8TIda\ntbi3YSM+3F1E0aFD3LjmEN71UU2MIgFhs3F6vUjo8VIZ0JPQaswB4CbcTXkVCRD73p2WbwL1CeW0\nN7AY2APx/fADlXMkO6pOLf2S4pGNzKtVIz4mhuJDh3ijsACA6+vWV4kRCRi/FEvnHuAy4Axol0Zo\n1+AZKjEiAeQtzNK5hVCB6QndkoA9wFmVVmKgiq3I/KP7tjzqsWa0rFGDWxo+oklRJMDMtjnMI/S2\n36O7UVBEjj6zMg+toIL7RXqPTGV4rskoyy8rY/3Bg+Eeioj8gtCz9QOEHi+JSNCFPi+yE5IvqvS/\nq8quyIiIiEjkq7IrMiIiIhL5VGREREQkYqnIiIiISMRSkREREZGIFbFF5qtPHvUfvv+d28eTtFtZ\nJOBsIr6KFj6PrsqrSNB9Zr7azPnMIiKvEVtklmX/wJJP86D+l9QcuysiLrZIlXXTUtot/J6Lv5uP\nPYDyKhJg1s+50L+n2wVzsInBz2vEFplbR//eEupW59OTmlPSYWa4hyMiP8O9o3EN0GY1LAn3aETk\n5/guLP/mE/ncLgz3UP5HIrbIAFx41e/s0Ox8JnX6Eft4ktcbsTHwzVGkqvItGCyEbLC++A4SlVeR\ngPI3MFgG48EW4jwZ3MdMUfFCvLfH3OpFBSU8f30S2z96lOKRjfQKc5GAMnvD07wbeXGtoA34t5X3\nDRYR+b8xm+M3+vdMbDIEehwpOMES0SsyRwwa/pqd3fM67p5Yzhn9x4Z7OCLyM9xvtinWmoll10LL\ncI9GRH6Oe09Ls7vwoQZdwj2afy0qigzAB20WcsGTJ3Hh6M1cu2F45C8ziUSxri2dG+0S/vJhV+yy\n4G8mFKnKMgc6NmwR/q1hjwYvr1FRZB7/McNLKirILi7mpoxfccEfYhib/7g3untl4C64SFVnttjJ\nAziLixfOhyVgttStbfAmSJGqzi7CKQE4lUbPbAr9zHAbEZy8RkWReey4DCs6dIii8nLyUmMZcO9x\nVBu1mWlPfIzNmRKYiy0i4H6W1S/+gXRaQ5c9+DSDmzpC0eFJU0QCw+dhGe8/RDq12ZF4AhtGG2me\nA4vAXghGXqNis+8Rv9k4wuNjYkiLi2NggwbMeGAltZPKWHpzY974SwZ+c1zgNimJVFX2DU53YBbQ\nGXydMa9VFy5+YD7M1SZgkSB5gVt96Muv0fbOZaw5/wzWZRuttoFtdZgM/rvw5TWqigzALZse8WoG\nDWPjuLROHVhQyKfv5ZL29Clcv6g19T4/jZ1Pt9QEKRIA+xPMax90BpT9ngmHbqNewS4qFtWGPvuB\nD4Ak3HsrryIBkGHmme2c9769jP63zWLQhNdZTgdWJp4J5UDJ0tA7o46xqCsyECozCTEx5JWVkRob\nywtpTVi7dDt7Tk4ge+9ehn9+OcmLqrPj+baaIEXCzGyBwxrgZmAmb/sfaMZGutpSYCPQHJaAd9QK\njUi4meU4tIJmQO4K/K3TKOkPNVs65O+BTklc/tUUZjDgmOU1KvbI/LPXmz5huaWllFRUUHToEB3X\nrqFl+2TysjZwe8OGLLz0U3YPeoH43+3xZgP/5q17z4q+NicSIdy7GItuBnKAFAa9No0UtsNggANA\nBkwGa4i/wK1uicF4Li9SFbm3sqt8Uug3XdrTbPAa4l/gcF5jYclUPlo1AOuKb6KRW9fKz2tUrsgc\nUWfkZi9NTCRu3z4qqlfHDh2i729e5FcTtnLJ4x3J3ruXG+efTeLXrUjIz+fH1zvpjk8kTLyBeczv\nHd4kdLc3fjrk9obRcaGf4cAP0CEttKfmVa3QiIRLhpmPyi+mevxBUupsZ1PfNlzw4Sw+m90LbgIS\ngBKgQ+iXj6q8vEbliswRu0efYIMHPEOfG17mkmufo//1L/HO7EweG1CL9QcP0r9aAoWXrWH47RPZ\neuZpJIze7vaae6u+c6K33YkElBW6AYz402OQDXAFNJsKawGyGeUPQJc0GAckgFmZX80kN8tTXkWO\nsQx3q5WwnwN/asCmP7aBGfCZ1Qh9Sy0/D3/NoBd0+dM8GBg6sr2DRDc7+is0Ub0i84/u3jzSE6pV\nI6+0lBlv3Un3619iWEoKJ+0NlcTCetWYtmsXj8x8jIarVrH9xXa62xMJk2Wc7GdcsQbb4jAShlz+\nNOPHPhj6wwLgmW3E5CdQ0aU2pILP1+qMSNicap6xEjJvcuKf2cnwOmPJWDWWOi3zqRZ7iJ2xx/Og\nZ/L0oHRIA3/q6Oa1yhSZI5KHrfMT/vY3tnbsyKG4OAquf4nuiYncmJxMYXk5dy48F4o64DfU08Qo\nEmZ2GZ4+y8gkH9qkwAQgmdCR7eErgMXAAdzvVl5Fwsym4unXGJk4dIC4T/ZwSYPZfDR7APQCyAJq\n4n61iszRYFPWOPW/pN6n59IsO5vV/a/iYOcXYN19+I21NSmKBIg1x8kF2oGPNS65ZDpzrDsw56hP\niiLyf2Pn46wHmkH6AqO996Rf4p+heA7uPY96XqtskTniwbzHHGDCjh2ceNfZHExMZM2cvpoYRQLI\nDKcXkAzpWUbGeUC2K68iAWQNcYYCeZA+IRTTDD/6ea3yRQagwl/zHeXl5JaWclbtezUpigRYhpnH\n+v08ao/jXkt5FQmwDDO/rzyOOsml+K7K2csWWxn/0UgTY7daShzEbH/OqR3u0YjIzzlyR9cWc/RK\nGZFAO5LX9KLKy6tWZERERCRiRfV7ZERERCS6qciIiIhIxFKRERERkYilIiMiIiIRS0VGREREIpaK\njIiIiEQsFRkRERGJWCoyIiIiErFUZERERCRiqciIiIhIxFKRERERkYilIiMiIiIRS0VGREREIpaK\njIiIiEQsFRkRERGJWCoyIiIiErFUZERERCRiqciIiIhIxFKRERERkYilIiMiIiIRS0VGREREIpaK\njIiIiEQsFRkRERGJWP8P5v9aaPANY9YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x285ce27eb38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dx = .1\n",
    "dy = .1\n",
    "isub = 0\n",
    "# declare a plot size\n",
    "plt.subplots(figsize=(12,6))\n",
    "\n",
    "for i in range(nf) :\n",
    "    fn = bnames[i]+\"_DEM_10cm_trimmed.tif\"\n",
    "    fn_slope = bnames[i]+\"_slope.tif\"\n",
    "    fn_aspect = bnames[i]+\"_aspect.tif\"\n",
    "    fn_tri = bnames[i]+\"_tri.tif\"\n",
    "    fn_range = bnames[i]+\"_range.tif\"\n",
    "    print( fn )\n",
    "    \n",
    "    tstart = time.time()\n",
    "    with rasterio.open(fdir_in+fn) as src:\n",
    "        z = src.read(1)\n",
    "        z[z <= -99]=np.nan\n",
    "        ysize, xsize = np.shape(z)\n",
    "        profile = {'driver': 'GTiff', 'height': ysize, 'width': xsize, 'count': 1,\\\n",
    "           'dtype': 'float32', 'crs':src.crs, 'transform': src.transform}\n",
    "        print(profile)\n",
    "        src.close()\n",
    "    \n",
    "    s=scipy.ndimage.filters.generic_filter(z,slope_filter,size=3,mode='nearest',extra_arguments=(dx,dy))\n",
    "    az = scipy.ndimage.filters.generic_filter(z,slope_dir_filter,size=3,mode='nearest',extra_arguments=(dx,dy))\n",
    "    tri = scipy.ndimage.filters.generic_filter(z,tri_filter,size=3,mode='nearest')\n",
    "    ptp = scipy.ndimage.filters.generic_filter(z,np.ptp,size=3,mode='nearest')\n",
    "    \n",
    "    # write the tiff files\n",
    "    with rasterio.open(fdir_out+fn_slope, 'w', **profile) as dst:\n",
    "        dst.write(s,1)\n",
    "        dst.close()\n",
    "    \n",
    "    with rasterio.open(fdir_out+fn_aspect, 'w', **profile) as dst:\n",
    "        dst.write(az,1)\n",
    "        dst.close()\n",
    "    \n",
    "    with rasterio.open(fdir_out+fn_tri, 'w', **profile) as dst:\n",
    "        dst.write(tri,1)\n",
    "        dst.close()\n",
    "        \n",
    "    with rasterio.open(fdir_out+fn_range, 'w', **profile) as dst:\n",
    "        dst.write(ptp,1)\n",
    "        dst.close()\n",
    "  \n",
    "    print(\"Elapsed time:\",time.time()-tstart)\n",
    "    \n",
    "    # small plots of elevation\n",
    "    plt.subplot(nf,3,isub+1);\n",
    "    plt.axis('off');plt.imshow(z,cmap=plt.cm.terrain,vmin=-2,vmax=8)\n",
    "    plt.subplot(nf,3,isub+2);\n",
    "    plt.axis('off');plt.imshow(s,cmap=plt.cm.jet,vmin=0.,vmax=.5)\n",
    "    plt.subplot(nf,3,isub+3);\n",
    "    plt.axis('off');plt.imshow(ptp,cmap=plt.cm.jet,vmin=0.,vmax=.1)\n",
    "                   \n",
    "    isub += 3\n",
    "    \n",
    "    print(\"Elapsed time:\",time.time()-tstart)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}