Sandwich_maps_local_coords_v2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sandwich_maps_local_coords_v2\n",
    "\n",
    "Read DEMS and convert to rotated alongshore/cross-shore coordinate systems.\n",
    "Modified to include some of the masked channels.\n",
    "\n",
    "This produces the 1-m gridded 3D elevation map saved as an .nc file. It has been run twice,\n",
    "once with a 3x3 median filter on the maps before interpolation, and once with no filter. Resulting files are:\n",
    "\n",
    "`Sandwich_1m_DEM_v2.nc`\n",
    "\n",
    "`Sandwich_1m_DEM_v2_nf.nc`\n",
    "\n",
    "This version differs from the original because x is now the alongshore coordinate (positive to SE), and y is the alonshore ccoordinate (positive to NE). This happens after the maps have been interpolated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\csherwood\\AppData\\Local\\Continuum\\miniconda3\\envs\\IOOS\\lib\\site-packages\\xarray\\core\\merge.py:17: FutureWarning: The Panel class is removed from pandas. Accessing it from the top-level namespace will also be removed in the next version\n",
      "  PANDAS_TYPES = (pd.Series, pd.DataFrame, pd.Panel)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from IPython.display import display, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate, signal\n",
    "\n",
    "import xarray as xr\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fdir = \"D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/\"\n",
    "\n",
    "#fdir = \"C:/crs/proj/2017_Karen_Sandwich_maps\"\n",
    "\n",
    "fnames = (\\\n",
    "\"2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif\",\\\n",
    "\"2016-01-25_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2016-02-11_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2016-03-30_SandwichTNB_AS_DEM_10cm_trimmed.tif\",\\\n",
    "\"2016-09-21_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2017-01-09_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2017-01-25_SandwichTNB_DEM_10cm_trimmed_channel_mask.tif\",\\\n",
    "\"2017-02-14_SandwichTNB_DEM_10cm_trimmed_channel_mask.tif\",\\\n",
    "\"2017-03-16_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2017-04-28_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2017-05-04_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "\"2017-09-18_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "#\"2018-01-03_SandwichTNB_DEM_10cm_trimmed.tif\",\\ # needs to be padded...not very good control\n",
    "\"2018-01-10_SandwichTNB_DEM_10cm_trimmed.tif\",\\\n",
    "#\"DEM_10cm_selfcal_lokicampos_mid_denseNAD83UTM19.tif\",\\\n",
    "\"2018-03-09_SandwichTNB_PT_v2_DEM_10cm_trimmed.tif\")\n",
    "\n",
    "# list of dates for use as titles / column headers\n",
    "titles = ([\\\n",
    "         \"22-Jan-2016\",\\\n",
    "         \"25-Jan-2016\",\\\n",
    "         \"11-Feb-2016\",\\\n",
    "         \"30-Mar-2016\",\\\n",
    "         \"21-Sep-2016\",\\\n",
    "         \"09-Jan-2017\",\\\n",
    "         \"25-Jan-2017\",\\\n",
    "         \"14-Feb-2017\",\\\n",
    "         \"16-Mar-2017\",\\\n",
    "         \"28-Apr-2017\",\\\n",
    "         \"04-May_2017\",\\\n",
    "         \"18-Sep-2017\",\\\n",
    "#         \"03-Jan-2017\",\\\n",
    "         \"10-Jan-2018\",\\\n",
    "         \"09-Mar-2018\"])\n",
    "\n",
    "nf=len(fnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def box2UTMh(x, y, x0, y0, theta):\n",
    "    '''\n",
    "    2D rotation and translation of x, y\n",
    "    Input:\n",
    "        x, y - row vectors of original coordinates (must be same size)\n",
    "        x0, y0 - Offset (location of x, y = (0,0) in new coordinate system)\n",
    "        theta - Angle of rotation (degrees, CCW from x-axis == Cartesian coorinates)\n",
    "    Returns:\n",
    "        xr, yr - rotated, offset coordinates\n",
    "    '''\n",
    "    thetar = np.radians(theta)\n",
    "    c, s = np.cos(thetar), np.sin(thetar)\n",
    "\n",
    "    # homogenous rotation matrix\n",
    "    Rh = np.array(((c, -s,  0.),\\\n",
    "                   (s,  c,  0.),\\\n",
    "                   (0., 0., 1.)))\n",
    "    # homogenous translation matrix\n",
    "    Th = np.array(((1., 0., x0),\\\n",
    "                   (0., 1., y0),\\\n",
    "                   (0., 0., 1.)))\n",
    "    \n",
    "    # homogenous input x,y\n",
    "    xyh = np.vstack((x,y,np.ones_like(x)))\n",
    "\n",
    "    # perform rotation and translation\n",
    "    xyrh=np.matmul(np.matmul(Th,Rh),xyh)\n",
    "    xr = xyrh[0,:]\n",
    "    yr = xyrh[1,:]\n",
    "    return xr, yr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make small box with local alongshore/cross-shore coordinates that will be our region of interest (ROI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xrot:\n",
      "(280, 1300)\n",
      "[[5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]\n",
      " [5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]\n",
      " [5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]\n",
      " ...\n",
      " [5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]\n",
      " [5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]\n",
      " [5.0000e-01 1.5000e+00 2.5000e+00 ... 1.2975e+03 1.2985e+03 1.2995e+03]]\n",
      "yrot:\n",
      "(280, 1300)\n",
      "[[  0.5   0.5   0.5 ...   0.5   0.5   0.5]\n",
      " [  1.5   1.5   1.5 ...   1.5   1.5   1.5]\n",
      " [  2.5   2.5   2.5 ...   2.5   2.5   2.5]\n",
      " ...\n",
      " [277.5 277.5 277.5 ... 277.5 277.5 277.5]\n",
      " [278.5 278.5 278.5 ... 278.5 278.5 278.5]\n",
      " [279.5 279.5 279.5 ... 279.5 279.5 279.5]]\n"
     ]
    }
   ],
   "source": [
    "# make a small box w/ cell centered grid with spacing dxy\n",
    "dxy = 1.\n",
    "ymax = 280. # cross-shore\n",
    "xmax = 1300.  # alongshore\n",
    "nx = int((1./dxy)*xmax)\n",
    "ny = int((1./dxy)*ymax)\n",
    "\n",
    "xcoords = np.linspace(0.5*dxy,xmax-0.5*dxy,nx)\n",
    "ycoords = np.linspace(0.5*dxy,ymax-0.5*dxy,ny)\n",
    "# Y is ascending\n",
    "xrot, yrot=np.meshgrid(xcoords, ycoords ,sparse=False, indexing='xy')\n",
    "print(\"xrot:\")\n",
    "print(xrot.shape)\n",
    "print(xrot)\n",
    "print(\"yrot:\")\n",
    "print(yrot.shape)\n",
    "print(yrot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert ROI to UTM coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of xrot, yrot:  (280, 1300) (280, 1300)\n",
      "Shape of xu and yu:  (280, 1300) (280, 1300)\n",
      "Corners in rotated coordinates\n",
      "Upper left [  0.5 279.5]  Upper right [1299.5  279.5]\n",
      "Lower left [0.5 0.5]  Lower right [1.2995e+03 5.0000e-01]\n",
      "Corners in UTM coordinates\n",
      "Upper left [ 376418.75924416 4625460.38726746]  Upper right [ 377442.3852131  4624660.64301101]\n",
      "Lower left [ 376246.98969254 4625240.5322672 ]  Lower right [ 377270.61566148 4624440.78801075]\n",
      "Corners: \n",
      "[[ 376246.98969254 4625240.5322672 ]\n",
      " [ 377270.61566148 4624440.78801075]\n",
      " [ 377442.3852131  4624660.64301101]\n",
      " [ 376418.75924416 4625460.38726746]\n",
      " [ 376246.98969254 4625240.5322672 ]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAABgCAYAAAAesrRvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQ10lEQVR4nO3dX4yc5XXH8e+xjb0ONvaaBuTYJjhgrUpyERILQVJVFTRA2yjmppKrJnVVqtygKmmrttBc9QIpaauoqipaoZDiphSEElJopFZGTiSrbQoyaRqMzda2tjFrXNtRjGOnNjb49GLexeNldnfO7Mz77/w+0mpn351ZjmH88PUzz67N3RERERERyW5J1QOIiIiIiNSBwlhEREREBIWxiIiIiAigMBYRERERARTGIiIiIiIALKt6AICrll/tYyvHqx5j6N4es6pHKM3bK6qeoFxLV7xd9QiluWb5+apHKNW1S8+GH7P/5Ys/cvf3jmCc2lq7bqmv3zia/4WcvLh6JF+3rc68OVb1CI2y5Hye/zcPw9JzVU8wfD99Y3rONbsWYTy2cpyP3vE7VY8xdKcmllc9QqnO3HSp6hFKs2rz6apHKNXdN0xWPUJpPjP+7+HHfPj90z8cwSi1tn7jMh7/p/cN/es+cvzOoX/NNttz8Gb0x4j+jb2qP0REXbu/fRtB//bMH8y5ZusoxQiNT16oegQZkbNTa6oeQUbka6c+VvUIaSmKY/YcvLnqERpFURxz7f63WxnFC1EYj1imOF59ONfTKVMc7zoyUfUIpVIcl09RHKMojlEUx2QM4hm1OErRduOTF9Icq1h9eEmqIxVnp9akOVax68hEqiMVUg4FcYyCOE5RHJM5ikE7xjIC2XaOM8m0c6xd49FTFMcoimPGXh1TFAdkPToxmwqmJJmOVECuOM50pAIUxzIciuIYRXGMgjhGQXxZnnqpgWxxnEm2OM5EcTx8iuIYRXGMojhGUXwlhXHJMsVxpl1jyBXHmXaNQXE8LI8cv1NRHLDn4M2K4gAdnYhTFL9bLcplybmLVY8gI5ItjjPJFseyOAriGAVxjII4RueJ51abalm572jVI5Qm064x5IrjTLvGkCuOtWs8OEVxjKI4RlEcoyCeX62KRXHcXopjaQPFcYyOTsQpimMUxTGK4oXlqZUayhbHmWSK40y7xqA47peCOEbniWN0njhGRyf6V7swzrRrDLniONOuMSiOJS9FcYyCOEZBHKMgjqllqWSL40yyxXEmmeJYu8ZzUxTHKIpjFMUxiuK42lZKpjjOtGsMueI4064xKI4z03niGB2diFMUxyiKB7NgoZjZJjP7jpkdMLNXzOxzxfV1Zva8mR0s3o93PeYhMztkZpNmds8ofwFtkS2OM8kWx5nUMY6rWLMVxDEK4hidJ47ReeLF6Wfr7i3g9939Z4HbgQfM7BbgQWC3u28BdhcfU3xuO/BB4F7gETNbOshwmXaNIVccZ9o1hlxxnGnXGGoZx6Wu2YriGEVxjII4RkG8eAvWibsfc/fvFbfPAAeADcA2YGdxt53AfcXtbcBT7v6mu08Bh4DbBh1Qcdxe2eI4k2xxXCdlrdknL65WFAcpimMUxTGK4uEIlYmZ3QjcCrwAXO/ux6CzEAPXFXfbALzW9bDp4trsr/VZM9trZnsvXDo37z83WxxnkimOM+0ag+K4Doa5Zhdf7511+/wb50c1duvoPHGMjk7EKYqHp+8qMbNVwDeAz7v7T+a7a49r/q4L7o+6+1Z337p8ycoF//mZ4jjTrjEojkVGYdhrNly5bo+tVbj0Q0EcoyCO0XnimNUvHWX1S/P3ZF9FYmZX0Vlgn3D3Z4rLx81sffH59cCJ4vo0sKnr4RuB1wNzC/niOJNMcaxd42poza4HRXGMojhGQRyzUBDP6OenUhjwGHDA3b/c9anngB3F7R3As13Xt5vZCjPbDGwBXuxz7nll2jWGXHGcadcYFMcyOnVas7PS0Yk4RXGMojim3ygGWNbHfT4OfAZ42cy+X1z7Y+CLwNNmdj9wBPhVAHd/xcyeBvbT+e7oB9x9aP8FV+47yrkP9Tz+Jg23+vASztx0qeoxZAR2HZng7hsmqx4ji1qt2dkoiGMUxDEK4rhIFEMfYezu/0rvM2gAd83xmIeBh0OTBGSK4/HJC5yaWF71GKXJFMdnp9awavPpqscojeK4HHVcs7NQFMcoimMUxTHRIJ7R2NevMx2ryHSkIptMRypE2kxRHKMojlEUxwwaxdDgMM4mUxzrvHF76byxtI3OE8foR7HFKYpjFhPF0PAwzrRrDIpjaQfFsbSFgjhGQRyjH8UW08+PYutH4+sjWxxnkimOM+0ag+JYmk9RHKMojlEQxwwjiGe0ojwyxXGmXWNQHLeZ4liaSEcn4hTFMYrimGFGMbQkjLPJFseZZItjkSZREMfoPHGMjk7EDTuKoUVhnGnXGHLFcaZdY8gVx9o1lqZQFMcoiGMUxDHDOk/cS6uKI1scZ5ItjjNRHEvdKYpjFMUxiuKYUQXxjNbVRqY4zrRrDLniONOuMSiOpZ50njhGRyfiFMUxo45iaGEYg+JY2iFbHIvUiYI4RkEco/PEMaM8OjFbK8M4m0xxnGnXGHLFsXaNpS4UxTGK4hgFcUxZQTyjtZWRadcYFMdtpjgWKY+iOEZRHKMojik7iqHFYQz54jiTbHGcieJYqqDzxDE6TxyjoxNxVUQxtDyMIVccZ9o1hlxxnGnXGBTHUi4FcYyCOEZBHFPmeeJeFiwLM/uqmZ0ws31d19aZ2fNmdrB4P971uYfM7JCZTZrZPaMaXHrLFseZZItjGZzW7f4pimMUxTGK4pgqg3hGP1tujwP3zrr2ILDb3bcAu4uPMbNbgO3AB4vHPGJmS4c27YAy7RpDrjjOtGsMueJYu8aL8jgNX7dHTUcn4hTFMYrimDpEMfQRxu6+B/jxrMvbgJ3F7Z3AfV3Xn3L3N919CjgE3DakWRclWxxnki2OM1EcD6Yt6/aoKIhjdJ44RueJY6o+OjHboEVxvbsfAyjeX1dc3wC81nW/6eLau5jZZ81sr5ntvXDp3IBjxGSK40y7xpArjjPtGoPieIiGum6ff+P8SIcdFUVxjII4RkEcU6cgnjHsmrAe17zXHd39UXff6u5bly9ZOeQx5qY4bi/FschABlq3x9Y2L5gUxTGK4hhFcUwdoxgGD+PjZrYeoHh/org+DWzqut9G4PXBx5PFyhbHmWSKY+0aD0XadVvniWN0dCJOURxT1yiGwcP4OWBHcXsH8GzX9e1mtsLMNgNbgBcXN+LwZdo1hlxxnGnXGBTHEtLodXtQCuIYBXGMzhPH1O08cS/9/Li2J4HvAhNmNm1m9wNfBD5hZgeBTxQf4+6vAE8D+4F/AR5w91o+Y7LFcSbZ4jgTxXF/2rpuRymKYxTFMQrimLoH8YxlC93B3X9tjk/dNcf9HwYeXsxQZVm57yjnPtTze0xaZ3zyAqcmllc9RmlWH17CmZsuVT1GKc5OrWHV5tNVj1GaXUcmuPuGyarHqLU2r9v9UBDHKYpjFMUxTYliSPA338llmY5UZJPpSIXIfBTFMTpPHKOjE3FNimJQGKc7UpEpjrMdqcgUxzpSIb0oimMUxDEK4pgmnCfuJVc5zEFx3F7Z4jgTxbF0UxTHKIpjFMUxTQziGaqGQrY4ziRTHGfaNQbFsehHsUXp6EScojimyVEMCuMrZIrjTLvGoDgWaSMFcYyCOEbniWOaenRitjy1IO+SLY4zyRTH2jXOSVEcoyiOURDHtCGIZyiMZ8m0awy54jjTrjEojqWddHQiTlEcoyiOaVMUg8K4p2xxnEm2OM5Ecdx+CuIYnSeO0dGJuLZFMSiM55QpjjPtGkOuOM60awyK4zZTFMcoiGMUxDFtOU/cS55CGIDiWNogWxxL+yiKYxTFMYrimLYG8QyFsbwjUxxn2jWGXHGsXeP20HniGB2diFMUx7Q9ikFhvKBMu8agOJZ2UBw3n4I4RkEco/PEMW0+OjGbyqAP2eI4k0xxnGnXGBTHTaYojlEUxyiIY7IE8Yw8VbBImeI4064xKI7bTHHcPIriGEVxjKI4JlsUAyyregCpp/HJC5yaWF71GDICZ6fWsGrz6arHELmCgjhGQRyjII7LGMWgHeOQTLvGkGvnONOuMeTaOdaucf0pimMUxTGK4phM54l7yVUDQ5AtjjPJFseZKI7rS1EcoyiOURTHZA7iGSqBAWSK40y7xpArjjPtGoPiuG70o9jiFMUxiuIYRXGHuXvVM2BmJ4EfVj3HgH4G+FHVQwxIs1dDs5dv1HO/393fO8KvXzsNXreb+hwGzV4VzV6NUc4+55pdizBuMjPb6+5bq55jEJq9Gpq9fE2dW4avyc8FzV4NzV6NqmbP87qxiIiIiMg8FMYiIiIiIiiMh+HRqgdYBM1eDc1evqbOLcPX5OeCZq+GZq9GJbPrjLGIiIiICNoxFhEREREBFMYiIiIiIoDCeF5mtsnMvmNmB8zsFTP7XHF9nZk9b2YHi/fjXY95yMwOmdmkmd1T3fTvzLPUzP7TzL5VfNyI2c1srZl93cxeLf7939Gg2X+3eL7sM7MnzWysrrOb2VfN7ISZ7eu6Fp7VzD5qZi8Xn/tLM7OKZv+z4jnzAzP7ppmtrePsMhpas6ujNbu0WbVmj3p2d9fbHG/AeuAjxe3VwH8DtwB/CjxYXH8Q+FJx+xbgv4AVwGbgMLC04l/D7wH/AHyr+LgRswM7gd8ubi8H1jZhdmADMAWsLD5+GvjNus4O/DzwEWBf17XwrMCLwB2AAf8M/FJFs98NLCtuf6mus+ttZM8JrdnVza01u5x5tWaPeHbtGM/D3Y+5+/eK22eAA3R+E22jswhQvL+vuL0NeMrd33T3KeAQcFu5U19mZhuBXwG+0nW59rOb2TV0fgM9BuDuF9z9DRowe2EZsNLMlgHvAV6nprO7+x7gx7Muh2Y1s/XANe7+Xe+sWn/X9ZhSZ3f3Xe7+VvHhfwAb6zi7jIbWbK3ZA9KarTX7HQrjPpnZjcCtwAvA9e5+DDoLMXBdcbcNwGtdD5surlXlL4A/BC51XWvC7B8ATgJ/W7yk+BUzu5oGzO7uR4E/B44Ax4DT7r6LBszeJTrrhuL27OtV+y06uwnQvNllkbRml0prdrW0Zg+RwrgPZrYK+AbweXf/yXx37XGtkp+HZ2afBE64+0v9PqTHtap+lt8yOi+3/LW73wr8lM7LQ3OpzezF2a5tdF76eR9wtZl9er6H9LhW15+hONestfs1mNkXgLeAJ2Yu9bhbLWeXxdOaXTqt2fXUmHWvTmu2wngBZnYVnQX2CXd/prh8vNjOp3h/org+DWzqevhGOi/JVOHjwKfM7H+Ap4A7zezvacbs08C0u79QfPx1OotuE2b/RWDK3U+6+0XgGeBjNGP2GdFZp7n88lf39UqY2Q7gk8CvFy+1QUNml8XTml0JrdnV0po9RArjeRTf6fgYcMDdv9z1qeeAHcXtHcCzXde3m9kKM9sMbKFzSLx07v6Qu2909xuB7cC33f3TNGP2/wVeM7OJ4tJdwH4aMDudl+NuN7P3FM+fu+icc2zC7DNCsxYv3Z0xs9uLX/NvdD2mVGZ2L/BHwKfc/f+6PlX72WXxtGZrzR6A1myt2VeKfrdepjfg5+hs0f8A+H7x9svAtcBu4GDxfl3XY75A57snJ6nJd7cDv8Dl73BuxOzAh4G9xb/7fwTGGzT7nwCvAvuAr9H5rtpazg48Sedc3UU6fxK/f5BZga3Fr/cw8FcUf6tmBbMfonMubeb369/UcXa9jew5oTW7upm1Zpczq9bsEc+uvxJaRERERAQdpRARERERARTGIiIiIiKAwlhEREREBFAYi4iIiIgACmMREREREUBhLCIiIiICKIxFRERERAD4f9/9BBtbJHgdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the ROI with bottom left corner and rotation of axes\n",
    "# original with origin in upper left\n",
    "e0 = 377270.; n0 = 4624440.\n",
    "rot = 52.\n",
    "\n",
    "# _v2 with oring in bottom left\n",
    "e0=376246.28785643\n",
    "n0=4625240.44609256\n",
    "rot=-38.\n",
    "\n",
    "print('Shape of xrot, yrot: ',np.shape(xrot),np.shape(xrot))\n",
    "shp = np.shape(xrot)\n",
    "xu, yu = box2UTMh(xrot.flatten(), yrot.flatten(), e0, n0, rot)\n",
    "xu=np.reshape(xu,shp)\n",
    "yu=np.reshape(yu,shp)\n",
    "print('Shape of xu and yu: ',np.shape(xu),np.shape(yu))\n",
    "\n",
    "# Corners in rotated coordinate system\n",
    "LLr = np.asarray( [xrot[0][0],yrot[0][0]] )\n",
    "LRr = np.asarray( [xrot[0][-1],yrot[0][-1]])\n",
    "URr = np.asarray([xrot[-1][-1],yrot[-1][-1]])\n",
    "ULr = np.asarray([xrot[-1][0],yrot[-1][0]])\n",
    "\n",
    "# Corners in UTM coordinates\n",
    "LLu = np.asarray( [xu[0][0],yu[0][0]] )\n",
    "LRu = np.asarray( [xu[0][-1],yu[0][-1]])\n",
    "URu = np.asarray([xu[-1][-1],yu[-1][-1]])\n",
    "ULu = np.asarray([xu[-1][0],yu[-1][0]])\n",
    "\n",
    "print('Corners in rotated coordinates')\n",
    "print('Upper left',ULr,' Upper right',URr)\n",
    "print('Lower left',LLr,' Lower right',LRr)\n",
    "print('Corners in UTM coordinates')\n",
    "print('Upper left',ULu,' Upper right',URu)\n",
    "print('Lower left',LLu,' Lower right',LRu)\n",
    "\n",
    "corners = np.asarray(  [LLu, LRu, URu, ULu, LLu ])\n",
    "# [[xu[0][0],yu[0][0]],\\\n",
    "#                        [xu[0][-1],yu[0][-1]],\\\n",
    "#                        [xu[-1][-1],yu[-1][-1]],\\\n",
    "#                        [xu[-1][0],yu[-1][0]],\\\n",
    "#                        [xu[0][0],yu[0][0]]])\n",
    "\n",
    "# write the UTM coords of the corners to an ASCII file\n",
    "print(\"Corners: \")\n",
    "print(corners)\n",
    "fn = 'sandwich_corners.csv'\n",
    "np.savetxt(fn, corners, delimiter=\",\")\n",
    "\n",
    "fig = plt.figure(figsize=(12,6))\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "plt.contourf(xrot,yrot,xu,5)\n",
    "ax.set_aspect(1.)\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "plt.contourf(xrot,yrot,yu,5)\n",
    "ax.set_aspect(1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ksize=3\n",
    "do_filter = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop through the maps and interpolate onto the rotated coordinate system\n",
    "Some of the DEMS have holes w/ NaNs.\n",
    "\n",
    "Here is one option for smoothing the raw data before interpolating:\n",
    "https://stackoverflow.com/questions/18697532/gaussian-filtering-a-image-with-nan-in-python\n",
    "\n",
    "Strategy: Use a 2D convolution kernal that is sharp, but covers a 5x5 footprint. Using the"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2016-01-22_SandwichTNB_PT_DEM_10cm_trimmed.tif\n",
      "(10226,) (11987,) (1, 10226, 11987)\n",
      "(11987,) (10226,) (10226, 11987)\n",
      "1 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2016-01-25_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10156,) (11420,) (1, 10156, 11420)\n",
      "(11420,) (10156,) (10156, 11420)\n",
      "2 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2016-02-11_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10204,) (11673,) (1, 10204, 11673)\n",
      "(11673,) (10204,) (10204, 11673)\n",
      "3 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2016-03-30_SandwichTNB_AS_DEM_10cm_trimmed.tif\n",
      "(9858,) (11504,) (1, 9858, 11504)\n",
      "(11504,) (9858,) (9858, 11504)\n",
      "4 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2016-09-21_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(12379,) (12221,) (1, 12379, 12221)\n",
      "(12221,) (12379,) (12379, 12221)\n",
      "5 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-01-09_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10491,) (11301,) (1, 10491, 11301)\n",
      "(11301,) (10491,) (10491, 11301)\n",
      "6 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-01-25_SandwichTNB_DEM_10cm_trimmed_channel_mask.tif\n",
      "(10175,) (11191,) (1, 10175, 11191)\n",
      "(11191,) (10175,) (10175, 11191)\n",
      "7 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-02-14_SandwichTNB_DEM_10cm_trimmed_channel_mask.tif\n",
      "(10544,) (10865,) (1, 10544, 10865)\n",
      "(10865,) (10544,) (10544, 10865)\n",
      "8 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-03-16_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10578,) (11474,) (1, 10578, 11474)\n",
      "(11474,) (10578,) (10578, 11474)\n",
      "9 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-04-28_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10975,) (11993,) (1, 10975, 11993)\n",
      "(11993,) (10975,) (10975, 11993)\n",
      "10 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-05-04_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(11815,) (12142,) (1, 11815, 12142)\n",
      "(12142,) (11815,) (11815, 12142)\n",
      "11 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2017-09-18_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10236,) (11762,) (1, 10236, 11762)\n",
      "(11762,) (10236,) (10236, 11762)\n",
      "12 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2018-01-10_SandwichTNB_DEM_10cm_trimmed.tif\n",
      "(10000,) (11425,) (1, 10000, 11425)\n",
      "(11425,) (10000,) (10000, 11425)\n",
      "13 D:/crs/proj/2015_Sandwich/2017_Karen_Sandwich_maps/2018-03-09_SandwichTNB_PT_v2_DEM_10cm_trimmed.tif\n",
      "(11721,) (13544,) (1, 11721, 13544)\n",
      "(13544,) (11721,) (11721, 13544)\n",
      "Wall time: 11min 15s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "dslist=[]\n",
    "for i in range(nf):\n",
    "    fn = fdir+fnames[i]\n",
    "    print(i, fn)\n",
    "    \n",
    "    # open the tif with XArray as a DataArray\n",
    "    da = xr.open_rasterio(fn)\n",
    "    \n",
    "    print( np.shape(np.flipud(da['y'].values)), np.shape(da['x'].values), np.shape( np.flipud(da.values)) )\n",
    "    x = da['x'].values\n",
    "    y = np.flipud(da['y'].values)\n",
    "    \n",
    "    # Not sure how da.values got a singleton dimension, but squeeze gets rid of it.\n",
    "    # However, make sure to squeeze before flipping\n",
    "    z = np.flipud(np.squeeze(da.values))\n",
    "    print(np.shape(x),np.shape(y),np.shape(z))\n",
    "    \n",
    "    if do_filter:\n",
    "        # smooth with 2D running median\n",
    "        zs = signal.medfilt2d(z, kernel_size=ksize)\n",
    "    else:\n",
    "        zs = z\n",
    "    \n",
    "    f = interpolate.RegularGridInterpolator( (y, x), zs, method='linear')   \n",
    "    \n",
    "    # Array for interpolated elevations\n",
    "    zi=np.NaN*np.ones((ny,nx))\n",
    "\n",
    "    # this is a slow iteration through all of the points, but allows us to skip ones that are outside\n",
    "    for ij in np.ndindex(zi.shape):\n",
    "        try:\n",
    "            zi[ij]=f((yu[ij],xu[ij]))\n",
    "        except:\n",
    "            zi[ij]=np.NaN\n",
    "            \n",
    "    da = xr.DataArray(zi,dims=['Cross-shore','Alongshore'],coords={'Cross-shore':ycoords, 'Alongshore': xcoords})\n",
    "    # dslist.append(da.chunk()) # this does not work...can't concatenate a method\n",
    "    da = da.chunk()\n",
    "    dslist.append(da)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsa = xr.concat(dslist, dim='map')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (map: 14, Cross-shore: 280, Alongshore: 1300)>\n",
       "dask.array<shape=(14, 280, 1300), dtype=float64, chunksize=(1, 280, 1300)>\n",
       "Coordinates:\n",
       "  * Cross-shore  (Cross-shore) float64 0.5 1.5 2.5 3.5 ... 277.5 278.5 279.5\n",
       "  * Alongshore   (Alongshore) float64 0.5 1.5 2.5 ... 1.298e+03 1.3e+03\n",
       "Dimensions without coordinates: map"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsa.to_netcdf('Sandwich_1m_DEM_v2_nf.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAELCAYAAAAcM9T7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29e/RdVXnv/X0gCCr3IPRFJBxRCEILhFJ6ahEcSbUeX6sWy6jGQkotUU/air42nJZwKeoBz3FQT+MlOUXDJViVothTFRpqULSNEBjgwUQHqQYQkJtyF1Se94+5Zpi/mbnWmnNd596/72eMPZLfXre59/ruZ37nMy9LVBWEEEIIIYTkyg5jF4AQQgghhJAqaFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYSQxojIQSKiIjJn7LIQQqaXTg2riOwsIheLyFYReUxEbhGR1xXbflNE/kVEHhaRB0Tk8yLy/1Scq7cgKCKHiMjVRTkeFpFrRORQb58zROQ+EXlERD4lIjs725aJyE0i8rSIrAmc/wUi8nERebA4/utdfwbSLTXatVp83HmtqDgXtUsGo+yeisjzRORKEflhoccTa85zoog86+n8n3ooLzVMAAAicrmI3Csij4rI90XkHc62hSKyWUSeFJGvici8ivOcWGj8Ku/9I4v313dYZup3JLrOsM4BcBeAEwDsAWAFgM+JyEEA9gKwGsBBAOYBeAzApzu+fix7AvgSgEMB7Afg2wCuthtF5LUAzgSwEKa8LwVwnnP8PQA+AOBTJedfDWBvAIcV/57RaelJH1Rp17Knqu5avM4fvoimDKB2yUyq7ukNAN4O4L7Yczka31VV39BVIR2oYWL57wAOUtXdAfwegA+IyDEisg+Aq2Di8N4AbgLw2ZpzPQDgt0RkrvPeqQC+37RwJUkH6ncsVLXXF4DbAJwUeH8BgMcqjjsIgAKYU/z9egC3AHgUxlicG9j3VAB3AngQwF8nlHHv4vi5xd9XAPiQs30hgPsCx30AwBrvvUOLMu7e93fLV78vq11fixHHUbt8Df4K3VNn290ATqw5/kQAd5ds2wGmEt4C4CEAnwOwd7HNavh0mMr4XgDvSyg3NcyXvXf3Aji50NK3nG0vBPAUgPklx55YaPyTAP5r8d6OxXtnA1jv7PvRIg4/CmAjgOOdbecCuBLA5cX2d0SUm/od6NXrGFYR2Q/AIQBuD2x+Vcn7ZTwB4BSY1s3rAbxLRN7k7fPbMDd8IYCzReSwyHO/CkZQDxV/Hw7gVmf7rQD281puZRwHYCuA84qU/ndE5KTIcpBMKNHuVhG5W0Q+XWQAYqF2yaTz5wDeBNMDsT+AnwD4mLfPqwG8HMBrAJwpIosiz00Nz2KK7u8nAWyGMaxfhqcBVX0CprF0eM3pLoWJtQDwWpj4fY+3z40AjoIxmlcA+LyI7OJsfyOMad0TwNqIj0D9DkRvhlVEdoK52Zeo6mZv26/BtHreH3s+VV2vqt9R1WdV9TYAn4EJni7nqepTqnorjEiOjCjnATCB973O27sCeMT52/5/t4iiHgDgiOKY/QEsA3BJggEhIxPQ7oMAjoUZynIMjA5iAhkAapdMFPuLyE+d18nF+0thMv93q+rTMJmot3hdpuep6hOq+h2Y4V5vrbsYNUxU9d0w9/Z4mGEAT2N7DaD4u1IDqvotAHsXY0pPgTGw/j6Xq+pDqvoLVf0IgJ1hkgWWf1PVLxbx+qmq61G/w9KLYRWRHQBcBuAZmC/c3fYyAF8B8Beq+o3ivQPdgf4l5zyuGHj9gIg8AuCdAPwslztO60kY4UBmTiI40DnniwBcC+DjqvoZ59jHAezu/G3//1jEx38KwM8BfEBVn1HV6wF8DSbrQDInpF1VfVxVbyoC3I+L918jIrtTu2QSqdDtPaq6p/P6XPH+PABfsEYWwCYAv4QZw2e5y/n/VpgKlxomtajqL1X1BhjD9i5srwEUfz8WEXMvg4nRrwbwBX+jiLxPRDYVE5p+CjNnwY3Hd3n7U7+Z0LlhFREBcDFMIDtJVX/ubJsHYB2A81X1Mvu+qt6pzkD/klNfATPQ+SWqugfMWBWJKZPOnERwZ1GWvWDE9iVV/aB3yO2YmeE6EsCPnZR/FbfFlInkR5V2PdQeQu2SSSRSty53AXidZ2Z3UdUfOfu8xPn/gSi6YqlhksAcAAfD04CIvNC+H6HdywC8G8CXVfVJd4OIHA9gOcw42b1UdU+YbKYbj9U9hvrNhz4yrJ+Amd32BjedLiIvBvCvAD6mqp9scN7dADysqj8Tkd8A8LamBRSR3QFcA+CbqnpmYJdLAfyJiLyiEOZZANY4x88pxrzsCGBHEdnF6Rr7Oszkmf9W7PdKmAHh1zQtLxmMMu0eJyKHisgOxTik/wUziN/vsiqD2iW9UXVPxSzXZsfnPa/YFtVY8vgkgA8WSQeIyItE5I3ePivELMlzOIA/RsmsbmqYAICI7Csifygiu4rIjmJm178Vxid8AcARInJSca/PBnCbP7wwhKr+AGbI1V8HNu8G4BcwKwrMEZGzsX0mt67c1O9YxMzMUlXs/Zd3qP/y94HpNlIAP4NJi9vXYgDnFNvc9x8vux6em3m6Y/H3W2C6mR4D8H8ArARwubfvHOf49SiZ4QczI1thJsO45TnQ2ee9AH4MM2Pv0wB2dradWxzvvs51th8O4N+K838XwJtjv+cuX2X3KeUeu6+jv/ueyu12H3+/MT576qtGu28F8IPift4LE5B+hdqdjJfVZd1+MfsMXe66mFt3TwH8MLDtoJLznIjqVQLeC+B7hY63oJgJje1XCbgPwF+Wfa5p13DbmOfe72ef/eR2993/O7eYG6PZ4j69CMD1AH5a3OfvAPhTZ/simIlYTxUxMajbCO2+A8UqATAm8eLievcC+MviN7LI0dblVZ9vmvXblXZSYm5K3JXiC6pkwaYzdOuaZdu9/9CFBzdpqUchZmLW19Wk7EkD5i7fsu3mVt0rd78q5i1ZiZAO/H0AzNivT53UsWDTGXrzYRcNen1qNy8WbDpDAaBKB3OXb9F5S1ZW7jMktsyWXH5PpB4bT5vcJz8Wh+KpT44x1/7flomanQzmLt+ibe+V1bD1C1Xns1qJ1W70k3hizEpXFBNfToZZLJh0QKwp7YKhdJIj1G5eGCNav5+t9HPALfNs/i3NJoaMz31Spt0ujBDJmzINu/febcxZraTEuFrD2uSkHXAnTLfSHw950Wmhj+CXU4Uei5+lGghqNwNCv4FpqDQnvfyzCb+Ha1pMaRWTWE+QmaTGyRhd+/s0/S1UGlbbTTY0qnrA4BedcGZDMJwEqN1xqfsdhLprY7OwORA7zIcMT5X2ZkN8ruoZoFYni7Z67cs3tlolYO7yLWpfXRVoNtPkuxzi+3fFN0kt6I3z549dBNICN740acUD9XrNMXbF/sZyLPtspet7MW/JyqS4m/vQEWo1b/ryEV37hcpJV26GNeYH0Ucrahq68WIIicXPAg1bomqxpQTIMe7fgk1n6Mb583HM5s3ZTKQhcfTZuvcH9/vXyiHWqK7SYzZvbmRCupow0cW5JoWmk6T6MKl1xGpi7JhbVs7Zoqkh6OK3OpSGu/ILlYY1NNuvzcWqiPnipkns9vM+eME6AMA+Z5rHbqc0EPpk0oOn6ioFQMM6IG1XZBiqG6pKtznEGGtYgWZxoIvZ6W3PN0n4M5VDn7erir2LjNNYyaM63JgL5P87mwZiVkDxGaOhBXTjFyrHsNouVSvAGFJbqylfXpvlQsak6jNao9oXdnWHPrry3ZUjcjHaIcZY2mq2ktIj0qQLvwty1qpPk9VZYu9BbOyd1LgbQ2j8ctMKfajhUn7czVXHVWXLsWdjkin7zXfd0PLr+9RztNVqb8taxaSrm36ZOQRQvwyh9fOGDiRlImoaSG33TpNr5kBKQ4t0g9HDRaXbxzKpIdzr5lrpN6WPCiyHuNuUqs/dpEE/dtzzx7dOun6bDP2r0/Ik6jQG93O7ja3U33ZT09kG//hU3dYa1r4q/Zgvty4T4i5Qu3H+fIgs7UWgsbM/3cAxRADpI2j6E5WaZNlzZOgs62zO6lb9Xsau6MvIudJvWrYm5jQm+1w33j62bP4xMcMS+hzLX/c956rdnEnRbkyDKCa22OtN03jsLnQ/hF9o6xPq7lN0hhXoriut7Mtv2n29dc0y4IJ10YHUX66rq4rKP0/sjzX0uYcMjjGz6afBuM5mE9kXodY+K/buydlQA80q1CbDwfpgUldBmRTa9M52xdjmte9xo0PrdqwVeGoNa6hbuM3KAWWPniv7u+w9//plY0H77IKs+h7stthrDSW4NkKbZLNqaTJIvcn5cyf2e6j7/UzK+qXTQNdjb6tiTs4Guc33QEM6Du4961JbTe5nV0MJmnbBp372sTTregW/7h/LC1SuErDqrNP09POPqyxcF8Lr44Y0/WH0OUlpDLpsCcWINPR95zBjtYy2xjXWoI6V2a3qzQDCvQK50VTDqUG17rnXQxGj3aErvVzMa0zP2KRpuIvKP4eF+leddZoCQJ1n6JtctOoSM3YzJ92mxNw29zpVt5WG9eb11+mN69Zu+3vV4j2iL5xCTjdqGujapFa1tHxyMaw2eAL1ATTVTPaRRbVl8JeSa/rdTUqmN8TQjSxLLoY1Vrs5jbfMNXExBm3128QA5GBYU/zCmPRpaCdVw13F3CZJAkuMZmszrMcuWgwAsELsMuM6qTc3N4YaTzJJGdab11+3TdhWuzEBNGReczJ/Nx92keRUnq7oU8OxQTQXw+prd+yMVVP8nqo2S+LkTtf6bVPxA+MlCXy/kKtpTWFaddtXzE1NElhiNJs06QoAVq/YgKUwxjVUuK5uaujLnMSg3SV9CGz1ig3b/u8GG3t/LU3Nak4sXftIbeWfuxnMvXx1dKnhWO3GLM+Ws3ZXr9iAjRNoWseeGNI1XWk3pFtgpnZjNGvJWbs25gKTW39Psm679gxl2j36hDtmXLOve12bYfUL56b8LW4F0aag/pfrfjmha7l09QX1mRpPNeF9Z05D3++xixbPyKa7NK30xx4SENJuV5r1qbtnkxq0U2mr3dUrNpT+1u12ny60m0OG1Wo31LtlmS06GoM+tVumW6D5fc6lV2usmAuk37Mhfj9DrqzTZWMqVbtd6bZ1hrXKpLq4H2RjRSUTS1ULNPSldXXdpoRusi/WWLNaJTz72asq8raUVfh15NbKd39IIWI169+jusBgz+vqFnjue11a/F31/TYNcGN2q8foFojX1ennH1drWrfff63z//hjc6NMu+7vf6xlZVz6qoi7+GyxZYv9PfeppaNPuGObXsvqt0nAj3k+XfuENsRqzJ/DUbZP1fnH+q0O7RfcmFuXZW06hKV20lXMSeoMrSXWAIeOaXP9uhvWVZAIXafOsMZW9LHX65Mmkz+A8cewltFEh1XHuBVOjG7Lztf0vrbRcZdaiinHkNpN6RnIIcOaql37XZZ973WNN//42Jg71D2s0lPXZZgk7VbF3VxjLpAed1O0W8UYDdiycvVRlkn2C70aVisgP/tZh7+/awJCQTJ0nVRSzXTba7QVRYrp6FKAVZnENpPsxjIAVdoNZT9tYCwLkCHtuvuVnbMN1G4crj7brGiRg1kFwtoti4VVmrX4ja2yv93ztaHOzKY26qquY5lE7fYVc4G8DGsffqFsny5J6Rqvour3OYm6Bep7HlOTW50YVnfWn0+Myewi+PnXqjIJY+GLMUUYOWTFXHHFrEGb0v0/lgko025VIyimEZVCjEEYm2nRbsxM3lwbVyHKTGuVpi0pBrSsAeeeK8eYC4yj3b5iLtBOuz5jrxLg04dfqNJuk97ctjTtDZ7UmAt05xd6N6wxxApwwYkLxQbomADqvxdi6MBaJsSuBJYyAS2FMvHVLcodG0BzM6wuVUE0VT8xWYS6TJh//FD0FdSruuv60G5bcjKslpSegqr3y6iLrXU9Dy6zQbs56fahCw+WsifRjf0I6tjhAS6pfqGq18ElRrshYrO3TXXfh3Zzj7llmu3EsALbz1htQlUGwN+nanZh6sD0kFDd99qakq7G1vjHDIWftm8aPKtaWGMFztCMVZ+qrqWyexyT4bL4urX/r6OqN6Fp5jc2uNty56xboDpwuk+4q3vanavdsSv5MvzGV1VWKka77j52wo9PH9rtmmnSbllFHkMo/uakZavfGD3av0P7+fi/A0uXfqGuDHXE9ooA+Wu3qV+wZtTXt9VtilZrDWuIqu7WlArdJVSRDpUq74omY16qJmrV0dXj/vp6JvfG+fMhsjSbwFk1NtDim4CYCX23XP+yGfumBpJp0W6K5sq023bB9BjKtJtTJR+izLgC4co/RbvuOoopzDbtdrlsXR8JAkuuWnYTCUBaAiB03GyLuaHjxvQLKQ2tLuJuI8PqE7teaxVtBQgMI8KyVlDKmpGpJrVujFPsgOeqc8cSWwZLTobVx3+UYFNc7aY2sppk4rumL+3GjM1LmSQVOncKdRV9rpV8iK6129SwApOh3Sa6tYR00+bBNrHZf586/bplyjnuWnzzmkKTmBta6m4s7ebkF+rOFzp3F0bVMrhh9Vl11mlal94Opd/dY0I3KXZtxVhDWWUYYpbC6npMCNDvo23bjvPzr1eXbZiEoOmjukpj1q8LVYZNNZGagY/twq263tjjmVK0W2YgYq85ib0BsZQ9ZMDS95g2n5QYHRqjX5V16mt2fhvt1jW4uo65IUJxeNK07MbdWNqst+zqrck60anXzdEvVGm3bsm0GLMa28BK0WovhtXHbUm5Iqur9KvWMPXxz5VacYcM6+oVG7YZ61DLrI0Au54wYgmJpM35Y7NjISYtaJahumq7H0nd/U/Rbgp15naIBfP9mfldEavd2MxUSkYKmHy9uua1rKeqzhiWLRHmHt/VOsFdx9Q6uo65TWb1d1XZWyaxZ6sO3y9Y2mq3K7pIqqUyrdpN1ekghjWEbVGVDQWIeWTpJDzNxn84QKzw/Jve95OkUrO1MQFhkoNmHaqrtGz8X5PH7bq/hTbds13iajclYE6qdqdNr7bi95MD7rhBV2ttnmaWM02127VuQxV+0+fUz7Zkgb92bxO/0JahkgFN/AIwGdptmyQYzbBa3OWsXDGkBk83a+u3ypp2GcRk0JpQt2RU0+ObkCK61B/+tATNMqx2Uyp9oP/Wf53em2o3pYGVenwKNsPadSNrmvXqZ61uuf5lyTHX0pd+pznm+hV+F5V9GdOqYzdJ0ES7Y8bdSdZuk2v6dJUkmFjDClRnq4bqWkqhrfBizpVCn2YVmN7A6XLz+us01bAC42SrctBt3Tlj6UO7067XUMa1qWENPVa6z7ibi3aH1C3AuFuGnSdjGduwVjFbtQt0nyQY3bASQgghhBBSxQ5jF4AQQgghhJAqaFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYQQQgjJGhpWQgghhBCSNTSshBBCCCEka2hYCSGEEEJI1tCwEkIIIYSQrKFhJYQQQgghWUPDSgghhBBCsoaGlRBCCCGEZA0NKyGEEEIIyRoaVkIIIYQQkjU0rIQQQgghJGtoWAkhhBBCSNbQsBJCCCGEkKyhYSWEEEIIIVlDw0oIIYQQQrKGhpUQQgghhGQNDSshhBBCCMkaGlZCCCGEEJI1NKyEEEIIISRraFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYQQQmYdInKQiKiIzBm7LKSeTg2riOwsIheLyFYReUxEbhGR1xXbniciV4rIDwuBnBhxvr8SkR+IyOMicreIfLbL8jrXeb2I3CAiPxWR+0Tkf4vIbt7n+pSIPFpsf693/GoR+Z6IPCsiSwLnf6mI/J/iO3lQRD7cx+cg7RCRvUXkCyLyRKHhtznb3iEidxRa/KqI7F9xnt6CoIgcIiJXi8gDIvKwiFwjIod6+5xR6PSRQrc7O9uWichNIvK0iKzxjltcfD77erL4HMd0/TlId9TEXatF976uqDgXtUsGo+qeOvucU9zLRRXnObGof10N/FMP5aWGR6TrDOscAHcBOAHAHgBWAPiciBxUbL8BwNsB3Fd3IhE5FcAfAVikqrsC+HUA13VcXsseAD4AYH8AhwE4AMD/cLafC+DlAOYBeDWAvxSR33W23wrg3QBu9k8sIs8D8C8A/hXArxTnvrzzT0C64GMAngGwH4DFAD4hIoeLyAkAPgTgjQD2BvADAJ8ZqYx7AvgSgEOLcn4bwNV2o4i8FsCZABYCOAjASwGc5xx/D4zWP+WfWFXXququ9gWj6f9AQNckK+riLgDs6dzb84cvoikDqF0yk9J7CgAicjCAtwC4N+ZcrgZU9Q0dltNCDY+Jqvb6AnAbgJO89+4GcGLNcSsB/G3F9j0AXAwj5B/BiGDHYtsSAN8E8HcAHgGwGcDChDL/PoDvOH//CMBrnL/PB/APgeNuALDEe+90AN/o+3vmq90LwAthzOohznuXAbgAwP8E8DHn/f0BKICDS851ULF9TvH36wHcAuBRGGNxbmDfUwHcCeBBAH+dUO69i+PnFn9fAeBDzvaFAO4LHPcBAGtqzv01AOeMfW/4Sn/ZuOtrMeI4apevwV9l9xTAVwD8FwA/hElelR1/IoC7S7btAGMitwB4CMDnAOxdbLMaPh3GTN4L4H0J5aaGB3z1OoZVRPYDcAiA2xsc/u8AThGR94vIr4vIjt72SwD8AsDLABwN4DUA3uFsPw6mdbIPgHMAXCUie0de+1W2zCKyF4xBudXZfiuAwyPP9ZsAfigiXxEzHGC9iPxq5LFkOA4B8EtV/b7znr3PUrws9v9HRJ77CQCnwLTOXw/gXSLyJm+f34ZptS8EcLaIHBZ57lfBBMSHir8Px/Za3U9E5kaeDwAgIvOKc1+achwZn5K4u1XMsKpPi8g+CaejdskoiMgfAHhGVb/c8lR/DuBNMD0Q+wP4CUxvmsurYXpRXwPgzKrhBx7U8ID0ZlhFZCcAawFcoqqbU49X1csB/BmA1wK4HsD9InJmce79ALwOwHtU9QlVvR/ARQD+0DnF/TAZ2p+r6mcBfA8m4NaV+3dgMgZnF2/tWvz7iLPbIwB2QxwHFOX6XzA/ln8GcHUxVIDkw66YeY+B5+7zlwGcLCK/JiLPh9GGAnhBzIlVdb2qfkdVn1XV22CGE5zg7Xaeqj6lqrfCBLkj684rIgfABF53TLX/Oez/Y/VqOQWmZ+AHiceREQnE3QcBHAsznOkYGB2sjT0ftUvGQER2hRmG9Z6Ew/YXMw/Fvk4u3l8Kk/m/W1Wfhhni9xZvnPZ5hZf4DoBPA3hrRBmp4YHpxbCKyA4w3anPAFgWsf+B7mBj+76aMR2LYFr37wTwN8UYkXkAdgJwrxUngFUA9nVO+yMt8uoFW2EEfbxzrRmZXxH5TZiU/lucTJstz+7OrrsDeKzucxU8BeAGVf2Kqj4D0708F2asLMmHxzHzHqP4+zFVvQ4mS/+PMDr6Icz9v7tMuy4icpyIfK0YqP8IjJb9LJc7rvtJFA0lbxD+gc45XwTgWgAfV1V3PK3/Oez/Y/VqOQWmF4NMCKG4q6qPq+pNqvoLVf1x8f5rRGR3apdkzHkALgsZtgrd3qOqezqvzxXvzwPwBccrbALwS5gxqJa7nP9vhUkuUcOZ0blhFRGBGVu6H8zY1Z/XHaOqd+rMwcb+9p+r6udhxmUdASOupwHs44hzd1V1u+lfXJTFciCMoL/hXGvb/iJyNMxg6tMKg2Kv/ROYcS1u1uBIxA9zuA0mG0fy5vsA5ojIy533tt1nVf2Yqr5cVfeFMa5zAPzfOu0WXAGjrZeo6h4APomZQwxK0ZmTCO4Etg1TuRbAl1T1g94ht2N7rf7Y6bKqRUReCROwr4w9hoxLQty1sUioXZIxCwH8uZjZ9vcBeAnMRMLlkbp1uQvA6zwzu4uq/sjZ5yXO/w+EGc9KDWdGHxnWT8BkD9+gqk+5G8Qsv7JL8efzRGQXz1S6+y4Rs9zUbiKyg5hlWg4HsEFV74URzEeKTMEOInKwmNncln1hBL9TMRbmMJiu3dC1jgDwVQB/pqqhpTAuBXCWiOwlIvMB/CmANc7xzys+lwDYqfhc9ru9HMBvisgiMeNw3wPTTbcp+O2RUVDVJwBcBZPFf2ERNN4I4LLifh4hhgMBrAbw0aIxE8NuAB5W1Z+JyG8AeFvdAWWIyO4ArgHwTVU9M7DLpQD+REReUQTWszBTq3MKre4IYMfis/lLGJ0K4B9VNTUrQMYjGHeLDOmhRYycCzM0ab2q+sNfyqB2SW9U3NOFMMmpo4rXPTBd+/7Y0xg+CeCDYsaFQkReJCJv9PZZISIvEJHDAfwxgOASmtTwuESttTd3+ZbSDOFDFx68zXAWglgKk/28z/GiS1V1Lcw40nnFe9cU//4nmC5Wn0cB/BWM4dsRJk3/LlW9odh+CswM7u/CBNX/AHChc/wGmEHUDwL4MUw3f1kr530AXgTgYhG5uHhvq5OBPQemQtgK08V/oap+1Tn+Wjw3ruu3YAzNq2Eqhu+JyNthfjT7wixR8XvF8IBesffNvUdl+8Qwb8lKbF1TPcJj3pKVADBjv6rr983c5Vt03pKVuPmwi2LK8G6Y5Ubuh5lN+i5VvV1E9oTJNB0M07XzaZilg+qw3+27YRpXK2HGY38OZphLE94MMybxcJm55u8riszDV8Ws8/s1AM+HyQaf4+x3lvf322G6384FgCKYngwzw3xU5i7foinaKdPyvCUrsXH+fOxzZvU8Cl/fY+oWiP/8VXEXwLMwYwH3hYmp/4KI8XmgdrNgwaYzdOuaZVExfGy9NiR4T1X1XHcnEfklgJ+oanDoSg0fhUkmXStm/ez7YQzp1c4+1wO4AyaJ9z9V9dqSc1HDIyIzh3luj63wLSHDktsPpRDSO1T1t8cuy5j4Fbi9Tykm1WUSDeuCTWds+6y2TH2XR0R+DcDXVbVpxT7rSamE6/Qc0mTMPmPHNWtWhiwLtZsHVtM25lbdfxvjItKBh2AAACAASURBVBvlvZNbecj0UJlh9SuCOrNC8qapUU0lF53YwDlkeYqhICcDuGmwi04xQ2k2N9yGFpCebW4CtZsHIc3799/dZ96Sme+N3dAipC9qhwS42dUyhgimJJ4+KvkYHRAAZgH1+2DGQREySVC7A+OazKohLVvXLIuO6znUx8yukj7o7HnRZd3PY6Cqa+AMdJ4U2gSa2ZqJyg1VPWDsMkwyXeuYDa14qN3uSNVxW93n1qtFZjdu72aZp0mcXwIgwrBunD8fx2zenPyDYPdEM1JM6xgmNWYcaw4s2HSGWu36uN8b9Tk+Q+k4ZixrLpRpl+QLG1uEbE+dp0nxPJXLWrk/mKF/PHOXb1H/NWgBRiT0mXP6Pqq0kKMZYODvn1RNxuh43pKVre7dpN531VXO+MTnPsNsioFD0jae9hGPfe1OQsy1SQIye5m7fIu6WXarW//3Yf9O1W5lhnVo8dX96HMYm9MHMwfQz8wCjV1JlQXK3DOtsdplT0B3zFuyEnOXL9vu+2yq4ab6qjOqZdrNSc9l2dVpjYFj8pxeLoo+pq+4XKXdaYm5ZHIpiz/WpNoJgC7uGOyq4QExcS1pDGuTH0xMOjj1fPb/kxi4xzagdaRkpVw9TFJ3a4hJ19WY+A2ulAkibWibffW1OymwkdUPQw9NmWYWbDpDOfFqcgglzULMWwIs2JT+5M7nMq3LSo+NMa2V67Darim/tZ+jKRk6eNdNMmuz1inQ7DvOJRD6ZR+jYlVdpW1128Wi9U3OlTN1XfhAf5nRvslBt0B44krMdzotGhsad63xNnXbkPqtK+dYMRd4zi/QsObNWBPk2mi3s1UC+iAliDQxiKk/6qprjDUBqk9sF0+byR85VaKpPQR1XRgp18otKxZbnjG6PptS1SU56ROYYrTbRmPuQvXAdJuNNpoeu1FlqXuYTw7462DnEvumDfcBIy7uUmm56BZol9ioNayhQD/p3b+WnLrnUwbZ98m0j0NqYlrrzgfE/xZyGHrgT+arWgOyCUNqN1avXTS+huLmwy6SUPYjVrsxXWv+k5SanGMsumxs5RJ3ge21PAla9SnzC5PuFXLDjw8h3S7YBA2NKc2FJrqIXtaqqwtOE02M+9gtnVRDmho0c9JDmXZT71vZ/g9esA7A9teIvcehijdUyXZpHPwWd9PxpmxQjUPZ9+5r09dW6B7bc5XpO7deAWBm933VeDifseOuT4yO3X0mxbxWx9z4SW1NmC3jZnNY63as3qzKMayrzjpNAeD084+rLURo8o0/scF/L0QosORkgmIqjJyCY5sKvonwchkHuOqs0zRGt0BYu/Zvq9mm5KRdlxTTPoaeY3S7esUGnH7+ccnnLtNELtoFwmOwU2kzqaysi7FVgRrgG+26rvBpib0hYnULjHOvYv1CH6bSNXE5m1Zbzqoyusbb//8wpTR0od9Y3+hSpd1Kw3rz+usUAG5ctxYAsGrxHnGl7JE+DUBOwa4tQ4gtRC6Vvg2eQLMGV1eMbVgfvGAd9jlzEYD6xtZsyZpOgmFNSRYMwdg6zj02D6HhSTCsxy5aHO0XujKWvpHLwbC6PQEh3GE/ZUOAhqRP/VbFr04NqxWg5cZ1a7MwrWQmY4ktRC6Vvm1sAXk1uFyaDCWpOyb3ih0Yt1s/90rf4iYLctNtDKFGELWbhtVq3fjrXLTbxi/EmsxYY1eXwUzZ126vWvopB9OZwqQmCZJXCVi69hEA+bT8p5XQWKBUka1esQEASrtN7XaLDTY3rlu77ZiUR0TmYlbLWLr2kax026SCzqFSj9Gh+x03CY5V2vV1C4S120U5cmASY25Ipzlo1yU0yamtRprGXGB77aaWyw5hyiHu2iSB1W6dce3a7JV1vYeuU3ftqkXxY88xJH3GOVe/rnaPPuGO3q7pUmtYrfCOXbQYbrp/9YoN2FhSMXTBNMyYbMvQFeyN69ZuE6EbeCfxueZWp8BMM7N6xQYsBUpNTRNC9ynH76tMT11U1LHXiiFkSOuI0W6O96QM+3mGjLnTRBP9dfEbOP3845L069avQPm47Nix+GObVf/z2L/HmkCWk5FsylA+ILWxZWPUjevW4sZ1qDw2hhjtRo1hLSukS5uCDoXfxTJJhIJg2+/cP2co2NprdD14egjsOEA3eOas3dRAPgk6Lqu8677z1ODZtXaBPIYEuPi6BfLQblMDkrt+22i3ah/3vHY/e4yv+xTtjh1v3WEsvmH1Kft+2vYqlp1naMbW9lB+wb5vtZuSCGqyPm+yYQ1RJkqfOhEDaS3UoYN1Vbn6KEvd99D1NasybakDp4HxA6iLP74KiNctMNnaHVpHMd9BzDXrTGssMdr90E4f3NblnsNC53Wx19dhqCEWsw2YWfHEMM1xt0vtdlG2GO3mFGdjPENq3G0Tc8ds0OXmF7q+blO/ADTXbiPD6qf9Q++F/vb39w1E6Dp1gh2LsnJ13YqpomvxucR04eT4eMA6qoxr19rNhbKydp1pbqvdUCUfU/GXabeusveX3fPfz82whnRXls1y41OZrqvoQ7+pxtgnVN/4525CTHnc7zNWu3WE9Bkz/CpHo+rjxtkqv2Dfj9Fumc5D5+wS97620e40+4VY2mq3cYa1qtJ2jWZMcKyjLFA3ISSc0A8kpWw+KcJo+gNIvU4Vrvjq1h1Nmdmec0Ct63LtQrf2nKlBtatrx9Cm4u9au20q/Ziluew+D16wDrdc/zIcfcId2x3vkoN+/dUu6mJul/hxvE3Mtecrey/VVLeJu01027Vh9WMuEP9QiLrHReegW8vN66/TqnsaqotjKdO8b267TnqV+YWm2p1NfqFMmzHr0wI1htUndohA11QFuhSqMmtNDWuojF0RyjZVbW9CKHCmUvZoz5wCZxm+pts2XKrMaax2hzSsQL+ZiabjAGMIadc8fQzb1p71efCCdTPKtGrxHtvp3gbgHNZzdFFdpbdc/7LKfZoa17qGVdNKP9ZQNKXvuJtLzA3F0rKHKuSmW2D7NVrrermqWHDiQqkywn1oty4DnMps8Qt1BtWnU8Pq4k9oAapbPHbfKpGVvXf0CXfMGM8W28rwr9c2UMaYEpfUGaNDj7dpY1b9J5aF1qVrV7phKctkxWQA6iYWpGq3q8o8li6C59jatWNQLa5hdYPmgk1n6Mb587fdh6rldnLXcF0MjtVuaH9/QoX9fx113bapTHPMBZrN7Pcfr+yTs25VV6m/NFJZlt2lypT6hrRsIlsddcaWfmF7v5Ci3TrdAj0aVhe39dQ0RR86LnVtr67Hl6QSO6O56piUsSFNBzz7NMms1q2pmHPQrCJkAixVk1rc7Sldk6lBqk/qWugxx8Ror8vnUNsuKbvGbllm1ZKyFugkabjJpMLTzz8Obsb2xnXhdUCryEG7sWv1lu3f1eocKecpq/BDPVUuMfqdBN26TyJsMtzPjblNjRv9Qvo5gPJGVlnGP/ZJijG67cSwusRMaqnbFlvx+8sp9Bk867IMscuYlO0bM7Gp6oY3Weeu6RCAUDn8H4PI0uyDZhVuQA1RlhkILe8Rw5gVf4p2Q/uH9JbynPehtFsXMKdJwyH9xkySBer10MUklBjqrhO7NFrZvl3qtux8IUIVfluj6pZlknS76qzTtEndbXVrk1oxMbdsUmdfVHmGIfwC0E/MdbXbRQPLlidGt50bVh8/cMZ2X3U9EDlWqGX7hcQ39hjSWDF09Rzw2MANTFbQrKLOuALh7tMmpDaGUpd9mjbtlo2dTrlWVbZ3WjRsu2FDQ0384VYp+HqKqWhD643GXL+rJc6A9uP222i3K6M6bbq1Go1p8PvDVJrqtqxBlKrDpma5CUPE3LJzp2g35lpNEgS9G1aXqlZ/qIXf5rGKsQKJXeS5yXI7VXQx2akMXyxtz5+amQImM2hWobpqhnZDlbWrhyEeUtFlJZ5C1cz6tsQ+ez7GrDbRrcu0adjiTt7ye7OGeJRtrEHtK+b2oVugXrtdVPgx92SSdWsnU4Viap1uctBuV8bUZ2y/EJtVbRNzszOsLlVPIQoFTxf7ofsQR9fndBfXDQ24L6Nr05lCSkusTICTHDTr8Ie9VFX6Pl0F0r4CYxmxATOknS6160/2i7l+GbPVsPrYbtm6rsEU7Y7VmAqRUtn3FXfrKvwmGakypkG3ZT1b/lAAyxAxdwiaPKDHEtvYT6GukTW0bkczrD5+Bit18kZZBV41biW1hZ9StjK6qGy7NgBNmLYuqSa4jyKMMauWsbOubRd99hlDtynX9ZlNFX8sKc9cb6rflK7eEH3E3SG0GzKrTbQ7G3XrrirgThL0J1vlEHenzS+UmdUukwIuE2VYXVRXaRsB9p19GqvCjz1fLKnXnY0BM4a2jS2X3LTrz9x0W9h96XbekpXYOH9+6Yz/Pit7l9mmYyDOuNqnMfkJA6A+o1r1OMc6msTdIWOub0y7MKrUbXj94S4Na24x1zKUdkPGdAztTqxhJYQQQgghxLLD2AUghBBCCCGkChpWQgghhBCSNTSshBBCCCEka2hYCSGEEEJI1tCwEkIIIYSQrKFhJYQQQgghWUPDSgghhBBCsoaGlRBCCCGEZA0NKyGEEEIIyRoaVkIIIYQQkjU0rIQQQgghJGtoWAkhhBBCSNbQsBJCCCGEkKyhYSWEEEIIIVlDw0oIIYQQQrKGhpUQQgghhGQNDSshhBBCCMkaGlZCCCGEEJI1NKyEEEIIISRraFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYQQQgjJGhpWQgghhBCSNTSshBBCCCEka2hYCSGEEEJI1tCwEkIIIYSQrKFhJYQQQgghWUPDSgghhBBCsqZzwyoiy0TkJhF5WkTWeNteICIfF5EHReQREfl6xXlOFBEVkau8948s3l/fYZkPEZGrReQBEXlYRK4RkUO9fc4QkfuKcn9KRHZ2tlV95sUi8rjzerIo/zFdlZ+0R0R2FpGLRWSriDwmIreIyOuc7SeLyKZi23dF5E0R5/xPIvKsiHy8h/LuKyKfEZF7Ck1+U0SO8/Z5W/F5nhCRL4rI3t7n+Vahx/Xeccd7mn280OxJXX8O0h4R2VtEvlDc560i8rbAPucU93BRxXkGi7mBa79eRG4QkZ8WcfZ/i8huzvadi7j7aLH9vd7xq0Xke8XvbYm37ZOelp8Wkcf6+iwknpq6MynmishficgPint8t4h8tqcyU6sj0UeG9R4AHwDwqcC21QD2BnBY8e8ZNed6AMBvichc571TAXy/aeFEZE7g7T0BfAnAoQD2A/BtAFc7x7wWwJkAFgI4CMBLAZznHF/6mVV1rarual8A3g3gPwDc3PQzkF6YA+AuACcA2APACgCfE5GDROTFAC4H8F4AuwN4P4ArRGTfmnOeAuAnAP7QbeCkUqLZXQHcCOAYmN/SJQD+WUR2LY45HMAqAH8Eo+knAbjG+WEAfwvgAv/EqvoNT7P/L4DHAXy16WcgvfIxAM/A3OfFAD5R3H8AgIgcDOAtAO6NOFfnMTeSPWBi6P4w9cMBAP6Hs/1cAC8HMA/AqwH8pYj8rrP9VpjYul1cVdV3enr+DIDP9/EhSDLBujM15orIqTCxblFxj38dwHU9lZlaHQtV7eUFc0PXOH8fCuBRALtHHn8igLsBfBLAfy3e27F472wA6519PwpjNh4FsBHA8c62cwFcCSP+RwG8I+LaewNQAHOLv68A8CFn+0IA99V95pJzfw3AOX1973x19wJwG4CTABwH4H5v2wMA/nPN8VsAvAvAjwG8xdumAP4cpvHyIEzA26HYtgTANwFcBGMsPxBZ3kcBHFP8/0MArnC2HQxjanbzjnmH+1sqOe+nAXx67PvBV/DevLC4r4c4710G4ALn768A+C8AfghToZedq3XMBfArMI2juc6+xxS/l50SPtfvA/iO8/ePALzG+ft8AP8QOO4GAEtqvq/HAJww9r3ja8Z98f1CUswFsBLA31acfw8AF8M02n5UXG/HYpuNt38H4BEAmwEsTCg7tTrQa8gxrMcB2ArgPDFDAr4T2cV4KUymCgBeC+B2mFaZy40AjoIxmlcA+LyI7OJsfyOMad0TwNqIa74KxpA+VPx9OEyryHIrgP28LEQtIjKvOPelKceR4RGR/QAcAqO3mwBsEpHfE5Edi66pp2EMbdnxx8O0vP8BwOfwnIZd3gyTCVgAo9HTnG3HwZjZfQF8MKK8RwF4HoA7irdmaFZVt6AwNnXn8s77Apjs3CUpx5HBOATAL1XVzYDeCnP/ISJ/AOAZVf1ywjkbx1xVvQ/AegAnO/u+HabC/nlCGV5VXBcishdMNsuPwYcHjqvjJBjjUzocjWRBasz9dwCniMj7ReTXRWRHb/slAH4B4GUAjgbwGpjGusXG230AnAPgKncIVQ3U6kAMaVgPAHAETAtmfwDLAFwiIodVHaSq3wKwt5gxpacgYPZU9XJVfUhVf6GqHwGwM0xG1/JvqvpFVX1WVZ+qup6IHADTxeaOO9m1KLfF/n83pHEKgG+o6g8SjyMDIiI7wTRsLlHVzar6SxjdXQETNK8AsFRVn6g4zakAvqKqPyn2f12gO+tCVX1YVe+E6Z5/q7PtHlX9u0LTdZrdHSardp6qWm36mkXxd6pmT4LJAF+feBwZhtL7XAwP+RCA96ScsIOYewmMSUVhHN4Ko88oROR3YH4/Zxdv7ep8Ljj/T9UyivNeqkUKi+RJasxV1csB/BlMA+t6APeLyJnAtuTD6wC8R1WfUNX7YXqv/tA5xf0wGdqfq+pnAXwPwOvrykmtDsuQhvUpAD+H6d58RlWvh+kef42IHOgONA4cexmMwX01gC/4G0XkfcXg7EdE5Kcw6f99nF3u8vZ3BzUf6Lz/IgDXAvi4qn7GOeRxmHE0Fvv/1MHQp4CZqqwRkR1g9PYMjOYgZqLKh2G6TJ8HM87170XkqJB2ReT5AP4ARTZfVf8NwJ0A/Mkwri63wjTkQtsgIrc71zneef/5AP4JwL+r6n93DvE1i+LvVM0yaOZN1X0+D8BloQZyzzH3agCvEJGXAvgdAI+o6rdjPoyI/CaMOXmLkzW25fNjcJKWReQlML9d9nBlTmrMBbbNF1kE05P6TgB/U8w/mQdgJwD3FhOlfgozvt9NIPzIi3FbAewvMyeg3u6VkVodmCENa2n3qareqTMHGvtcBjNI+cuq+qS7oai8l8N0Qe2lqnvCtGjEvYR3vV2d153FefaCMatfUlW/C/Z2AEc6fx8J4MfOkIFaROSVMIbkythjyLCIiMCMc9oPwElOF+ZRAL6uqjcVWfobAWyAGQ8Y0u6bYYLUx4tZovcBeDG2HxbwEuf/B2Jmt6uv2cOd63yjKO/OAL4IM2ZqqXfuGZotzMPOSJg8UwTNE8GgmTPfBzBHRF7uvHckzP1fCODPHQ2+BGYi4fI+Y66q/gxmGMximIkwUdlVETkaZvLraaq6bcJM0UtxL7aPwbcjjVMAfEtV/yPxODI8qTF3G0WW9PMwnuMImMb/0wD2UdU9i9fuqup207+4iP+WA2F6udwJqO5ERmp1BPpY1mpOMX50RwA7isguxSznr8Nkmf5bsc8rYSrDa+rOWWQITgDw14HNu8GMTXkAJnCfje0zDnVl3r0oxzdV9czALpcC+BMReUVhbM8CsMY5vuwzu5wK4B9VlUtU5MsnYGZ9vsHrhr8RwPHFOFEbrI5HeSPsVJhZr78KE3iPAvBKAEeJyK86+71fRPYqjOFfAIhehqUYtnAlTM/FKar6rLfLWgBvKDIELwTwNwCusvorxoXtArM6wg6FZnfyzvFHMEFzS2y5yLAUXaRXwWSTXljE1TfCmMSFMBW21eA9MA2bj0Wct23MvRRmMsvvwUx4rUREjoBZheLPVPWfArtcCuCs4vcyH8CfYmYMfl6hZwGwU6Fnv347xT2GjE9F3ZkUc0VkiZjlpnYTkR3ELEl4OIANqnovTDLqIyKye7H9YBE5wTnFvjCNu53EjPs+DEBw3De1Oh6h5XK2Y+7yLTpvyUoAwNY1y7a9/9CFB0tg97NgBi1b3g4ztu5cEXkjgL+HWSJqK0xFuzmmDKp6Q8mma2BmwX4fwBMwY1PuKtm3jDcDOBbA4TJzXbRXFK25r4rIh2GGMDwfwD9i5mcMfmaYFQpQiPNkmPGAgzF3+RYFSu9T5TEh5i1ZOeP+x+yTcu0+WLDpDLXlqSqLmAlxS2Fa4vc5je2lqrpWRM4FcKWY8VAPwKwacW3gPC+GMQpHFxNQLPeJyFdhzOz/V7x3NcwM6z1ggtPFCR/tt2CWm3oKwE+d8r6uyArcLiLvhDGucwGsA/DHzvF/BDP73/IUzHCVJc57p2Dmci2DkaLdKs1aJlm7obJ7ZXs3TAPpfgAPAXiXqm6X0RGRXwL4iaqGhgBsR5uYq6rfFJFnAdysqj+MuNz7ALwIwMUiYn8HW52s1jkwDcqtMFq9UFXdZdauhTHYgPltrIYZzrAeAETkP8PMo+h1iSBXi031454jVrdAVN08GNYzbF2zrK4sVX7hXETE3IJHAfwVTONoRxidvMvR8CkwS/h9F6bB9R8ALnSO3wCzFNWDeG5Vl7Ie1KnQqk8Tv+Af69JEu3XXlpihaQs2naHuSV3G/mGQcnwRVd2rvir9uuv2TajSz0GzIqIAXq6qd9TuPAuJ0eNDFx4ssboFwvHL3y+XSt+NuWW/uxx0XIWI/CvM0mp/P3ZZhqTLRMEkGlarXQAxhnV0iiTVO1T1t8cuy5g0aWx1keAC4g1rbYbVtJSqt+cuSGKIqdynCT9wkumiSz3bwJkDfswt027OsVdEjsVzy7URh1TdxlT8OWLLnLNOSZi6ZFeXsTdF25WGtWwoQGg/CjIvZps5JdPBWLq1mcxJo003Xl+IyCUA3gTgL9wx+yLySRTLXXlcrqrvHKp8feLqdzbGYDdJQKaHJo2sPogaw0omh76C5CRW5jmjqtkYjDHJpVLPSd9NypKTcVXVU0vefyfMckNTyZhazj0Dm3NSS1XXYJZPcMolDtfppDPDmlPAnFSafoe5iC0n2NLPi9w1mluGdeP8+Thmc9R8VNKCLiZJ+efpipz0GMPc5Vv0wQuo2yFJbQiE9N5nbO56OEulYZ23ZGVy4OwqAMxmQiLMpcKvGyKSU0uflf6wuLrNTa/2/1X6zEW7G+fPb3xszpmsHImduAeUj+sbQutV2s1Ft+a3Vq5darMfYr7XkEaHitFdarfSsLYJnAAF2gRbwc5dvmz0Cr+qhR/bchrr/ldpl7rsHjvefe7yZdl8tyH9TsoEljYNLeo7nSrtllXsQw+/yl27MX6BCa1uyckvAP1rN3pIQOzSMD5Ng2dZMJh0keeSefJp0v2UewC1hMrJSr0/xtB4in4nRbfAeGWdrcZiaO2mxt1J0y5QvcqF/f9s0lgXPPfdjTeUaQztZjfpqm7R5Ek0GmOaVP877FrcOQdQN0tVpiX/mK60NYk6bcNDFx4sCzZB+9RDE+2WDQupK+dsunf+78B8z8u22z5p30mdIQp97ly0WzWcKeeY25QuNDapOi0jp+RWrH5tlr1Ou00ToLWG1b9wHz+W1MVnUxbE75OqH8iQYqsTU5Ng6VPVRemeP+dAGqPdroKe+U4uanOKiWBmA7P5ebpoSJV1SVaZVpdctdskuDeNTVX3IQdDkGpC697vgrbardItEI69k6JdIE2/sRore7KS2bZsdJ02ZcwMf5OJp11oN4VGGdZUg1Imwq5uTt15unhqQ9W5uvgc/k3sO81vK/G245QnjdgGV9un1bQxb7lTpvfYSrNPc1q3/yRPwmuSLOijAhwrYeBfd+jFy8fQrXvMJGi3LrERq98u7m0ODSyXvsxoV3G3bntT7Xap21rDWnfBFPM6Voq7y+v6QxaA5q3bIceehMQWK8BJCJQhYjNrQPo4q5gsVc7DAuyyXzcfdlFt+ep+P33quOsG1aSvHJFjZq1PnVdl0rr47H1pt49EQO7aXb1iAzaef1y0XwCeu4d9DnMY27ja63epW8sk6LdL3VYa1tUrNpgL1ojQ0kUwnbS154C0Z+aORZtsatMW/phmra12fbauWYYFm85Q/z7XraSQU/dUKAOcsoTaEMFxNmb9fWzFD9T/5mJ02+b4WFKTAinLSbk0LW8f2p3tOvU5dtFirF6xNjrmAtsvO2eJGeeYanKbrlnaZHhCFW3MeZc6nkT9imr5d73qrNNmbFy1eI/eCxTLmJmFsc1nDH2JsS4QufdlTJN28/rrFABuXLcWQB7ardOsX4l3OZTFpS7DPEkZ01hiKlD7HYzduHC1e3pC5d8HOWRwXeqSIpOQcWpClQZC92gMDa866zQ9dtFiAEa7Y8bcPnTbdt3dut7ovhpVOSQByvSbqt1aw+oKEMDoATSGqhZZ1XZ3n0lhaCFOmmEF8jKtsZRVxEONtasjdThJ3zpdvWIDTi+yknVlCZGLboHttTsJMddSVhGnGIghxu/HMFRF36V2LWNo2NUt8Jx2gXyGlTXRYW49xUPV+baX0sXVaYxuLaMY1mMXLc5ShNNMW3HWicoVpXuvy45JCZy5GFZgshpcuTJk48jqcgzt5mRYAWq3Dama7cKkxlTkVrtWt0B3cRfII8Nq/QIw83dMDcczpjmtw9Xu0SfcUblv6J6XNQSqtFs76coKzi2cOz6wrDAp+DdlWgU9dDa0ToT+dnuv3W2xLajccD/LsYsWbwueq1dswNLi/dluAMoGw4/dfRRDl9rNrds7FHMBbNPubEwY+OPoQ9rtQrdD9AS4uLq121Ni7iRpFzDarfuOp03XOcfTMo8Q0mCVdm9cFz6u63tZmWEFnhvHWtZysrQxNu4X0VUw7mopkBzEVmU8u/je/XPYoOkGz5ShAJYcMqxWr+7/fe0Czb5H/76EzjF08I3Ra2ql2IaUgNjkvH4XVdU1Ylv5uWVYGVkLpAAAFXpJREFUXbqKu6GY6zOUdnOIsSGG1q6Nuf62lO5USw5DAly6irlAnHbbUvad56pVl1S/ENvAr/ILZdfoWru1htVSJUYgLMgQZYbB4v5oy7ZXUTZmrm0mKSZl3uWPJzZF3+U1q7rDJqFLyqVOr5YU3dbtX6XdMTPVOWq36++D2i0nJuYC/Rm0NtRpqeuyDa3dproFJsuw+sRolzG3+2vWXTd1XGqVh2piVoGODCuw/aoB/vgb+57fHetvt3+HMl9+V4J77JiEym3/dumq9R1D14EzldwCpktZ8AxprExfVfc8pNPQdXLSrk9fvSQxdKXdplnAHI2qpSvthmKuv83dXnfOsan6Pc127frkZFhDwwTKtJwaZ6uuNxYxdYPLUNptk131KRuOk5ocCNGZYbWortLVKzZUmssqA5pKSqCuoqzFFgruqeXziRVAk8HOqdeowxVZ1ePZUsZLjV35p1T67jagmQbKqGvgxFB1/FjabaPb2GvU4esWKJ+d20WwHIpU7YY0EUtZsiAUx7s0AjHZs7JkR9VxfWu3r5gLpD3QBKhey3ZMHbuTrywpjfy2xPiFKm21uW7MuegXypm3ZGXlw2waGVbLzeuv0zoDkBo8gfoWf9PgWRaQu6CPVl3KWJGmhAJnF4xd8YeCZixNAmdMpdpEt139tqqYTdqtW7B+bN36hFYMqLrvTWMusP2s37aNkirKei9S6Vq7uerWUqVf1zjEPMFuKOp8gqWrmOZqN6Zh5O43FLMp5lrK1rH1DW9vhtXiDhWI7dIvM6UhgbljrOrGW/nXKLt2F6akTnR143FD+w9Jl2Y1RXRjEMpa+ZmjsuEqPmWm1D8+tNRHrG5D125TyZd1B4fIXbdAWLtuQIx9klLuurX4DbAyoxkaauX/3983pPey7sM6bdTFxrLYW5c9nsaYC7SLu1VrfeaqYyA8IdYl1i+UvWfjblO/0HUiIFa7qUyCXwjF4bo1ans3rEDcMIEqYir7OvrMCMRQJqCUSWSp45pCg55Tz+FW+LFP8WgjulyoGpNtqTN6p59/HG65/mXb9h1Su1120w6t3aYTTXz89VObmFSXSdCt/xQ3oFy7deM+rXbt/l2slhFDV9qNXae36rgxY24qVdp1yyWyNHsdW8rWHfaHo9QNf4qNuzENryr61u6k+IVUyrTrl6lKu50ZVh9rYGOpyqzG0odhrUuzl+1XdUys6FIe3+Y/jz2Gvp7gMQmVfhm+iY1hkrUbu5RJ2f4pC0L3pd0ujOrG+fMnqpIHmmkV2H7cc5tJIEPotuw6KdqNXXIu9ZGZTZZPTE0QxJpUl0nTsktVj21ZjxYwzDjQWKoyuznEXCBdu136harVBUYxrCF8IYYCZRdrRKaKJNR1EGNA25azi+6h2EfA9dX15GJFOMnB0qfKFHRhVC1lrf42lXLot9SHdsfSbUyFn6LdSdVt1ROGLGUZpbK1P5sQo92qbtoYrY89Dg/IS7vTaFh9yiZxucTEPh9X9ynardrX/y2F1jRvyqT4hbrr1C0nmo1hDaG6alsB/Jvq3qCyDznWE5naCLCviU4ubZ6DHDpPFSkp/WnAH/4SauHHaLdLfD3WPRSiKTlot4vMlGVSu1HLsJW734hyK047bKXugS1DPhbS1eQt17+s0bCaOvrWrqu3NtptY1Qt06DlKtwe3FjtdvnAlLKHP9jfXpeJghxiLtCNdmNiStaG1WXVWadpSHwh+gimfZjfpovnpi5x0hWp2dQQ0x4sQ8Rqt08TEBuQY56X3nQ9vbF0W3V9n9mSlbL/r1p31W9sDa3dVBMxjdrtIub6TJOWy7AaP/38+kdsb5w/v3PD6pvVunO30W6qbuuO6ZIuGlguE2NY3WxrzLiKNo+hrBNY1fY2jyxsMtYk9hxNSLku0L6FNM2ortJYbbjfY1nXfZsJMDlqt+sA2rV2p1G3ZUNarGl1M5mpMbfrR/z2FXOBvLTbR8z1mUYt+zT1C1UxsukwqtniF1Kv3bV2szSsTSp9l9jWTtU5Yo5NoUkLqc05Y+grcM6GYBmiqWF16aqlXnZcE7oImjHnS4FmNZ3QZJZUwwqUd7cCcY/O3jh/frKGp0W7Q5hVYHboGUjzDG1irr0G/UI8fWg3K8NKCCGEEEKIzw5jF4AQQgghhJAqaFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYQQQgjJGhpWQgghhBCSNTSshBBCCCEka2hYCSGEEEJI1tCwEkIIIYSQrKFhJYQQQgghWUPDSgghhBBCsoaGlRBCCCGEZA0NKyGEEEIIyRoaVkIIIYQQkjU0rIQQQgghJGtoWAkhhBBCSNbQsBJCCCGEkKyhYSWEEEIIIVlDw0oIIYQQQrKGhpUQQgghhGQNDSshhBBCCMkaGlZCCCGEEJI1NKyEEEIIISRraFgJIYQQQkjW0LASQgghhJCsoWElhBBCCCFZQ8NKCCGEEEKyhoaVEEIIIYRkDQ0rIYQQQgjJms4Nq4gsE5GbRORpEVnjbVsoIptF5EkR+ZqIzKs4z0EioiIyp4cyHiIiV4vIAyLysIhcIyKHevucISL3icgjIvIpEdk55jMW208WkU0i8piIfFdE3tT1ZyDdIyJ7i8gXROQJEdkqIm9ztr1DRO4QkcdF5Ksisn/FeU4stHuV9/6RxfvrOywztUwIIWTq6SPDeg+ADwD4lPumiOwD4CoAKwDsDeAmAJ/t4fox7AngSwAOBbAfgG8DuNpuFJHXAjgTwEIABwF4KYDznOODn7E49sUALgfwXgC7A3g/gCtEZN8ePgfplo8BeAZGE4sBfEJEDheREwB8CMAbYbT7AwCfqTnXAwB+S0TmOu+dCuD7TQtX0nijlgkhhEw9nRtWVb1KVb8I4CFv0+8DuF1VP6+qPwNwLoAjRWR+zHlF5PUicouIPCoid4nIuc42m409VUTuFJEHReSvK8r4bVW9WFUfVtWfA7gIwKGOuTgVwMWqeruq/gTA+QCWRHxGADgAwE9V9Stq+GcATwA4OOZzknEQkRcCOAnAClV9XFVvgDGCfwTgDQA+X+jhGRg9vEpEqu7pMwC+COAPi/PvCOBkAGu963600POjIrJRRI53tp0rIleKyOUi8igcDVqoZUIIIbOBIcewHg7gVvuHqj4BYEvxfgxPADgFJqP0egDvCnRP/jZMpmkhgLNF5LDIc78KwH2qaivtGWUt/r+fly0r4yYAm0Tk90Rkx6KMTwO4LbIsZBwOAfBLVXUzoLfCaEGKl8X+/4iac14Ko1kAeC2A22Eymi43AjgKJnN7BYDPi8guzvY3ArgSRvdrUQ+1TAghZOoY0rDuCuAR771HAOwWc7CqrlfV76jqs6p6G0yX7Anebuep6lOqeitMxXxk3XlF5ACYruD3VpTV/r+2rKr6SxijcgVM5X4FgKWFQSf5UqXPLwM4WUR+TUSeD+BsAArgBVUnVNVvAdi7GFN6Cowu/H0uV9WHVPUXqvoRADvDNLos/6aqXyx0/1TV9ahlQggh08qQhvVxmHFwLrsDeExEDiwmszwuIo+HDhaR44qJWg+IyCMA3glgH2+3+5z/PwlTWcM9t4gc6JzzRQCuBfBxVXXHJPpltf9/rO5DisgiAB8GcCKA58GY6r8XkaPqjiWjUqpPVb0OwDkA/hHAVgA/hNHC3RHavQzAMgCvBvAFf6OIvK+Y1PSIiPwUwB6Yqeu7vP2pZUIIIbOOIQ3r7XAynsWYwYNhxrXeqaq72lfJ8VfAjCl8iaruAeCTmNlNW4p7blW9s7j+XjAV/JdU9YNVZS3+/2Onm7WKowB8XVVvKrJiNwLYAGBRTFnJaHwfwBwRebnz3pEwWoCqfkxVX66q+8IY1zkA/m+Edi8D8G4AX1bVJ90NxXjV5TBjW/dS1T1hMqCurtU9hlomhBAyG+ljWas5xRi8HQHsKCK7FLObvwDgCBE5qdh+NoDbVHVz5Kl3A/Cwqv5MRH4DwNvqDqgo4+4ArgHwTVU9M7DLpQD+REReUZiBswCsifiMgBmTeLzNQonI0QCOB8f9ZU3RzX0VgL8RkReKyCthxo9eVtzfI8RwIIDVAD5aTGKqO+8PYDKToUmAuwH4BcyKAnNE5Gxsn+WthFomhBAyG4ha43TBpjMUALauWbbtvYcuPLgsu3kWTPep5e0wY0vPFZGTAKyEWSpnA4oZ1DXYDNO7AXxERFYCuB7A52AmojThzQCOBXC4iCxx3n9FkTH7qoh8GMDXADwfJqPmfqbgZwRwrqpeX6xgcKWI7AdjRj6kqtc2LGtj5i7fokDlvYo63jJvycoZGggxb8lKAM9ppem1u8Dq1rJ1zbK68rwbZnmn+2Fmzb9LVW8XkT1hMvwHw3SlfxpmebYoihUHQlwD4Csw2d0nYGb431WybxlTqeW22vXPNW/JStx82EWl54rZhxBCyHiIqlbuYAM5EG1YuymYyK/BdEc2NaWzHtdwptwv36hamhjW1Gt3RUi3Y5pnkkaXhtU2XGhYCSFkcqkcElBmXPpGRHaAGdd30xjXn81UmdVY6kwtIbG0jUGxx7v6nrt8i44V+wghhISpHcMaypgBvZvZOwH8LoAzerzGrMJWwlX3bZoq6XlLVmLj/KhnUpApo0zrZfp235+m3wAhhEwTUWNYh0ZVDxi7DJPObDGmZbhm1Q5lmLt8i3JYQD+M/d3Gano2aJ8QQqaRVoa16RhJEqaLcXt9V8gx41jJ7MJqri/T6p+3D41T04QQkjeVk64WbDpDN86fj2M2b44O6DSuzel6klQXlI1dTangx9CE6ioFgGM2P7dqGidf9UeXjdcu9DypuiWEEBKmMsPaZAwgs67tsBXt3OXLSrNVfRnUlIlVk5xpHbv7elppq4kmuk7RrN1/UnVLCCGzmcoMa1WWKgaagnhsZV02ya1LUiv5KmLKOVaG1dWtZezltqaRkNGM+W67WJGiKbnqlhBCSJheJ10xkzWTIcaXDs2kZazc8vY55tJeZ5r1X/VgiRStU7eEEELqSDasfQX6UAU36Qu/NzGoW9csi6rAx6jky3DLkpMJsOOvfVxd9TGEJad7E6LPsdJ18SGn74amlRBCJofkIQGW1EAfqhjrKsDY7vGhjWwfmdKcKvI2hO5VLpOuyuh6mICbYQ0xVsOrTLdtfpvTQG73iRBCyPZEZVhDmarUsZax2ZkmpHTtNskuxYy1SzHwOVb2VRPsYkxfrpRlWbvG1Ujd/a36LbTVZCyp3fc5arYrmGklhJD8qcyw3rz+Oj36hDu2/R1b8fvBv6pCiKkI20yQGGpGfaiMuVbybZ4AVaeBXDKsq846TQHg9POP2/ZeVdlTsqxVmsp1eIRPVYMzV932jf9dMMNKCCH5UGlYV511mh67aDFuXLfW/L14j6HKVctQZmDSK+8+Hk8aa/wsY1T8N6+/bpuw3UYXkJYx7svQjW1mXcM66RrvChpWQgjJlyjDCmCbaT39/OMmoovYVsShCnkas0p9GNMyUgzrWJW+a1hvXLd2RqbVJTctV/VEtDG5k67vIchFu4QQQrYn2rBa3Mo/t8p+GunTiK5esSH4fpm5s3TVtd4nrmEFZja4qshN036Dy36/ZY0x0hwaVkIIyZdaw2r/H8q0+uRW2efMkBnRMsoMqyU1K5nT8mN+Y8vq1lJnXGOg3qcLGlZCCMmX2lUC/Ayrfe+W67cfG9inCRvDHORgKi2+uYw1XKtXbCjd9/TzjwturzOyZeRUwZfpFjDm1X7GlO/RZ6nz/0kZKtMHod/JpH0XY48pJoQQUk1lhhXYvmu1DD/zWmZ63ElcIbrIfA1Fqulpcu46urz2MZs315r0sjV5czKrQLpugecMfAhft/7fkzJGNoR7z/3yDj1pr46qsqaS2/rOhBBCyqk1rEB592rVezF/h8bHhjJjMcagD6x5CZXJpavu5iaZzTG6tnNZCaCKkGEN6da+bw1o3b0uOycQdy+GMLBVJnP1ig3bPmuXv6WmDayU76NtJreLh50QQggZhyjDCoQnsZRV7n4mKgX/vFUmo4rU67alrDxVpqBp13vMuVNwK/2yST515Fy5+ysGxDRAUvRT1YDzs7ZDNLiqdOU3IEPUlbGtbv1rxJrOMiNeN6a6KTlrmhBCZhvRhtVSlrnyK+smGasyc3r0CXc0qiSHNK11BroJdWNLu86uNq3gJ6Fij9Ft3ftlNG1U9UFdI3KoMqZqN8a01mVYuxqHahtsNx92Ufa6JoSQ2ULUo1ldFpy4UPyMlftvqMu/rFJ0xwEeu2jxjL/98zWdGFRWUTehbJhDGVVjIsv2b7M9lVBl/9CFB0tfTwcbkwUnLhQgzbiGtoV0naJbS9ssaNl+VeNsq87parVr3cbs0+QRunb/LowqlwcjhJC8Sc6wWsrWuQSwXXY11rC6+zYxZ110VbahiTlpM6YPGH5cn8+8JSsnLhMVanAB4UZNTAayTOdlqzCEGFO7TYetNNVu1RjbqnP4E6660G4Vk6ZrQgiZZpIzrBY/0+qbVFNZx3c3bn/shqQKvy9iMrup5Yup6KueeOQfmzqTu++KPndstrXswRhN6GOd17LzNX3gg/ubitk/5hqx2nU1Y49J1e3MjOqipGPLykIIIWQyaJxhtYQqfYs7ljVEVTdkmwq/C1PpnqvMtLY1JV2Mv7OVb91yP312m056JqpqCawqAxvSdl22MnZCU2xjzd3H1+VQ2gXS9BXSbNW5m14ndM063DKJLJ1oXRNCyDTR2rAC9WteVq11eeyixTMmVZXNHg5VbkNkX/u4Rp8TRXzKzp8yVnU2dJ3evP469XsNgPJGV11DDCjPfnepKT9raq9vl6/yH+7RlhTt+rrx95+3ZGXSkJaY34p9ZG1ZGcoIlYOGlRBC8qETw2q5ef11GrvclaVs7GrdeLi+n0IVYyp8QxL6u4zYLv+YY8qObZL1SmUazKpLbLa17LGvdV3mfeo21gj7DcFQF33X2i079sEL1lWWr4trllF3L2hYCSEkHxqPYQ2x4MSFcuO6taWmNfR0oFuux4y/Y2kydrMr3ArVr1zrjHabitfuU3aOJtlajuebSdVqAl3g67bMZDbJwsaMT63Sa5/aLWOfMxeVZlqrrtunUSWEEJIfnWZYLavOOm3GZCyXqvUgUxcTL6t4YrOjMeeKOTaGKjPZpPIdeykf+31NcxYq9IS3umXbQroL6aWtdu3xKQ23pk/Z6tqohoYGpFy76XVTjeo0a5sQQiaNXgwrIYQQQgghXbHD2AUghBBCCCGkChpWQgghhBCSNTSshBBCCCEka2hYCSGEEEJI1tCwEkIIIYSQrKFhJYQQQgghWfP/A2Xo/nwKMj4OAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 14 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# declare a plot size\n",
    "clevels = [-2,0,2,4,6,8]\n",
    "fig  = plt.figure(figsize=(12,5))\n",
    "for i in range(nf):\n",
    "# small plots of elevation\n",
    "    ax=fig.add_subplot(4,4,i+1)\n",
    "    plt.axis('off')\n",
    "    plt.contourf(xrot,yrot,np.squeeze(dsa[i,:,:].values),clevels,cmap=plt.cm.terrain)\n",
    "    plt.text(100,270,titles[i],fontsize=12)\n",
    "    ax.set_aspect(1.);"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}