Generic convert ESRI shapefile (.shp) to NMSIM trajectory (.trj).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert an ESRI shapefile *(.shp)* to a NMSIM trajectory *(.trj)*\n",
    "<font color=\"salmon\" size=4>**Note:** </font><font size=4>This tool requires use of the NSNSD `geo` Python environment.</font>\n",
    "\n",
    "\n",
    "----\n",
    "\n",
    "<font size=4> Before using this tool, <font color=\"salmon\">you must run</font> `NMSIM Create Base Layers` in ArcGIS.</font>\n",
    "<br><br>\n",
    "\n",
    "----\n",
    "\n",
    "### Expected Input\n",
    "\n",
    "The main input to this script is an ESRI shapefile of aircraft GPS positions (real or abstracted.) As such the input is expected in WGS84, and it must have the following $p \\ (x, y, z)$ columns:\n",
    "\n",
    ">**latitude** (in D.d, WGS84)<br><br>\n",
    ">**longitude** (in D.d, WGS84)<br><br>\n",
    ">**surface_elev_m** (surface elevation in meters above Mean Sea Level) <br> \n",
    "\n",
    "Using a surface elevation input allows aircraft altitudes Above Ground Level (AGL) to be computed quickly as needed. Aircraft altitude is defined relative to the surface as the variable `AGL`, below.\n",
    "\n",
    "### Output\n",
    "\n",
    "This notebook will create a trajectory and save it to the project folder you designated when using `NMSIM Create Base Layers` (i.e., to _\"..\\Input_Data\\03_TRAJECTORY\\somename.trj\"_).\n",
    "\n",
    "----\n",
    "\n",
    "<br><br>\n",
    "\n",
    "## [Step 1] Import and Define"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============= import modules (and IPython 'Magic') ====================\n",
    "\n",
    "# geospatial libraries\n",
    "import pyproj\n",
    "import geopandas as gpd\n",
    "\n",
    "# a convenient file-selector dialog\n",
    "%gui qt\n",
    "from PyQt5.QtWidgets import QFileDialog\n",
    "\n",
    "# other libraries\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import numpy as np\n",
    "import glob\n",
    "import os\n",
    "\n",
    "# ================= define functions ====================\n",
    "\n",
    "def gui_fname(dir=None, guidance=\"Select data file...\"):\n",
    "    \n",
    "    \"\"\"Select a file via a dialog and return the file name.\"\"\"\n",
    "    \n",
    "    if dir is None: dir ='./'\n",
    "    fname = QFileDialog.getOpenFileName(None, guidance, \n",
    "                dir, filter=\"Shape Files (*.shp)\")\n",
    "    return fname[0]\n",
    "\n",
    "def climb_angle(v):\n",
    "    \n",
    "    \"\"\"\n",
    "    compute the 'climb angle' of a vector\n",
    "    A = 𝑛•𝑏=|𝑛||𝑏|𝑠𝑖𝑛(𝜃)\n",
    "    \"\"\"\n",
    "    \n",
    "    # a unit normal vector perpendicular to the xy plane\n",
    "    n = np.array([0,0,1])\n",
    "    \n",
    "    degrees = np.degrees(np.arcsin( np.dot(n, v)/(np.linalg.norm(n)*np.linalg.norm(v))))\n",
    "    return degrees\n",
    "\n",
    "def get_utm_zone(project_dir):\n",
    "    \n",
    "    \"\"\"\n",
    "    Glean the UTM zone that was saved as a suffix to the elevation file's name.\n",
    "    \n",
    "    This assumes that the elevation file was clipped using the NSNSD ArcToolbox \n",
    "    packaged with this tool.\n",
    "    \"\"\"\n",
    "    \n",
    "    # the ArcPro tool will figure out the UTM zone associated with the western extent of the file\n",
    "    # (this is the native UTM zone of the NMSIM project)\n",
    "    elev_filename = os.path.basename(glob.glob(project_dir + os.sep + r\"Input_Data\\01_ELEVATION\\*.flt\")[0])\n",
    "\n",
    "    # this will glean the UTM zone from the filename as an integer\n",
    "    UTM_zone = int(elev_filename.split(\"_\")[-1][3:-4])\n",
    "    \n",
    "    return UTM_zone\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 2] Set global variables\n",
    "\n",
    "><font color=\"green\">**line 3**</font> `project_dir`: (str) path to the NMSIM project folder you created with `NMSIM Create Base Layers` <br><br>\n",
    "><font color=\"green\">**line 7**</font> `velocity`: (float) velocity of the aircraft to be modeled - assumed constant <br><br>\n",
    "<font color=\"green\">**line 10**</font> `AGL`: (str) altitude of the route in feet AGL<br>\n",
    "\n",
    "#### <font color=salmon size=4>You must edit this cell.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ================ Global Variables ======================\n",
    "\n",
    "# path to the NMSIM project folder you created with NMSIM Create Base Layers\n",
    "project_dir = r\"C:\\Users\\DBetchkal\\Desktop\\NMSIM_2014_local\\Data\\BAND002\"\n",
    "\n",
    "# aircraft velocity ...THIS CELL WILL UNDOUBTEDLY BE CHANGED... \n",
    "velocity = 74.7222 # m/s  (for a C-185)\n",
    "\n",
    "# what is the altitude of the flight in feet Above Ground Level?\n",
    "AGL = 800 # ft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 3] Read ESRI shapefile into a `geopandas.GeoDataFrame`\n",
    "<font size=3>To reverse track direction, uncomment <font color=\"green\">**line 11**</font>, `track = track.iloc[::-1]`.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAHMCAYAAABlf4QRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABl1klEQVR4nO3dd3wUdf7H8dcnCRA6SFGkGMCKSI2IvXvIKYrlrCf2A/UsZ/9ZTlHPhmdv6FnPglg5OxZEFMGA9CIioCDSpJeQ8vn9MRNYlk2DJDtJ3s/HY2FnvlM+3+/Mbj77nWbujoiIiIhET0qyAxARERGRxJSoiYiIiESUEjURERGRiFKiJiIiIhJRStREREREIkqJmoiIiEhEKVETERERiSglaoUwsxpmNtDMfjSzqWY23cweMLMaCaZtZWZjw/duZvXiypeaWUYFxX2QmQ2NiWVS+PrRzF4zsw4x095mZovNbIKZTTSz0WbWuQTryDCzi0sYzwgzO660ZUUsb66ZdYwZjm37uWY2I6zLT2b2npkdUJrlR4mZNTKz60ox/dlmdn/McO9wHzhxO2LY0cxeMLOfw/1kkpk9ZmYNw/IXzOyy8H3jcB+bYmaTw+mPiFve9LBeL5jZ/HCamWZ2zzbElmlmr2xH3bapfcP9PzeMfUJYp8sTTD/WzCbEjcsIt8mTceOWxgxv834cboMNZvZQOPynmDh/j/m8TzCzvuF3gJvZn2OWUc/M1phZVsy4gu+SiWb2g5kdWUwcsXWYbmYXliT+IpZX4u+cEizr2YL6hvvQx2Y2x8zGmdl4M7smZtrYekwxs9PD8TXN7AkL/jZMCsvOLGKd78S0e0FbTjCzT8qgPofFbqtSzHezmV0azr8uJr4JZvZSzLILyqaY2Vdmtmcxyy12Hy9i3grZHuE8A83stOJi2hYW871YJtxdrwQv4L/AW0D9cLgGcDFQL8G0lwE3h+89fhpgKZBRRnGlFlM+CDg7PhaCpLw/sBxoG467DRgUM+8lwEcliOEwIKuE8Y4AjittWRHLmwt0LKTt48tOAlYA+yV7f9rGbZ0BLC3F9G8CB8UNfw4MK2KetCLK6gAzgZsK9jsgHbgF2C0cfgG4LHz/OHA/YOFwE6BNzPI6Ap8lmK8hMAfokyCGIvf3ZLRv/HxAM2ANsFPMuL2BecAsoFvcOlcDC4FdE8WxPftx+Hn4ElgE1Iwr2+LzHjNuPPBGzLjzgKzYzzhbfpecQPCdllJEHJvqEG73jcDO27GtDqOE3znFLCcFmB3ux/sAS4DeMeU7AY8UUo+uwHqgKXAtMCTmc1Gv4DNRghi2+hsRji/0s1gebQN8D7Quav74MuBeivkbUZJ9vJD5krI9yuNFzPdbWbzUo5aAme0G9AUudPfVAO6e4+6D3X1NgllOAN4rwXJTwqy/4BfBNzFlvc3sm/BXxGgz6xmOPyz8NfOomX0HHGtmF4e/Ugt6OGJ/4fwZ+DB+3e6e7+5PAZ8QJGSJNCRI5IqMieAPcodw/W+G0+5lZp+G8Uw2s34xyz3UzEZZ0CuTsOfEzE43s+8t6CE708zGhL/ci/v1Xmjbu/vbwFPANeE66pnZ8+GvrSlmdn3M+lua2Vu2uQfyxnD8Fr1+scPh+wfM7Gsz+9XMrjGzM8zs2/CX36kx8+1nZl+GbTnONv+iz7Cgx/WusK4zzeygmHZuFLbzt+H0/wz3nwnh9I3C8bWAbkDBdE2AI4EzgQPNbKe4OvzLzD4vaDszuzHcbhPD+FPCef9w97vcPS9s0w3ufoe7z0rQ5K2A37zgr5H7Mnf/pbht5e4rCf5o7GFm54a/qF82s3HAPmbWK6zrJDP73Mx2DWPeoiehiP0VMzs/rNvEcD/bcXvaN059gkRkfcy4C4CXgBeB8+Omzyb4QXVHgmVtJX4/Lsb5wJ3AZKBPSZZPkNh1MrPG4XA/gj80hRlOkIQ3KcnC3X0KwfdKSwAz2zfcPpPC//cNx8dvz9jhRN85e5jZR+H2nGhm54Xj65jZUDObFo5/Iyac/YEJ7r4BuB541t03fV+6++/uvlXvaFj2A0EC0pZgX/895nOxppDPRJHiP4tmlmZmn5hZlgW9Q8+bWc2Y6RN9TmOX18jMvjCzK63ovzetCH5Q/VrKkEcAbUowXan28VC5bA8zO8CCnrkJYZueEY6PPRpwm5m9bmYfWtCLPcTMuoZtOdtijlTEKmzZoY7h/LPM7CUzs5j1PmNmn1nwff9M7DZOKFkZZ5RfwF8IPsyJynaOLQMaATNjhgvtUSP4BfAj4S9RoHH4f3tgNNAgHN4b+CV8fxiQB+wfs7yVQOvwfS2gTsx8XxQTy1XAh+H724DFwATgF4JfQLuXMKbYX1lpYb1OjRnXJPx/BMEvnRSCRHApm3tjRgDHAdcBnwENC+Zlc6/MHsD8mOXOZfMvqvi231QWM64vMC18fy/BH04DGgBTgWPDsi+Ba2PmaxobY8z4TcNxddsZWAfcFZb1KIg7jPMHoEU43AKYH47PCLdTwTLPAr4J32ewZU9LY4IvptrhcH3CX+EECfpzMdNeCbwUvn8GuC6uDsNi5u0Xt60Ltt0TwEPFfFZeYHPP2JEEf5DHAg8BR8RNOxbYJcF8OwMLgKOAcwl6p9qHZc0JfmV3CIcvAMbE74cUv7/+RNjjRfBrO31b2zecL5fgczOF4I/SlTHLqUHQo9WeoMdiKVArdpuG658HdEkQx1yK2I+L2Badw3mNYD/6MK78NhL3qA0CBgIDCP7ojWXrz3hsj9q5wK/FxLKpDsCBBJ+1muHrF+ComH3ml3B8/Dpjt298WRowDtgzZlvNBPYM2+qz2O0a8/5+4Jzw/TTgxFLU43BgFcHndh/gN4KE+KnilhO3zNi2HMGWn0Vj8+fPCJL9/sV8Tg8j6AHdJWyTU8LxCf/ehO8vBW6NmX8dwf5c8IotK9gGKWFdbymmfhkUs48XMl+5bA+CH4d/jWnTRgm+g24j6P1uCKQCEwk6NWoBdQn+Tm7VQ1fMskeFbVCTYP8/OqZsEsH3UBrwKcX0vqUhiVhhBe7+G8GOV6A3CXqwEs0K/EywE/zHzL4A3g/L/kTwpT4yTLoB0sJf/QCz3H10zLK+AJ43s/eAD9z953D8CRTfsxdft5fcvaDH6RyCxKNrCWKKtQfBF83QTZV1XxZTPtTd84GVZjY9XG7Br8/bCL6oe7v7xnBce+A1M2sJ5AA7mdlO7v573HpL0vax9T0KuMKDT8sqM3sNOMrMvgYOAI6Oib/Y8yni6vabmS0D3gnHjwNamll6uOy2wEcxbenArgRfaGvcvWBf+A54oJB1rSL4Y/RfM/sYeN/DHl+23vbnESTlEHwx/Ae4L6b8VXfPDd8fBzzp7qvCusduu03C/eMfBF+M17v7kNhyd//czNoQfIEeBAw1s/vd/Z5wW6a5+7yYWW6w4NylXOA+d//MzM4FRrn77HCa/YCJ7j4tHH4eeMLM6seFV9T++meC/fz3MM41YX3iq1ia9l3h7l3C5ewMjDazUe6eBRxP8ANidlj+A0Hy8HpMW20wszuAf1F4D3esQr+TYlwQ1tPN7G3gUTNr6e4LSjDvC8ArBIeaXipkmm/NLJ8gCT2xBMt8M+xFaA/0dfeNZrYPsNHdP4NN+8xGgu+Q0tgd2At4PWY71grHTQT2NLPHCRKhD2Lm+zNQWK/+I8AhBD8O9vPNvU1vmtkGgv3jZHdfAawws3bh9AcRtHUvd+9fynrAlp/FFOAaMzuW4G9FY4IkCor+nLYg+LF5jruPCscV9vcGgv352pjhae6eWUh8HSw417IlsIzg+6xY27CPb6EMt8eXwI1mtgsw3N3HFLLKTzzo3cfMJhF872QD2WY2ky3/bhUoatnvetBzi5mND+cfHpYNifkeehE4GXissLbQoc/ExgO72eZDAUU5kS2/wJcQc0jAzNIIsvQl4U6wN0Ey1AmYasEhKQM+dvcuMa+d3X1RuJj4w60nAf9HkOl/GX6oIfjwDSsm3n0JegESGQp0MbNmJYgpVnF/RDbEvM+DLX4gfEfQJrvEjHsNeMLd9yY43JRL8Msk3okUn5jG1tcIEqRY8cPxctnycxIfR3zdNgB42AVPUFcDJsW1ZevwjzoEPTKxy0j4AypcZk/gEYKu/nFm1in8Y3gU4ZeAmWUCHYDnzGwuwR/gdrblCemx+1Rh2+8Hgp7BgvW/FCYnWUDtQmJc7e7D3P06gi/nghN6E+2b94RtkenuDxcRW3HbqGC6wvbXkiQ5JW7fBPP9RrAfHx6OOp/gj9vcsP27sPXhTwiSzoI/LsUp6nNLeOjkTOD8cJ3TCXr2+pVg2YQ/9jYCFxF8/hI5IGzXP7n7uBIs9hR33wM4gyD53ZHCt6dT/GctlhH00MRu7wx3fyesy14E2+soYKKZpVtwIdXimATnB4J2DQJwvzzcv2sSJDix9eji7oe4+/CY6Te4+6fufitwCpv39dKK3d/PJEg0Dnb3fQh6tQvaoaj9eDnBj4zeMfEl/HtjwYVAGe4+sYTxTQvbpSVBr9eTRU++hdLs4+WyPdz9IYIfT0sIErg7C1l/wu/ymOGtvpeLWXax84eK/Y5TopaAB8e2hwFPF/xyN7NUM7vCYq7oDL8cexB0cRYYDvwtZvhi4Dt3XxcmQLXd/WPgBoJDmO0Iuj57mdneMcvelwTCxK+du49193vCebuGv+prufucQuZLCXsvelH4B+1wgl9My4qJaRVB8llgBpBrW56TVaLzV4CPCQ65fBizrkYEJ5dD0EtQK0F9ErV9/DQnhMv+dzhqOHChBeoDpxMcIllDcO7RVTHzNg3fzib88gi/6LuUsF6xviVI/Av+kBecp1NcArEKqBNuc8KYm7n7V+7+T4I/3B0JkovJ7l7wy/t84N7wD1eGu2cA/yRxsgDwP2BAzL5esO1eBZqZ2Q1mlhqWGYUkaWZ2tJk1iJmuK5u3Y0l6exMZTfDjoeA8zH7ADzE9XQWK2l//B5xT0BtswbmKtdj29o2vd32gO/CjmbUg+KPUNqbtWwOZYW/jJmFieDNwe1ENkGA/TuREYIa7t4pZ7zEEPasldQPBIfKEParbyoOe9k/D5c8AahV8FsL/axAcoptD8IOicbj/xJ7vE/+dMxNYZ2Z/LRhhZnuaWQMLzr/Kc/d3CT7TzYAd2HofvA+42Mx6xSyjFlsmBQmZ2cG25dGFbmze17dHI4IEdHWYUMUmf4V9TiFICk4A9jKzh8PvuML+3vQGPiptYB4c8RhAcJ50lxLOU6J9PFQu28PMdnf32e7+NPAwMT8+t9d2LPtUM6sbfvecTdAzVygd+ixcP4I/buMs6JpPITjM1sSCQxxdCM6vGBnTbQ3BuUEPh12n+cCvQMGXSWvgmXDjpBF8WL5z93wzO5ugi7o2wS+IbwhOsI6XCrxgwUnOBcu/geC8ukS9ad+G+UA6QU/hgXHJ3DlmdhRBVp9DcJ5ZPjCriJgmATPNbArBH4dTwj8mj5nZrWFcg4CXC2/ezdz9CwsOeQ0Lk70rgXfNbAHwFUHiGC9R20PQLZ5N0Ns4jeCQ6ndh2R0E3cuTw+GXwy8xCD4sj1twEUQeQZJyb/gaakGv5SSCX32l4u7LzawPcL8Ft02oSXBY4vhi5vvDgttPTDaz5QTb+K1we6QQbM+3Cb4ECy4KSCdIQA+MW9yrwCQzuyLBql4i+LX8nZnlAqvN7BB3X2tmhxAcKvrJzFYQnDD/PUGCHa8T8O+YBHQWcFmYvLXz4OTfUnH3JeEf41fDz80Sgm21aZJwukL3V3f/yszuBj6z4NBdNnC8uy8qbfvGaGSbb72RDrzu7u9ZcIHKR7GJpAeHgN4lOLdri8OK7v5mOE/8odyi9uNEziPoOY1d9ujwB9oh7j6yiHk3TU+QGJeHGwlOB7iX4DDPI2ZWF1hL0EOyEVhgZg+E080h2M8KEu9E3znHAw+Z2bUE34uLCLbhPsA94W6YCtzt7r+F31GbEh93n2jBhUF3WHAriSUEvYp3EZzvVJSMsA41Cb4vFrPlfrmtXgJOMLOpBOdtfs3mH0YJP6cx9dkYfn/+FxhM8IN8q783wBXA03Hr7WBb3krmN3fvHTcN4WdmEMHfxr4lqVAR+3j8dOW1PS4PfxBsJPjs/70kcRfGzAYStM9T27HskcC7BBdmjCTYXoWv070kRxUkETN7iuC49FsRiOVj4KYSHpKo9KLU9slmZtOAw9x9cbJjScSCex0d6O7b9QWZYLl/Ac53917FTrx964l0+0rxwp7OT9y9U7JjSaYwkfmR4LYZ8T9ypQKY2QsEF2gUek7aVvMoURORysbM/kZwu4pL3f3TZMcjIlISStRERKqw8NSCkxIUHVPRPX5RiiUqLDgPONEd6c919wkVHE65MrNhbH1PtV/cvdD795lZc4LzFeO97e4DyzK+qkSJmoiIiEhE6apPERERkYhSoiYiIiISUUrURERERCJKiZqIiIhIRClRExEREYkoJWoiIiIiEaVETURERCSilKiJiIiIRJQSNREREZGIUqImIiIiElFK1EREREQiSomaiIiISESlJTuAqGnatKlnZGQkOwwREZEKMW7cuKXu3izZcUhiStTiZGRkkJWVlewwREREKoSZzUt2DFI4HfoUERERiSglaiIiIiIRpURNREREJKJ0jpqIiERCTk4O8+fPZ8OGDckOpUpKT0+nVatW1KhRI9mhSCkoURMRkUiYP38+9evXJyMjAzNLdjhViruzbNky5s+fT9u2bZMdjpSCDn2KiEgkbNiwgSZNmihJKwdmRpMmTdRbWQkpURMRkchQklZ+1LaVkxI1ERERkYhSoiYiIhJatGgRZ555Ju3ataN79+7sv//+vPPOOwCMGDGChg0b0rVrV/baay9uv/12ANatW8dZZ53FPvvsQ8eOHTnooINYs2bNdsdy66238tlnnxU5zYgRI/j222+3e10SXbqYQEREhOCE+xNPPJF+/frx6quvAjBv3jyGDRu2aZqDDz6Y999/n7Vr19KlSxeOO+44Pv30U3bccUcmT54MwMyZM8vkysqBAwcWO82IESOoV68eBxxwwHavT6JJPWrlyPOWkb9hJO4bkx2KiIgU44svvqBmzZr0799/07hddtmFv//971tNW7duXbp3787s2bNZuHAhLVu23FS2xx57UKtWra3mqVevHldffTXdunXjyCOPZMmSJQBMmDCBnj170qlTJ/r27cvy5csBOPfcc3nzzTeB4PGG//znP+nWrRv77LMPM2bMYO7cuTz11FM8+OCDdOnSha+//pqhQ4fSsWNHOnfuzCGHHFKm7SPJoUStHPmKy2DFhfiijuQvPpT81Y+Rn78q2WGJiFQJ+fn5TB4/l+wNOWWyvKlTp9KtW7cSTbts2TK+++479t57b84//3zuvfde9t9/f26++WZmzZqVcJ61a9fSrVs3xo8fz6GHHrrp0Ok555zDvffey6RJk9hnn302jY/XtGlTxo8fz4ABAxg0aBAZGRn079+fq666igkTJnDwwQczcOBAPvnkEyZOnLhFT6BUXkrUylPtM4D04H3+Qlj7CPmLu5H7eya5K24mP/eXpIYnIlKZjfpiOtdc/AIvPPlFuSz/0ksvpXPnzuy7776bxn399dd07dqVY445hhtuuIG9996bLl268PPPP3Pttdfyxx9/sO+++zJ9+vStlpeSksJpp50GwNlnn82oUaNYuXIlK1as4NBDDwWgX79+jBw5MmE8J510EgDdu3dn7ty5Cac58MADOffcc3nmmWfIy8vbnupLROgctXKUUqcP1OlDfs5sWPMU+dmfA6uAFbDhdXLWv0o+6ViNA0mrfwmpNbvo8mkRkRLas2Mr9jt4Nw46Yq8yWd7ee+/NW2+9tWn48ccfZ+nSpWRmZm4aV3COWrx69epx0kkncdJJJ5GSksKHH37IXnsVHVdpv+8LDqempqaSm5ubcJqnnnqKMWPG8MEHH9ClSxcmTJhAkyZNSrUeiRb1qFWAlBrtSWl8P2k7jSel2WiocxH5NMMBWI/nfMbGP/qy5vf2rF7Sl43rPsJdv4RERIrSfKeGDHzwLPbu3KZMlnfEEUewYcMGnnzyyU3j1q1bV+x833zzzabzyjZu3Mi0adPYZZddtpouPz9/0zlnr776KgcddBANGzakcePGfP311wC8/PLLm3rXSqJ+/fqsXr160/Ds2bPZb7/9GDhwIE2bNuXXX38t8bIkmtSjVsFSUpuR0uB6aHA9+fnryF07hLx1/8XzZ4Pnkp+TxboVWeSvAFLak173HNLrnkFKSu1khy4iUqWZGe+++y5XXXUV9913H82aNaNu3brce++9Rc43e/ZsBgwYgLuTn5/Pn//8Z04++eStpqtbty5Tp06le/fuNGzYkCFDhgDw4osv0r9/f9atW0e7du14/vnnSxzz8ccfzymnnMJ7773Ho48+yoMPPsisWbNwd4488kg6d+5cukaQyDF3T3YMkZKZmelZWVkVvl73fHI2fMKGNYPJz5lAPjkUbJlczyfPmlOzzsk0qj+AtNSmFR6fiEh5mz59erGHCyuzevXqlcn91bZHojY2s3HunlnILJJk6lGLCLMUatY+lpq1j8Xdydk4nvVrHmNj9kic9eT7IlatfoJFq59kvTehTnov2jQaQL2a7ZIduoiIiJQTJWoRZGbUrNWdmrWC7u/c3HksX/U469f9j7z8VeT7Spaue5v5a/7HBupTv9b+7NmoHzvV6aaLEUREIirZvWlSOSlRqwTS0nah2Q730WyH+8jJW8GCVS8xf/Xb5OYvIc+zWbx+FL+s/Z4calI3bS8673AqezQ4lNQUbV4REZHKTH/JK5kaqY3IaHw5GY0vJ99zmLv6A6YuH8q67Hm457Nk40w+/O1e3l3wb+ql7kznxr3o2bQ3tVJ1MYKIiEhlo0StEkuxGrRrcCLtGpyIez7z1mQxdtkQflk/A3dYmbuQ4Yte5sPfXyM9pSGdGh3MkTv2oVHNHZIduoiIiJSAErUqwiyFjPo9yKjfA4Df1/3MyCVvM23VeDb6Rtbmr+LrpR/x5ZJPSbU6dKjfnWNb/Jk2dcvm/kMiIiJS9pSoVVE71WnHX3a5BoBVG5fz+aL/MfaPUazKW8OGvA1krRjNt8u+B2qSUWd3+uz8J7rt0EEXI4hItbVixQpeffVVLrnkkjJZ3m233Ua9evW45pprSjXfhAkT+O233+jdu3ep5ps7dy7ffvstZ555Zqnmk2iLxJMJzCzdzMaa2UQzm2pmt4fjbzOzBWY2IXwl3GvN7Kpwvilm9pqZpceVX2NmbmbV8gZkDWo2pm/rc7i782Ae6PQcfVueQeMaO+FAjucybfWP3DH9cfp8fSX9v7+HDxd8R05+4seTiIhUVStWrOCJJ55IWFaRz82cMGECH374Yanmyc3NZe7cubz66qvlFJUkSyQSNSAbOMLdOwNdgF5m1jMse9Ddu4SvrfZcM2sJXA5kuntHIBU4Paa8NXA0oCegA+lp6fRq8Wfu63wfz+/7DP3bXUhG7QyMNPLJZ97633hw1uv0+upaTvvmDl78+VNW56xPdtgiIlsbNgw6dID09OD/YcO2a3E33HADs2fPpkuXLlx77bWMGDGCww8/nDPPPJN99tkHgBNPPJHu3buz9957M3jw4E3zfvzxx3Tr1o3OnTtz5JFHbrXsZ555hmOPPZb167f8Ph06dCgdO3akc+fOHHLIIWzcuJFbb72VIUOG0KVLF4YMGcLYsWM54IAD6Nq1KwcccAAzZ84E4IUXXuDUU0/l+OOP3/SQ+K+//pouXbrw4IMPbldbSHRE4tCnB49HKLjBTI3wVZpHJqQBtc0sB6gD/BZT9iBwHfBeGYRapaSmpHJI854c0rwn7s6kFT8x5JfPmbB8Ntnksjh7Oc/P+YTn53xKvbTaHNa8M2fschgt61TLjkkRiZJhw+D006Eg8Zk+PRh+/XXo02ebFnnPPfcwZcoUJkyYAMCIESMYO3YsU6ZMoW3btgA899xz7LDDDqxfv559992Xk08+mfz8fC666CJGjhxJ27Zt+eOPP7ZY7mOPPcann37Ku+++u+nB6gUGDhzIJ598QsuWLVmxYgU1a9Zk4MCBZGVl8dhjjwGwatUqRo4cSVpaGp999hn/93//t+nh8aNHj2bSpEnssMMOjBgxgkGDBiV8aLxUXpFI1ADMLBUYB+wKPO7uY8zsWOAyMzsHyAKudvflsfO5+wIzG0TQY7Ye+NTdPw2X2QdY4O4Tizr3yswuBi4GaNOmep5cb2Z0brwbnRvvBsCva5fw33nDGbVoGmvzN7A6dz3DfhvDewvGUDOlJpmNd+WvbY9g74ZtdF6biFS8G27YnKQVWL8ebrxxmxO1RHr06LEpSQN45JFHeOeddwD49ddfmTVrFkuWLOGQQw7ZNN0OO2y+sv7ll1+mVatWvPvuu9SoUWOr5R944IGce+65/OUvf+Gkk05KGMPKlSvp168fs2bNwszIycnZVHb00UdvsT6peiKTqLl7HtDFzBoB75hZR+BJ4A6C3rU7gAeA82PnM7PGwAlAW2AFMNTMzgbeBm4CjinBugcDgyF41mfZ1Khya123GTd2OBM6wPLs1Qz55Ws+WZjFko2ryc7PYdTS6YxaOoMUUtijfmvOyjiYQ3bcm1SLytF0EanSfv458fjZs8t0NXXr1t30fsSIEXz22WeMHj2aOnXqcNhhh7FhwwbcvdAfrB07dmTChAnMnz9/i4SvwFNPPcWYMWP44IMP6NKly6bevFi33HILhx9+OO+88w5z587lsMMOSxifVE2R+6vq7iuAEUAvd1/k7nnung88A/RIMMtRwBx3X+LuOQQJ2gFAe4LkbaKZzQVaAePNbKfyr0XV0rhWffrv1pt3DrmVzw67iyt260Ob2s0AyPN8pq76hf+b9AoHDf8/Thv1AEPmfcOG3I1JjlpEqrR2hTznuH37bV5k/fr1Wb16daHlK1eupHHjxtSpU4cZM2bw3XffAbD//vvz1VdfMWfOHIAtDn127dqVp59+mj59+vDbb79ttczZs2ez3377MXDgQJo2bcqvv/66VRwrV66kZcuWQHBe2rbGL5VTJBI1M2sW9qRhZrUJkq8ZZtYiZrK+wJQEs/8C9DSzOhb8pDkSmO7uk929ubtnuHsGMB/o5u6/l2ddqrr0tJqcustBvHbQtXx91N38q8tf2adhG1IJfk3+sm4JD838H4d/cSvHjbiTx3/8gD82rEpy1CJS5dxzD9SOe+JK7drB+G3UpEkTDjzwQDp27Mi11167VXmvXr3Izc2lU6dO3HLLLfTsGVzz1qxZMwYPHsxJJ51E586dOe2007aY76CDDmLQoEH8+c9/ZunSpVuUXXvtteyzzz507NiRQw45hM6dO3P44Yczbdq0TRcTXHfdddx4440ceOCBRV592qlTJ9LS0ujcubMuJqhCLDiPP8lBmHUCXiS4YjMFeMPdB5rZywRXgTowF/ibuy80s52BZ929dzj/7cBpQC7wA3Chu2fHrWMuwZWhW35K4mRmZnpWVlYZ1q56cHemrfyVl+Z8yZils8j24PYehmMGdVNqcWjzjly46zHsWLtRcoMVkUiaPn06e+21V8lnGDYsOCdt9uygJ+2ee+D448svwCogURub2Th3z0xSSFKMSCRqUaJErWz8tm4ZL88dwZe/T2J13oYtyuqn1uaEVvtxdtvDqF9DzyAVkUCpEzUpNSVqlY8StThK1Mre6px1vDHvG95bMIZlG7c8f6JFeiPO2OVQ+rTsQY3UyFzbIiJJoESt/ClRq3yUqMVRola+Vm5cy/OzP+ODheNYl1dwdDo4PLpr3RZc0O5oDmq+t275IVINTZ8+nT333FOf/3Li7syYMUOJWiWjRC2OErWK8+vapQz+6WNGLZ1Kjm8+QTYFo2vj9vxt12Pp0LB1EiMUkYo0Z84c6tevT5MmTZSslTF3Z9myZaxevXqr24QoUYs2JWpxlKglx6Tlcxj808dMWjmX/PChFO5QKyWNQ5rvTf9de7NT7cZJjlJEylNOTg7z589nw4YNxU8spZaenk6rVq22uvGuErVoU6IWR4lacuV7Pl/+PpkX5nzGnLWLcKDgh3W91HSOb7kv57Q9kvo16iQ1ThGRqkKJWrQpUYujRC06cvJzefvXb3l93kiWbCy4F5uTYtCsVgNOa30oJ7Y+gJopughBRGRbKVGLNiVqcZSoRdPanA28MOcz3v9tDGvjbvfRru6O9Gt7NIc234cUPcJKRKRUlKhFmxK1OErUom/h+j945qeP+HrJFLI9N7ipLsFFCPs0bMuF7Y9ln8ZbP1NPRES2pkQt2pSoxVGiVrlMXzmPZ376iAkrZm+6CAEgjVQObNaBi3f9My3rNE1ihCIi0aZELdqUqMVRolY5uTujlkzl+Tmf8POahVuU1UmpyZE7daVfuz/RtFaDJEUoIhJNStSiTYlaHCVqlV9ufh7D5o9myC9fsih75RZl9dNq06vFvvw140ga1KybpAhFRKJDiVq0KVGLo0StalmXm83rc7/kfwtGszx37RZljWvUo0/Lnpy2y2HUSUtPUoQiIsmlRC3alKjFUaJWda3MWct/f/6cjxd+z+q89TElTvNaDenb6iBOan0QtVJrJi1GEZGKpkQt2pSoxVGiVj0s2bCCF38ezheLfmBdfvamm+oazs61m3Bq60PovfP+1NA92kSkilOiFm1K1OIoUat+FqxbyvM/f8w3SyaxkVwgSNjMoE3tHTmjzZEcvmM3UlNSkxypiEjZU6IWbUrU4ihRq97mrFnIC3M+ZOyyaeSSv2m8Ae3q7MzZbY/hwKad9MBoEakylKhFmxK1OErUpMCMlXN5fs7HTFzxI3kx92hLIYU96reh3y696NZkdyVtIlKpKVGLNiVqcZSoSSI//PEjL839hGmr5pAXfmYcSCWVTo3ac367Y+nQMCOpMYqIbAslatGmRC2OEjUpirvz7dIpvDJ3OD+umU/spyeNVDIb78n57Y+lff2dkxajiEhpKFGLNiVqcZSoSUnlez5fLBrP6798wZy1v+NO+NRRqGlpHNh0Hy5o30uPsBKRSFOiFm1K1OIoUZNtkZefx4cLxzLkly9ZsH4ZAMFHy6idUpMjd+rCuW3/RNP0hkmNU0QknhK1aFOiFkeJmmyvjfm5vDv/G96cN5LFGzc/wsqBumm16bXjvvRrdyQN9QgrEYkAJWrRpkQtjhI1KUvrc7MZ8ssI3ps/mmU5a3A3HHA3GqbV54RW+3F2xiHUraFHWIlIcihRizYlanGUqEl5WZ2zjpfmfMEHC7JYmbuOfLcgcXOjSc2GnNpmf87I2J/0ND3CSkQqjhK1aFOiFkeJmlSEZdmr+M/sL/hkwQRW52aTHzwLAXdomb4DZ7Y9kFPaZFIjVY+wEpHypUQt2pSoxVGiJhXtt3V/8MSPn/PF71PZkB88wirfATd2qduE09rsR99dulNbPW0iUg6UqEWbErU4StQkmX5evZjHf/yMkYtmsjE/eIRVwdWjTWrU47hWnTl314Nokq4LEUSkbChRi7ZIJGpmlg6MBGoBacCb7v5PM7sNuAhYEk76f+7+YYL5rwIuJLiwbjJwnrtvMLP7geOBjcDscPyKomJRoiZRMWX5fAbP+pJvF89mo+dvSthwSE+pxeEt9uDCXQ9i90Y7JjtUEanElKhFW1QSNQPquvsaM6sBjAKuAHoBa9x9UBHztgyn7+Du683sDeBDd3/BzI4BvnD3XDO7F8Ddry8qFiVqEkXz1/zBC7NH8fGCqazI2UC+BzfWdYcalka3xm04b/f9OWSnXfXsUREpFSVq0RaJM5U9yBbXhIM1wldpMsg0oLaZ5QB1gN/C5X4aM813wCnbH61IxWtVbwdu7tyHmzv3YU1ONq/8PIahc8azYP0KNubn8d2yuYwePRfDaF+3KX/dtQcn7tKJ9LQayQ5dRES2QyR61ADMLBUYB+wKPO7u14eHPs8FVgFZwNXuvjzBvFcAdwHrgU/d/awE0/wPGOLu/01QdjFwMUCbNm26z5s3r6yqJVKuNubl8tH8abw4ezTTVizadBFCwWHSHWvV56SMzpy/x340qlU7ydGKSBSpRy3aIpOoFTCzRsA7wN8Jzk1bStC7dgfQwt3Pj5u+MfAWcBqwAhhKcI7bf2OmuQnIBE7yYiqsQ59SmY1ZPJdnZ37L6MXz2JCXB2y+GKFBWjpH77w7f+vQk/YNmyU1ThGJDiVq0RaJQ5+x3H2FmY0AesWem2ZmzwDvJ5jlKGCOuy8Jp3sbOAD4bzjcDzgOOLK4JE2kstuveQb7Nc8AYO7qZTw9bTQfL5jBqpxsVuVs4M25k3lz7mRqWir7NWtD/w7703OnXXRem4hIREUiUTOzZkBOmKTVJki+7jWzFu6+MJysLzAlwey/AD3NrA7Boc8jCQ6TYma9gOuBQ919XXnXQyRKMuo34e79juNujmNl9nqemzmWN3+ezML1q9iYn8fI3+cwcuFcUjD2aNicC/balz7tOlAjJTXZoYuISCgShz7NrBPwIpAKpABvuPtAM3sZ6EJw6HMu8Dd3X2hmOwPPunvvcP7bCQ595gI/ABe6e7aZ/URwy49l4aq+c/f+RcWiQ59S1WXn5fLOnCm8OCOLH1cuIS8f8M09ai1rN+Avu3Xi3L0yaVhLzyAVqep06DPaIpGoRYkSNalO3J1Rv83h6WljGfv7r2zMz9uUtHk+NKxRm94Ze9C/Uw8yGjZOcrQiUh6UqEWbErU4StSkOvtpxVIenzSa4fN+YvXGjRTcYBc30lNSOaDFLlzStQeZO7XSeW0iVYQStWhTohZHiZpIYOn6tTwz+XvemTWNxevWhgkb4EYKsE+TnTi/U3eO220PUlNSkhytiGwrJWrRpkQtjhI1ka1tyM3hjRmTeXnqBH5a/geez6akzRzaNGjI6R06cXanztSvVSvZ4YpIKShRizYlanGUqIkULd+dL+fOZvAP4/hh4QJy8sPvkHwwjB1q1ebPu+3OxfvuS8sGDZIbrIgUS4latClRi6NETaR0pi1ZzOPff8dXc+awLie4yW7BYdLaKakcktGW/j160HnnnZIap4gkpkQt2pSoxVGiJrLtfl+zhv+My+LdadNZtm49BkFPWz6kmLF3k2ac0b0zffbek/Qaeg6pSBQoUYs2JWpxlKiJlI0NOTm8Nmkyr/4wkTl/LIe84Hw2C3vbmtauzZ/22p3ze3aj9Q6Nkh2uSLWlRC3alKjFUaImUvby8vMZ9fM8nhs7jnHzFrAxNx/YnLTVNKNr65b8db8uHLnnrqSk6NYfIhVFiVq0KVGLo0RNpPz9tmIVL4wZx0dTZ7F49dpNCVtwQQK0bFif4zvvydn7d6VJvbpJjlakalOiFm1K1OIoUROpWBtz8/hoykxeGTOBab8tIi+8HqEgeaudlsYBu7am30GZdG/bUjfaFSljStSiTYlaHCVqIsn14+9LeH7UOL6YNpvVGzaCs+mihBSDLq1b8HC/PuxQr06yQxWpEpSoRZsStThK1ESiY112DkO/n8SbYyYzZ+lyPOxtq5OWxit/P53dWzRLboAiVYAStWhTohZHiZpINLk73/w4j6tffp+12TmkmvHIuX04tEO7ZIcmUqkpUYs2PaBPRCoFM+OgPTL4/KaLaNW4AXn5zmX/eY9XR01IdmgiIuVGiZqIVCr1atfi/evPI7PtzgDc/faX3PXWF0mOSkSkfChRE5FKJy01hecv+Qsn77c3AK9/M5ELHh9Kfn5+kiMTESlbStREpFIyM277yzFcf8KhGPD9T/P5853Ps3p9drJDExEpM0rURKRSO/vQbjxx4YmkphgL/lhFr9ue5ZclK5IdlohImVCiJiKV3kEd2jL0mrOpXSONNRs20vdfL5I169dkhyUist2UqIlIlbBri6Z8dOsFNK1Xh9zcfC565E0++2FWssMSEdkuStREpMrYoX4dPrztAvbcuSmeD9c+8z5vjpyU7LBERLaZEjURqVJq1UjjtevPYv892gBw16uf88wHY5IclYjItlGiJiJVTkpKCk9efjJ9euyFOTz53rfc86rutSYilY8SNRGpsm4/txfn/2lfzGHolxO59rFhyQ5JRKRUlKiJSJV22UkHcfVfDgV3vvzhJ86/61Xy8nRjXBGpHJSoiUiVd+bR3bjzwl6kAJN/+p3TbnqRDdk5yQ5LRKRYStREpFo4tmcHHrqiLylmzFu4nJOu+Q8r16xPdlgiIkWKRKJmZulmNtbMJprZVDO7PRx/m5ktMLMJ4at3IfNfFc43xcxeM7P0cPwOZjbczGaF/zeuyHqJSLQc0Kktz990BjVTUli2Yh0nX/Usvy9ZmeywREQKFYlEDcgGjnD3zkAXoJeZ9QzLHnT3LuHrw/gZzawlcDmQ6e4dgVTg9LD4BuBzd98N+DwcFpFqrEP7nXjtX+dQu2Yaa9Zu5PRrnueneUuSHZaISEKRSNQ8sCYcrBG+vBSLSANqm1kaUAf4LRx/AvBi+P5F4MTtj1ZEKrs2LXbgzQfOp1HddLKz8zjvxv/yw9T5yQ5LRGQrkUjUAMws1cwmAIuB4e5ecIfKy8xskpk9l+jQpbsvAAYBvwALgZXu/mlYvKO7LwynWwg0L+96iEjl0LRRPd565CJaNKlPXm4+f799CF9992OywxIR2UJkEjV3z3P3LkAroIeZdQSeBNoTHA5dCDwQP1+YvJ0AtAV2Buqa2dmlWbeZXWxmWWaWtWSJDoGIVBd1a9dkyMMXsGubpni+83+DhjHs04nJDktEZJPIJGoF3H0FMALo5e6LwgQuH3gG6JFglqOAOe6+xN1zgLeBA8KyRWbWAiD8f3Eh6xzs7pnuntmsWbOyrZCIRFqNtFReuO8cMvdpA/lw31PDeWno6GSHJSICRCRRM7NmZtYofF+bIPmaUZBkhfoCUxLM/gvQ08zqmJkBRwLTw7JhQL/wfT/gvXIIX0QquZQU46FbT+VPB+0JDs+8MoqHn/k82WGJiEQjUQNaAF+a2STge4Jz1N4H7jOzyeH4w4GrAMxsZzP7ECA8l+1NYDwwmaBOg8Pl3gMcbWazgKPDYRGRrZgZt/7jOM7o0w0c3np/PLfc/S7upbmuSUSkbJm+hLaUmZnpWVlZyQ5DRJLo9XfG8sTzX4HDPnvuzCP3nEFqalR+14qULTMb5+6ZyY5DEtM3j4hInNP79uCmK3tj7kyZtoDzBjxHth45JSJJoERNRCSBPx2xN/fedgopGL/8+gdnnvc0q1bpkVMiUrGUqImIFKJnZjse//dZpKWmsGzZWs7q9zSLF+uRUyJScZSoiYgUocMeO/P8k+dTu1YN1qzZwDnnDmbOHN1vUUQqhhI1EZFitG69A/99/mIa1E1n44Zc/nbRc0ya9EuywxKRakCJmohICTRpUo9XXxlA82b1ycvN4+or/su3o2YmOywRqeKUqImIlFDdurV4+b/9aZvRjPw859b/e5OP/vdDssMSkSpMiZqISCnUqJHG4OcupEuX4JFT/77nA1576ZtkhyUiVZQSNRGRUkpJMQY9fDaHHbYnAM899SWPP/hJkqMSkapIiZqIyDYwM26+82ROPCW4ofu7Q7/njpve1COnRKRMKVETEdkOl17Vi/P7Hw7ufP3FdK4Z8CJ5efnJDktEqgglaiIi2+mMcw7k2lv6ADB5wi/87ayn2LgxN8lRiUhVoERNRKQMHNO7M3c+cDpmxi9zl9Kv76OsXr0h2WGJSCWnRE1EpIz0OGA3HnrmPFJTU1i2ZDX9TnyYpXrklIhsByVqIiJlaK+OrRj8Sn9qpaexZtUGzj/pceb9rEdOici2KTZRM7N0MzvFzB42s6Fm9pKZXWdme1dEgCIilU3rjKY8/9Zl1G+QTnZ2Dpec9RTTJv+a7LBEpBIqMlEzs9uAb4D9gTHA08AbQC5wj5kNN7NO5R2kiEhl07RZA1589wqaNq9Pbk4e/zj/OcaO+jHZYYlIJWNF3fPHzP7s7h8UUd4caOPuWeURXDJkZmZ6VlaVqY6IJFl2dg6Xnf00v8xZCjjXDTyJI3t3TnZYIpuY2Th3z0x2HJJYkT1qRSVpYfniqpSkiYiUtVq1avDU65ew1z6twOG+W97m7Ve+TXZYIlJJFHfocyczuzt87VxRQYmIVCWpqSk8+NwF9Dh4NwCeHvQJzz48PMlRiUhlUNzFBP8NXy+HLxER2QZmxh0Pnc0xx3cFYOgLoxh069tJjkpEoq64RG0t0CZ8ZZd/OCIiVdvVt53IX847GIDhwyZw82Uv6/mgIlKo4hK1MwmStLbAGeUfjohI1XfB5UfR/5peAGSNmsUVZw/W80FFJKHiLiZY6+5Pu/uT7q7ba4uIlJG+Z+3PDXefAjg/TpnPxX0fITs7J9lhiUjEFHcxwc1m1riI8iPM7LiyD0tEpOo7/NhO/OupczEzFsxdynm9/82a1euTHZaIREhxhz4nA++b2edmdn/4RIJbzexlM5sMHE9wI1wREdkG3Xq259HX+pOWlsIfi1fR70+DWLp4VbLDEpGIKO7Q53vufiDQH5gKpAKrCK4E7eHuV7m7HmInIrIddt1rZ5557wpqpaexdtUGLuj1AL/O0VeriBTzZILqSE8mEJFkWb50NRcd/xBrVq4nLS2FB175G3vs0zrZYUkVpycTRFuxD2WvCOGD38ea2UQzm2pmt4fjbzOzBWY2IXz1TjDvHjHlE8xslZldGZZ1MbPvwvFZZtajgqsmIlJijZvW58Xh19K0eQPycvK46rQnyPp6ZrLDEpEkikSiRnCPtiPcvTPQBehlZj3DsgfdvUv4+jB+RnefWVAOdAfWAe+ExfcBt4dlt4bDIiKRVbdeOs99cjUtM5ri+c4tFz7PqE8mJzssEUmSSCRqHlgTDtYIX9tyTPZIYLa7zytYNNAgfN8Q+G27AhURqQA1a9Xg6f9dye77tAJ37rrsZSaO/inZYYlIEpQoUTOz3cMrP6eEw53M7OayDMTMUs1sArAYGO7uBVeTXmZmk8zsuaJuFRI6HXgtZvhK4H4z+xUYBNxYyLovDg+NZi1ZohN4RST5UtNSeeiNS9mzc3CO2g3nDObnGfqtKVLdlLRH7RmCJCcHwN0nESRFZcbd88JDlK2AHmbWEXgSaE9wOHQh8EBh85tZTaAPMDRm9ADgKndvDVwF/KeQdQ9290x3z2zWrFkZ1EZEZPuZGYNeG0Crtk3BncuOf5Dff12W7LBEpAKVNFGr4+5j48bllnUwAO6+AhgB9HL3RWECl0+QLBZ1McCxwHh3XxQzrh9Q8NTjocXMLyISOalpqTw27Cp2aFYfz82n/5/uY/kS3WdNpLooaaK21MzaE543ZmanEPRwlQkza2ZmjcL3tYGjgBlm1iJmsr7AlCIWcwZbHvaE4Jy0Q8P3RwCzyiRgEZEKVCu9BoM/vY56DdLJXreRvx1zL2vXbEh2WCJSAUqaqF0KPA3saWYLCM79GlCGcbQAvjSzScD3BOeovQ/cZ2aTw/GHExy+xMx2NrNNV4CaWR3gaDb3nhW4CHjAzCYC/wIuLsOYRUQqTN366Tw9/Hpqpqex+o+1XHLMPeRs1LNBRaq6Ut3w1szqAinuvrr8Qkou3fBWRKJs/s+L6X/UPeTl5NF2rxY89vF1pKRE4gJ+qaR0w9toKzJRM7N/FDWzu/+7zCNKMiVqIhJ1M36Ywz9OfAjPczofuBv3DPl7skOSSkyJWrQV9zOsfvjKJDjU2TJ89Qc6lG9oIiKSyJ5d23Lbfy4Cg4mjfuTuAc8nOyQRKSfFPZT9dne/HWgKdHP3q939aoInALSqiABFRGRrPY7qyNX/PguAke9m8cRNQ5IckYiUh5Ke2NAG2BgzvBHIKPNoRESkxI46dT8uurUvAP/7z1e8+u8PkhyRiJS1tBJO9zIw1szeIbhFR1/gpXKLSkRESuSkvx3B8iUrefOx4bx87/s0atqA3uccnOywRKSMlKhHzd3vAs4DlgMrgPPc/V/lGJeIiJTQBTf35YiTg/t5P3rNK3z70YTkBiQiZaakz/psAywF3glfy8JxIiISAdc81o9uB+8JDnec8yRTx+gh7iJVQUnPUfsAeD98fQ78DHxUXkGJiEjpmBl3vPF32ndshbtz9Z/v45eZeoi7SGVX0kOf+7h7p/C1G8EzM0eVb2giIlIaKSkpPPTpDbRo0xTy8rnkkNtZPF8PcRepzLbpdtbuPh7Yt4xjERGR7ZRWI40nR/2TRs3qk5udy9963sLKZWuSHZaIbKOSnqP2j5jXNWb2KrCknGMTEZFtkF6nFs+MuZPa9WqxftV6+u93E+vX6iHuIpVRSXvU6se8ahGcs3ZCeQUlIiLbp37jugz+7g5qpKfxx6KVDDjgFnJycpMdloiUUkkTtWkFTylw97vc/RXg+PIMTEREtk+zVk145MtbSU01Fs5ezD+OuIP8/PxkhyUipVDSRO3GEo4TEZEIabt3a+7+3/WYwY/j53DLyQ8kOyQRKYUiEzUzO9bMHgVamtkjMa8XAPWhi4hUAp0P3pObX/k7AFkfT+L+C55MckQiUlLF9aj9BmQBG4BxMa9hwJ/KNzQRESkrB/XJ5IpHzwPgs1dGMfj6/yY5IhEpiSKf9enuE4GJZvaKu6sHTUSkEut9/uGsWraa528ZwlsPf0TDZg047Zo+yQ5LRIpQ3KHPN8K3P5jZpPhXBcQnIiJl6PRr+9DnkmMA+M9NrzP85ZFJjkhEilJkjxpwRfj/ceUdiIiIVIxL/92PZQuWM+qdsdx/wVM0aFaf/Xp1TXZYIpJAkT1q7r4wfHuJu8+LfQGXlH94IiJSHm55/Qo6HrRn8L7P/czM0kPcRaKopLfnODrBuGPLMhAREak4Zsb9w28mo0NL8HwuP/AWfv1RD3EXiZrizlEbYGaTgT3izk+bA+gcNRGRSiw1NYXHxvyL5q2b4vnOgO7Xs3TBH8kOS0RiFNej9irBEwiGhf8XvLq7+9nlHJuIiJSzmrVq8MykQTTYoR7Z6zZyUedrWLVcD3EXiYrizlFb6e5z3f2M8Ly09YAD9cysTYVEKCIi5ap2vXSemfJv0uvVYs3ytVzc+Rqy12cnOywRoYTnqJnZ8WY2C5gDfAXMBT4qx7hERKQCNW7WgKfG30eNWmksW7CcSzJvIC8vL9lhiVR7Jb2Y4E6gJ/Cju7cFjgS+KbeoRESkwrVsvxMPjhxISloKv0yfz9WH/RN3T3ZYItVaSRO1HHdfBqSYWYq7fwl0Kb+wREQkGfbIbM9d/7seSzGmfjOTf/a9P9khiVRrJU3UVphZPWAk8IqZPUwZPpTdzNLNbKyZTTSzqWZ2ezj+NjNbYGYTwlfvBPPuEVM+wcxWmdmVMeV/N7OZ4XLvK6uYRUSqqsxjunDjK8H9zkcP+54HLtRD3EWSpbgnExQ4geDB7FcBZwENgYFlGEc2cIS7rzGzGsAoMys4B+5Bdx9U2IzuPpOwd8/MUoEFwDvh8OFh7J3cPdvMmpdhzCIiVdbhpx3IqmWreeyy//Dxc1/QqHkDLvjXWckOS6TaKVGi5u5rYwZfLOsgPDgJouB68Brha1tOjDgSmB1eoQowALjH3bPD9Sze3lhFRKqLEy7pxfLfV/DKnW/x+j3v0minhpx8uZ4oKFKRirvh7erwUGL8a7WZrSrLQMws1cwmAIuB4e4+Jiy6LLzJ7nNm1riYxZwOvBYzvDtwsJmNMbOvzGzfsoxZRKSqO3fg6fQ6/wgAnrryRb56Q9eRiVSk4u6jVt/dGyR41Xf3BmUZiLvnuXsXoBXQw8w6Ak8C7QkObS4EHihsfjOrCfQBhsaMTgMaE1yxei3whplZgnkvNrMsM8tasmRJ2VRIRKSKuPrZAfTo3Q2AO09/iGmjf0xyRCLVR0kvJqgw7r4CGAH0cvdFYQKXDzwD9Chi1mOB8e6+KGbcfOBtD4wF8oGmCdY52N0z3T2zWbNmZVUVEZEq445h17NbZjsArjzkZub/tDDJEYlUD5FI1MysmZk1Ct/XBo4CZphZi5jJ+gJTiljMGWx52BPgXeCIcLm7AzWBpWUTtYhI9ZGSksLDo+5kx12a4nnOJd2vY8XSMj0DRkQSiESiBrQAvjSzScD3BOeovQ/cZ2aTw/GHE1x1ipntbGYfFsxsZnWAo4G345b7HNDOzKYArwP9XHdvFBHZJjVq1uDpiQ9Qr1Fd1q/ewICu15K9YWOywxKp0kx5y5YyMzM9Kysr2WGIiETW4l+XcO7ul5OTncuuXdvy+Pf3kJISld/9UlpmNs7dM5MdhySmT5aIiJRK89bNeHDkHViq8dMPc7jt5EKv8xKR7aRETURESm2PfXfl1jeuBkth9LAsHr/iuWSHJFIlKVETEZFtclDf/bj0oXMBePexT3jr4Q+SG5BIFaRETUREttmJfz+Wk/8RPK3g6WteZtR73yc5IpGqRYmaiIhsl/73/5UDTgjORR94yr/56Yc5SY5IpOpQoiYiItvtn29eTbsuGWDGZQfczJJflyU7JJEqQYmaiIhst5SUFB4ZdQdNd25Mfp5zYZdrWfXH6mSHJVLpKVETEZEyUbNWDZ6ZeD91G9Vhw9psLupyPdnrs5MdlkilpkRNRETKTN2GdXh24v2k1UxjxeJVDNjvJnJzcpMdlkilpURNRETKVJMWjXn8u7uw1BQWzFrElUcMRE/BEdk2StRERKTMZXRoxb0f3QgGs8bN4cbj70t2SCKVkhI1EREpF50P2Yt/vnEVmPHDl1O565zHkh2SSKWjRE1ERMrNAcd145pnLgJg5FtjefjyF5IbkEglo0RNRETK1dFnHswlg/4KZnz43Jc8e8uQZIckUmkoURMRkXJ3woCj6XfLSQAMffBDhvz7/SRHJFI5KFETEZEKceb1J3DSZb2wlBSev+0t3n/uy2SHJBJ5StRERKTC/O2eMzji9P0BeOzKlxjx1tgkRyQSbUrURESkQl03+GIyj+4EZtxz/lNkfT452SGJRJYSNRERqXB3vHkle/VoD8DNJz/EjKzZSY5IJJqUqImISIUzMwZ9fCO7dGwNZlx5zN3Mnb4g2WGJRI4SNRERSYrU1BQe+/IWdtqlKThceujtLJyzJNlhiUSKEjUREUmaGjXTeOrb29lhx4bk5eTT/+DbWLpwebLDEokMJWoiIpJU6XVqMXjMHdRrXIfsdRvpf8BtrPpjTbLDEokEJWoiIpJ09RrW4Znv7iC9Xi3WrFzHxQf+k3VrNiQ7LJGkU6ImIiKR0Lh5Q54c+U9qptdgxeJVDDh0IBs35iY7LJGkUqImIiKR0SKjGQ9/9n+k1Uxj0bylXH7EneTl5SU7LJGkUaImIiKR0rZDK+5992pSUo25M37j2uMH4e7JDkskKZSoiYhI5Oy9367c/urfMTOmjZ3NLac/qmRNqqVIJGpmlm5mY81soplNNbPbw/G3mdkCM5sQvnonmHePmPIJZrbKzK6Mm+YaM3Mza1pBVRIRke2075Ed+b9nL8YshXFfTOPui55JdkgiFS4SiRqQDRzh7p2BLkAvM+sZlj3o7l3C14fxM7r7zIJyoDuwDninoNzMWgNHA7+Ucx1ERKSMHXxCd6586K8AfP3eeB65+r9JjkikYkUiUfNAwU1zaoSvbenjPhKY7e7zYsY9CFy3jcsTEZEk+9NZB3LBbScB8NFLo3hu4NtJjkik4kQiUQMws1QzmwAsBoa7+5iw6DIzm2Rmz5lZ42IWczrwWswy+wAL3H1iMeu+2MyyzCxryRI9vkREJGpOufQYTvn70QAMffRThj72aZIjEqkYkUnU3D0vPHzZCuhhZh2BJ4H2BIdDFwIPFDa/mdUE+gBDw+E6wE3ArSVY92B3z3T3zGbNmm1nTUREpDxccOvJ/OmsAwB47va3+eSVb5IckUj5i0yiVsDdVwAjgF7uvihM4PKBZ4AeRcx6LDDe3ReFw+2BtsBEM5tLkACON7Odyit2EREpX1c+dA4H9O4CwENXvsy3H05Iajwi5S0SiZqZNTOzRuH72sBRwAwzaxEzWV9gShGLOYOYw57uPtndm7t7hrtnAPOBbu7+e1nHLyIiFefmF/5Gp4P3AIM7znuaH0bOSHZIIuUmEoka0AL40swmAd8TnKP2PnCfmU0Oxx8OXAVgZjub2aYrQMPDnEcDOsNURKSKMzPufvMKduuaAcBNf3mY6ePmJDcokXJiuoHgljIzMz0rKyvZYYiISDHycvO45PA7+eXH30mtkcpjw28kY6+WyQ6r0jGzce6emew4JLGo9KiJiIiUSmpaKo9+9n80b70DeTl5XHHsfSycqyv3pWpRoiYiIpVWzVo1eHLELTRq1oCNG3K47Kh/sXzxqmSHJVJmlKiJiEilVqdeOk99dQt1G9Zm3ZpsLjnqbjasy052WCJlQomaiIhUeg2b1OPxz2+iVp2arFiyisuOvoe83LxkhyWy3ZSoiYhIlbBjqx14YNjVpKalsODnxVx9woPogjmp7JSoiYhIldG+Y2vuePVSMJg5fg4Dzxuc7JBEtosSNRERqVK6HrwnNzxxPqSk8N3wKTx6w5BkhySyzZSoiYhIlXPoid254OYTAPjw5VG88fjwJEcksm2UqImISJV0yoCjOP68gwF4/u7/8eW7upm5VD5K1EREpMq65M6/sH+vTgDcd9lLTBg1M8kRiZSOEjUREanSbnn2QvbslgHATWc/yeyp85MbkEgpKFETEZEqzcwY9PYVtGzfnPw85x8nPsSiX/9IdlgiJaJETUREqrzUtFQe++haGjWrHzxq6tj7WL1iXbLDEimWEjUREakW0uvU4slPr6dOvXTWrFzPpb3uJ2djbrLDEimSEjUREak2GjWtz8MfXE1arRos+W05V/R5kPz8/GSHJVIoJWoiIlKttGrXnPuGXEZKijFn2gJu+utTetSURJYSNRERqXb26p7BLYMvgBRjwjezGHTVK8kOSSQhJWoiIlIt9TymI5fffSoAX7w7jufvez/JEYlsTYmaiIhUW8eecQBnXHY0AG888TkfvPJNkiMS2ZISNRERqdbOuaY3R5yUCWY8dvNbfPf5lGSHJLKJEjUREan2rnngTDrvvysAt1/4HNN/mJfkiEQCStRERKTaMzPuerk/GXu2AIxrT3+c+XMWJzssESVqIiIiAKmpKTz87lU0a9mYvJx8Lu/zEH8sWZXssKSaU6ImIiISqlkrjSc+/Af1GtZm/bqNXHrcg6xfm53ssKQaU6ImIiISo16DOjzxwdXUqlOTFcvW8PcTHiI3Ny/ZYUk1pURNREQkTrOdG/HAkEtJTUthwdylXH+Wnl4gyaFETUREJIH2HVpyx3MXYCnGtPHzuOPSl5IdklRDkUjUzCzdzMaa2UQzm2pmt4fjbzOzBWY2IXz1TjDvHjHlE8xslZldGZbdb2YzzGySmb1jZo0qtmYiIlKZdT1gd6574AwARg+fyuO3v5PkiKS6iUSiBmQDR7h7Z6AL0MvMeoZlD7p7l/D1YfyM7j6zoBzoDqwDCj5Jw4GO7t4J+BG4sZzrISIiVcxhx3fl4puOA+D9/45myOAvkxyRVCeRSNQ8sCYcrBG+tuVkgCOB2e4+L1zup+6eG5Z9B7Ta7mBFRKTa6XvuIfTpdxCY8cKgj/nyfz8kOySpJiKRqAGYWaqZTQAWA8PdfUxYdFl46PI5M2tczGJOB14rpOx84KOyiVZERKqbATf3Yf+j9gYz7rtmCBPHzk52SFINRCZRc/e88PBlK6CHmXUEngTaExwOXQg8UNj8ZlYT6AMMTVB2E5ALvFLIvBebWZaZZS1ZsmQ7ayIiIlXVLY//lT27tgGDG/o9y5yZvyc7JKniIpOoFXD3FcAIoJe7LwoTuHzgGaBHEbMeC4x390WxI82sH3AccJYXcm21uw9290x3z2zWrFlZVENERKogM+P+l//GzhlNAbjitMdZvHBFcoOSKi0SiZqZNSu4ItPMagNHATPMrEXMZH2BKUUs5gziDnuaWS/geqCPu68r06BFRKRaSquRyuPvXE6jJvXIyc7l0r6Psnql/sRI+YhEoga0AL40s0nA9wTnqL0P3Gdmk8PxhwNXAZjZzma26QpQM6sDHA28Hbfcx4D6wPDw1h1PVUBdRESkikuvXZOnhl1Jer1arF61nktOfpSN2TnJDkuqINOdlreUmZnpWVlZyQ5DREQqgQXzltL/xEfIycmj3V478fgbl2FmyQ6rVMxsnLtnJjsOSSwqPWoiIiKVTstdmnLvCxeQkmr8PON3bhrwoh41JWVKiZqIiMh26NB5F2556CwwY/y3P/HALfFn4YhsOyVqIiIi22n/w/fiiltPAOCz/03ghUeHJzkiqSqUqImIiJSBY0/Zl7MGHAHA68+O5P03xhQzh0jxlKiJiIiUkb8OOIKj+3QF4NF/fcCITyYnOSKp7JSoiYiIlKF/DOxL1wN2BeDuG95kzKiZSY5IKjMlaiIiImXIzLjr8b+yV5fWYHDrFa8yadzcZIcllZQSNRERkTKWkpLCA/+5gHZ77AjAtRc/z8yp85MclVRGStRERETKQWpqCo++9Ddah88FvfK8Z5n9ox7iLqWjRE1ERKScpNVI5YnXBrBjy8bk58Pl/Qbz69wlyQ5LKhElaiIiIuWoZs00nnnjUnZoWo/cnHwuPftpFi74I9lhSSWhRE1ERKSc1UqvwbNvXkbDxrXJ3pDLgDOeYuniVckOSyoBJWoiIiIVoG69dAa/cRn16tVi/bqN9D/9CVatWJfssCTilKiJiIhUkEaN6/LUkEtIr1OT1as2cNGpj7N2zYZkhyURpkRNRESkAjXbsSFPvtKfmrXSWLF8LRf/5Qk2bNiY7LAkopSoiYiIVLCdW+/Aoy9fTI0aqSxdvIr+ZzxJTk5ussOSCFKiJiIikgQZ7Zrz4HMXkJqawsJfl3PpWU+Tl5uX7LAkYpSoiYiIJMlue+3MfYPPxcyY9/MSrjjvWfLz85MdlkSIEjUREZEk6ti5Df967GwAZk1fyDUXP4+7JzkqiQolaiIiIknWbb/23Pbv0wGYOuFX/u/vLytZE0CJmoiISCTsf8ieXH/HSQCM/+5n7rhuSJIjkihQoiYiIhIRRxzbictvPA6Ab76cwf3/fDvJEUmyKVETERGJkD+fnMmFlx8NwGcfTOLx+z9MckSSTErUREREIubUcw7kzAsOBmDYkLG8+OQXSY5IkkWJmoiISAT1G3Akfc/sCcCr/xnJ0Je/SXJEkgxK1ERERCKq/z96cWzfbgBMuf0xVrdqC+np0KEDDBuW5OikIihRExERibArb+rDRe3Wc+PPb1J/wVzIzobp0+H005WsVQNK1ERERCLulCn/I93jngW6fj3ceGNyApIKE4lEzczSzWysmU00s6lmdns4/jYzW2BmE8JX7wTz7hFTPsHMVpnZlWHZDmY23Mxmhf83ruCqiYiIbL+ff048fvbsio1DKlwkEjUgGzjC3TsDXYBeZtYzLHvQ3buEr62uUXb3mQXlQHdgHfBOWHwD8Lm77wZ8Hg6LiIhULu3aJR7fvn3FxiEVLhKJmgfWhIM1wte2PDvjSGC2u88Lh08AXgzfvwicuD1xioiIJMU990Dt2luOq107GC9VWiQSNQAzSzWzCcBiYLi7jwmLLjOzSWb2XAkOXZ4OvBYzvKO7LwQI/29eyLovNrMsM8tasmTJ9lVERESkrPXpA6+/HlztWatW8P+QIXD88cmOTMqZRe2hr2bWiODQ5d+BJcBSgt61O4AW7n5+IfPVBH4D9nb3ReG4Fe7eKGaa5e5eZLKXmZnpWVlZZVATERGR6DOzce6emew4JLHI9KgVcPcVwAigl7svcvc8d88HngF6FDHrscD4giQttMjMWgCE/y8un6hFREREyl4kEjUzaxb2pGFmtYGjgBkFSVaoLzCliMWcwZaHPQGGAf3C9/2A98okYBEREZEKkJbsAEItgBfNLJUgeXzD3d83s5fNrAvBoc+5wN8AzGxn4Fl37x0O1wGOLiiPcQ/whpldAPwCnFoBdREREREpE5FI1Nx9EtA1wfi/FjL9b0DvmOF1QJME0y0juBJUREREpNKJxKFPEREREdmaEjURERGRiFKiJiIiIhJRStREREREIipyN7xNNjNbAswrdsLSa0pw896qpCrWCVSvyqQq1glUr8qkKtRpF3dvluwgJDElahXEzLKq2p2fq2KdQPWqTKpinUD1qkyqYp0kWnToU0RERCSilKiJiIiIRJQStYozONkBlIOqWCdQvSqTqlgnUL0qk6pYJ4kQnaMmIiIiElHqURMRERGJKCVqIiIiIhGlRK2UzOxUM5tqZvlmlhlXdqOZ/WRmM83sTzHja5rZYDP70cxmmNnJCZabYWbrzWxC+HoqpmxEuMyCsuZVpF7dzWxyuOxHzMwqQ51ipm1jZmvM7JqYcZV2WxVTr3LdVuVZLzPrEbM9JppZ35iyct1eSapTZd5WR5vZuDD+cWZ2RExZZd1WRdWp3LeVVAHurlcpXsBewB7ACCAzZnwHYCJQC2gLzAZSw7LbgTvD9ylA0wTLzQCmFLLOLdZVheo1FtgfMOAj4NjKUKeY5bwFDAWuqQrbqph6leu2Kud9sA6QFr5vASyOGS7X7ZWkOlXmbdUV2Dl83xFYUFGfrSTVqdy3lV6V/5WGlIq7TwdI8MPnBOB1d88G5pjZT0APYDRwPrBnOH8+EbyLdUXXy8xaAA3cfXQ4/BJwIsGXVZkozzqZ2YnAz8Dasoq3pCq6XhWxrcK4yqVe7r4uZjAdqLArqCq6TlVgW/0QMzgVSDezWuHyylVF1wnYgQrYVlL56dBn2WkJ/BozPB9oaWaNwuE7zGy8mQ01sx0LWUZbM/vBzL4ys4Pjyp4Pu/tvqeDu8fKqV8twWVsstywDL8J21cnM6gLXE/yaTqRSbqsi6pXMbVWw/u3aB81sPzObCkwG+rt7bkxxMrZXedWp0m+rGCcDP8QlaZVyW8WIrVOyt5VUEkrUEjCzz8xsSoLXCUXNlmCcA2lAK+Abd+9G8CtsUIJpFwJt3L0r8A/gVTNrEJad5e77AAeHr79WgXoVttxSSVKdbgcedPc1Ccoq87YqrF5lsq0gafXC3ce4+97AvsCNZpYeFm339opYnSr9tgrXvTdwL/C3mNGVdluF646vU5ltK6nadOgzAXc/ahtmmw+0jhluBfwGLAPWAe+E44cCFyRYZzaQHb4fZ2azgd2BLHdfEI5fbWavEnS7v1TaACNWr/nhsuKXWyrJqBOwH3CKmd0HNALyzWyDuz9WmbcVhdSL4Jy17d5WkLR6xa5/upmtJThXqEw+WxGrU5l8rsLlJqVeZtYqnO4cd58dE0+l3VaF1KnMtpVUbepRKzvDgNPNrJaZtQV2A8a6uwP/Aw4LpzsSmBY/s5k1M7PU8H27cP6fzSzNzJqG42sAxwFTyrsyMcqlXu6+EFhtZj3DQxjnAO+Ve20C21Undz/Y3TPcPQN4CPiXuz9W2bdVYfVK8raC7d8H25pZWvh+F4ITxucmeXuVS52qwLZqBHwA3Oju38SMr8zbqhEJ6hSBbSWVRWmvPqjuL6AvwS+hbGAR8ElM2U0EVwTNJObqHWAXYCQwCfic4FAgQB9gYPj+ZIITTScC44Hjw/F1gXHhvFOBhwmvOKrM9QrLMgm+bGcDjxE+KSPqdYpbx22EV0dW9m1VWL0qYluV8z7413B7TAj3wRMrantVdJ2qwLa6meBClgkxr+aVfFslrFNFbSu9Kv9Lj5ASERERiSgd+hQRERGJKCVqIiIiIhGlRE1EREQkopSoiYiIiESUEjURERGRiFKiJiIiIhJRStREREREIkqJmoiIiEhEKVETERERiSglaiIiIiIRpURNREREJKKUqImIiIhElBI1ERERkYhSoiYiIiISUUrURERERCJKiZqIiIhIRClRE5EtmNmaclhmHzO7IXx/opl12IZljDCzzFLO85CZHRK+39vMRpvZi2aWEo4bZGZHlDYWEZGKokRNRMqduw9z93vCwROBUidqpWVmOwA93X1kOOofQB8gCzgmHPcocEN5xyIisq2UqIlIQha438ymmNlkMzstHH9Y2Lv1ppnNMLNXzMzCst7huFFm9oiZvR+OP9fMHjOzAwiSpfvNbIKZtY/tKTOzpmY2N3xf28xeN7NJZjYEqB0T2zFh79h4MxtqZvUSVOEU4OOY4VTAgXzAANx9HtDEzHYq08YTESkjStREpDAnAV2AzsBRBMlVi7CsK3AlQc9YO+BAM0sHngaOdfeDgGbxC3T3b4FhwLXu3sXdZxex/gHAOnfvBNwFdIcgmQNuBo5y924EPWT/SDD/gcC4mOGHgQ+A/YFPY8aPD6cVEYmctGQHICKRdRDwmrvnAYvM7CtgX2AVMNbd5wOY2QQgA1gD/Ozuc8L5XwMu3o71HwI8AuDuk8xsUji+J0GC+E3YkVcTGJ1g/hbAkoIBd/8B2C/BdIuBnbcjThGRcqNETUQKY0WUZce8zyP4Lilq+qLksrl3Pz2uzAuJa7i7n1HMctcnWF4i6eG0IiKRo0OfIlKYkcBpZpZqZs0IerjGFjH9DKCdmWWEw6cVMt1qoH7M8FzCw5oE55XFrv8sADPrCHQKx39HcKh117CsjpntnmA904Fdi4i3wO7AlBJMJyJS4ZSoiUhh3gEmAROBL4Dr3P33wiZ29/XAJcDHZjYKWASsTDDp68C1ZvaDmbUHBgEDzOxboGnMdE8C9cJDntcRJonuvgQ4F3gtLPsO2DPBej4ADiuqgmZWgyCZyypqOhGRZDH3REcWRERKz8zqufua8CrQx4FZ7v5gEuMZBRzn7isKKe8LdHP3Wyo0MBGRElKPmoiUpYvCiwumAg0JrgJNpquBNkWUpwEPVFAsIiKlph41ERERkYhSj5qIiIhIRClRExEREYkoJWoiIiIiEaVETURERCSilKiJiIiIRJQSNREREZGIUqImIiIiElFK1EREREQiSomaiIiISEQpURMRERGJKCVqIiIiIhGlRE1EREQkopSoiYiIiESUEjURERGRiFKiJiIiIhJRStREREREIkqJmoiIiEhEKVETERERiSglaiIiIiIRpURNREREJKKUqImIiIhElBI1ERERkYhSoiYiIiISUUrURERERCJKiZqIiIhIRClRExEREYkoJWoiIiIiEaVETURERCSilKiJiIiIRFRasgOImqa2k29k45YjzbD4CW2rN1tMv/W4oqYvSVnMiASTFV9W2vWAF4wsaplFzE8R83uR85Vk2QVltnlZpZqviJiKm7ekyyjJckqwnjJrqwTlCduuJOsp4XReFjGXaDov+bJKPU2CVirJR8Ni3xa6lyb8uihYb2HfOxY3XXHLi11/YV8DCaeJHTbfelyCebeOwUs1/VZlMcvZKi5LFHMJ1he7zELaofj5t56uqJi2btNE64mfJsF6Nu0DFjdd4uXE//Wyrf6FcZOyP3H3XkgkKVGLs5GN7JdyNACWUvCJSNniPQCbhg1LSdn0PiiLGQ6nswTzJZx+qzKLK0tJMC7m/7gy37ScxNNvUV7CMreYhLBgXMw0HldWMK2bbVqWb1UWu3y2KCPFthoXu45SlaVsHo6PM4ghfhmUoMwKLdtyPVtPk2j6rdZH6Za5rbEUOt82LctLPN/m/33rxG6L9vBC5tu6zGOWuekPW+y4TfNtquEW89sW88WXbV6Gbfo/0XSb/8BagukBUsxLVJZCEWVx86eYx0y/5TQp5lu8j112wrKYaTaPyy+2LDUmhvjpU/Gtx21aVn7MdPHT5G9aZuqm+fO3WmZq3HpTEyxz0/x4zPvNsW+eLz9uvs1xpsaVxca2KYaY9aWyZZtuWmZM+xVMs3namPcQLiv8H4sZZ+E0FlNmcWUpm6ZJtZQtxqW2mNUUiSwd+hQRERGJKCVqIiIiIhGlRE1EREQkopSoiYiIiESUEjURERGRiFKiJiIiIhJRStREREREIkqJmoiIiEhEKVETERERiSglaiIiIiIRpURNREREJKKUqImIiIhElBI1ERERkYhSoiYiIiISUUrURERERCJKiZqIiIhIRClRExEREYkoJWoiIiIiEaVETURERCSilKiJiIiIRJQSNREREZGIMndPdgyRYmYfA02THUc5aQosTXYQlYTaquTUViWntio5tVXJbW9bLXX3XmUVjJQtJWrViJlluXtmsuOoDNRWJae2Kjm1VcmprUpObVW16dCniIiISEQpURMRERGJKCVq1cvgZAdQiaitSk5tVXJqq5JTW5Wc2qoK0zlqIiIiIhGlHjURERGRiFKiVomZWWsz+9LMppvZVDO7Iq78GjNzM2saM+5GM/vJzGaa2Z9ixnc3s8lh2SNmZhVZl/JWVFuZ2d/D9phqZvfFjFdbxbSVmXUxs+/MbIKZZZlZj5h5qmtbpZvZWDObGLbV7eH4HcxsuJnNCv9vHDOP2mrLtrrfzGaY2SQze8fMGsXMo7aKaauYcn23VyfurlclfQEtgG7h+/rAj0CHcLg18AkwD2gajusATARqAW2B2UBqWDYW2B8w4CPg2GTXryLaCjgc+AyoFZY1V1sV2lafFtQV6A2MUFthQL3wfQ1gDNATuA+4IRx/A3Cv2qrQtjoGSAvH36u2KrytwmF9t1ezl3rUKjF3X+ju48P3q4HpQMuw+EHgOiD2JMQTgNfdPdvd5wA/AT3MrAXQwN1He/DJfgk4sYKqUSGKaKsBwD3unh2WLQ5nUVuxVVs50CCcrCHwW/i+OreVu/uacLBG+HKCNnkxHP8im+uttgpsait3/9Tdc8Px3wGtwvdqq0DsfgX6bq92lKhVEWaWAXQFxphZH2CBu0+Mm6wl8GvM8PxwXMvwffz4Kim2rYDdgYPNbIyZfWVm+4aTqa3Yqq2uBO43s1+BQcCN4WTVuq3MLNXMJgCLgeHuPgbY0d0XQpD4As3DydVWW7dVrPMJen1AbbVVW+m7vXpSolYFmFk94C2CP6S5wE3ArYkmTTDOixhf5cS2lbuvAtKAxgSHYK4F3gjP4VBbbd1WA4Cr3L01cBXwn4JJE8xebdrK3fPcvQtBT1APM+tYxORqq0LaysxuIvj+eqVgVKJFFDG+SknQVp3Qd3u1pEStkjOzGgR/TF9x97eB9gTnKEw0s7kEH/LxZrYTwa+p1jGztyI4fDWfzYcbYsdXKQnaCoK6vx0eahgL5BM8N09ttXVb9QMK3g8FCi4mqNZtVcDdVwAjgF7AovCwE+H/BYfU1VZs1VaYWT/gOOCs8BAdqK2ALdrqBPTdXi0pUavEwp6f/wDT3f3fAO4+2d2bu3uGu2cQfFC7ufvvwDDgdDOrZWZtgd2AseGhmdVm1jNc5jnAe8moU3lJ1Fahd4Ejwml2B2oSPNxYbbV1W/0GHBq+PwKYFb6vzm3VrOAqRTOrDRwFzCBok37hZP3YXG+1FVu2lZn1Aq4H+rj7uphZ1FZs0VY/6Lu9ekpLdgCyXQ4E/gpMDs9lAPg/d/8w0cTuPtXM3gCmERxiuNTd88LiAcALQG2Cc0Q+SrSMSixhWwHPAc+Z2RRgI9Av/EWvttq6rS4CHjazNGADcDFU+/2qBfCimaUS/PB9w93fN7PRBIfRLwB+AU4FtRWJ2+ongqsVh4d3jvjO3furrbZuq8ImruZtVeXpyQQiIiIiEaVDnyIiIiIRpURNREREJKKUqImIiIhElBI1ERERkYhSoiYiIiISUUrURERERCJKiZpINWJmfc3MzWzPmHEZ4X3kMLMuZtY7pqyPmd0Qvj/RzDpswzpHmFlmKedJM7OlZnZ33Ph6Zvakmc02sx/MbJyZXRRfj2KWfaWZnVOKWGqa2cjw/nEiIhVKiZpI9XIGMAo4vZDyLsCmRM3dh7n7PeHgiUCpE7VtdAwwE/hLeEf1As8Cy4Hd3L0rwSOIdijpQsNk63zg1ZLO4+4bgc+B00o6j4hIWVGiJlJNhA9ZPxC4gASJmpnVBAYCp5nZBDM7zczONbPHzOwAoA9wf1jWPranzMyahs8fxMxqm9nrZjbJzIYQ3BG9YB3HmNloMxtvZkPDmBI5A3iY4K7+PcN52xM8X/Rmd88HcPcl7n5vKZrhCGC8u+eGyxxhZg+GPWbTzWxfM3vbzGaZ2Z0x870LnFWK9YiIlAklaiLVx4nAx+7+I/CHmXWLLQx7jm4Fhrh7F3cfElP2LcHzBK8Ny2YXsZ4BwDp37wTcBXSHIJkDbgaOcvduQBbwj/iZw2cbHgm8D7xGkLQB7A1MLEjSttGBwLi4cRvd/RDgKYLnIF4KdATONbMm4TRTgH23Y70iIttEiZpI9XEG8Hr4/nU2J0Bl7RDgvwDuPgmYFI7vSXDo9JvwGaL9gF0SzH8c8GX4gO63gL7hMw+3YGY3hb17v5UithbAkrhxw8L/JwNT3X2hu2cDPwOtw3rkARvNrH4p1iUist10cqxINRD2DB0BdDQzB1IBN7PrtmOxuWz+sZceV5boIcIGDHf34hLEM4ADCw6lAk2AwwkeON3ZzFLcPd/d7wLuMrM1pYh5fYJYs8P/82PeFwzHfkfWIngYvYhIhVGPmkj1cArwkrvv4u4Z7t4amAMcFDfdaqCwXqP4srmEhzXD5RcYSXg+l5l1BDqF478jSMB2DcvqmNnusSswswZhTG3CODMIDkWe4e4/ERwuvbOgh83M0gkSwJKaDuxaiukL4moCLHH3nNLOKyKyPZSoiVQPZwDvxI17CzgzbtyXQIeCiwniyl4Hrg1vi9EeGAQMMLNvgaYx0z0J1DOzScB1wFgITvwHzgVeC8u+A/ZkSycBX4SHHgu8B/Qxs1rAhQQ9bD+Z2TjgM+D6mGn3MLP5Ma9T45b/EcGh2dI6HPhwG+YTEdku5p7oCIWISNVkZu8A17n7rFLM8zZwo7vPLL/IRES2ph41EalubiC4qKBEwtuWvKskTUSSQT1qIiIiIhGlHjURERGRiFKiJiIiIhJRStREREREIkqJmoiIiEhEKVETERERiaj/B0a89cSVeXuHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>longitude</th>\n",
       "      <th>latitude</th>\n",
       "      <th>surface_el</th>\n",
       "      <th>geometry</th>\n",
       "      <th>elev_AGL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>745</th>\n",
       "      <td>-106.173019</td>\n",
       "      <td>35.701873</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17302 35.70187)</td>\n",
       "      <td>2426.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>744</th>\n",
       "      <td>-106.173184</td>\n",
       "      <td>35.701974</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17318 35.70197)</td>\n",
       "      <td>2426.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>743</th>\n",
       "      <td>-106.173656</td>\n",
       "      <td>35.702265</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17366 35.70226)</td>\n",
       "      <td>2426.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>742</th>\n",
       "      <td>-106.174128</td>\n",
       "      <td>35.702555</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17413 35.70256)</td>\n",
       "      <td>2426.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>741</th>\n",
       "      <td>-106.174600</td>\n",
       "      <td>35.702846</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17460 35.70285)</td>\n",
       "      <td>2426.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-106.549319</td>\n",
       "      <td>35.847982</td>\n",
       "      <td>3196</td>\n",
       "      <td>POINT (-106.54932 35.84798)</td>\n",
       "      <td>3439.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-106.549865</td>\n",
       "      <td>35.848082</td>\n",
       "      <td>3198</td>\n",
       "      <td>POINT (-106.54986 35.84808)</td>\n",
       "      <td>3441.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-106.550410</td>\n",
       "      <td>35.848182</td>\n",
       "      <td>3198</td>\n",
       "      <td>POINT (-106.55041 35.84818)</td>\n",
       "      <td>3441.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-106.550955</td>\n",
       "      <td>35.848282</td>\n",
       "      <td>3200</td>\n",
       "      <td>POINT (-106.55096 35.84828)</td>\n",
       "      <td>3443.839992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-106.551500</td>\n",
       "      <td>35.848382</td>\n",
       "      <td>3200</td>\n",
       "      <td>POINT (-106.55150 35.84838)</td>\n",
       "      <td>3443.839992</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>746 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      longitude   latitude  surface_el                     geometry  \\\n",
       "745 -106.173019  35.701873        2183  POINT (-106.17302 35.70187)   \n",
       "744 -106.173184  35.701974        2183  POINT (-106.17318 35.70197)   \n",
       "743 -106.173656  35.702265        2183  POINT (-106.17366 35.70226)   \n",
       "742 -106.174128  35.702555        2183  POINT (-106.17413 35.70256)   \n",
       "741 -106.174600  35.702846        2183  POINT (-106.17460 35.70285)   \n",
       "..          ...        ...         ...                          ...   \n",
       "4   -106.549319  35.847982        3196  POINT (-106.54932 35.84798)   \n",
       "3   -106.549865  35.848082        3198  POINT (-106.54986 35.84808)   \n",
       "2   -106.550410  35.848182        3198  POINT (-106.55041 35.84818)   \n",
       "1   -106.550955  35.848282        3200  POINT (-106.55096 35.84828)   \n",
       "0   -106.551500  35.848382        3200  POINT (-106.55150 35.84838)   \n",
       "\n",
       "        elev_AGL  \n",
       "745  2426.839992  \n",
       "744  2426.839992  \n",
       "743  2426.839992  \n",
       "742  2426.839992  \n",
       "741  2426.839992  \n",
       "..           ...  \n",
       "4    3439.839992  \n",
       "3    3441.839992  \n",
       "2    3441.839992  \n",
       "1    3443.839992  \n",
       "0    3443.839992  \n",
       "\n",
       "[746 rows x 5 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use the file selector to find the ESRI shapefile you'd like to process (.shp)\n",
    "track_path = gui_fname(guidance=\"Select ESRI shapefile (.shp) of track...\")\n",
    "\n",
    "# load the track, plot it for confirmation, and print the attribute table\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 9))\n",
    "track = gpd.read_file(track_path)\n",
    "\n",
    "# add a series with the altitude AGL (but keep it in meters for computation's sake)\n",
    "track[\"elev_AGL\"] = track[\"surface_el\"] + (AGL/3.28084) \n",
    "\n",
    "track = track.iloc[::-1] # to reverse track direction, uncomment\n",
    "\n",
    "# plot the points colored by elevation (scaled relative to the current track)\n",
    "track.plot(column=\"elev_AGL\", markersize=1, \n",
    "           ax=ax, label=\"GPS points\", cmap=\"viridis\", \n",
    "           legend=True,\n",
    "           legend_kwds={'label': \"Altitude AGL (m)\",\n",
    "                         'orientation': \"horizontal\"})\n",
    "\n",
    "# plot the start of the track as a large red point\n",
    "track.head(1).plot(markersize=35, color=\"red\", ax=ax,\n",
    "                   zorder=5, label=\"track start\")\n",
    "\n",
    "# the default aspect tends to stretch everything horizontally... adjust\n",
    "ax.set_aspect(1.7)\n",
    "\n",
    "# label the plot\n",
    "plt.title(track_path, loc=\"left\", fontsize=11, y=1.04)\n",
    "plt.xlabel(\"longitude (°)\", labelpad=15)\n",
    "plt.ylabel(\"latitude (°)\", labelpad=15)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "track"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 4] Add all NMSIM prerequisite calculations to the GeoDataFrame\n",
    "This allows a preview of results before saving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>longitude</th>\n",
       "      <th>latitude</th>\n",
       "      <th>surface_el</th>\n",
       "      <th>geometry</th>\n",
       "      <th>elev_AGL</th>\n",
       "      <th>long_UTM</th>\n",
       "      <th>lat_UTM</th>\n",
       "      <th>knots</th>\n",
       "      <th>climb_angle</th>\n",
       "      <th>time_elapsed</th>\n",
       "      <th>heading</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>745</th>\n",
       "      <td>-106.173019</td>\n",
       "      <td>35.701873</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17302 35.70187)</td>\n",
       "      <td>2426.839992</td>\n",
       "      <td>393879.684606</td>\n",
       "      <td>3.951517e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.012023</td>\n",
       "      <td>-52.273734</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>744</th>\n",
       "      <td>-106.173184</td>\n",
       "      <td>35.701974</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17318 35.70197)</td>\n",
       "      <td>2426.839992</td>\n",
       "      <td>393864.953629</td>\n",
       "      <td>3.951528e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.408696</td>\n",
       "      <td>-52.273506</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>743</th>\n",
       "      <td>-106.173656</td>\n",
       "      <td>35.702265</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17366 35.70226)</td>\n",
       "      <td>2426.839992</td>\n",
       "      <td>393822.638210</td>\n",
       "      <td>3.951561e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>6.805371</td>\n",
       "      <td>-52.273182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>742</th>\n",
       "      <td>-106.174128</td>\n",
       "      <td>35.702555</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17413 35.70256)</td>\n",
       "      <td>2426.839992</td>\n",
       "      <td>393780.323011</td>\n",
       "      <td>3.951593e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.202049</td>\n",
       "      <td>-52.272648</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>741</th>\n",
       "      <td>-106.174600</td>\n",
       "      <td>35.702846</td>\n",
       "      <td>2183</td>\n",
       "      <td>POINT (-106.17460 35.70285)</td>\n",
       "      <td>2426.839992</td>\n",
       "      <td>393738.008215</td>\n",
       "      <td>3.951626e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>9.598731</td>\n",
       "      <td>-52.272361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-106.549319</td>\n",
       "      <td>35.847982</td>\n",
       "      <td>3196</td>\n",
       "      <td>POINT (-106.54932 35.84798)</td>\n",
       "      <td>3439.839992</td>\n",
       "      <td>360090.619154</td>\n",
       "      <td>3.968196e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>2.269272</td>\n",
       "      <td>1068.218217</td>\n",
       "      <td>-76.425576</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-106.549865</td>\n",
       "      <td>35.848082</td>\n",
       "      <td>3198</td>\n",
       "      <td>POINT (-106.54986 35.84808)</td>\n",
       "      <td>3441.839992</td>\n",
       "      <td>360041.558346</td>\n",
       "      <td>3.968208e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1069.698726</td>\n",
       "      <td>-76.425362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-106.550410</td>\n",
       "      <td>35.848182</td>\n",
       "      <td>3198</td>\n",
       "      <td>POINT (-106.55041 35.84818)</td>\n",
       "      <td>3441.839992</td>\n",
       "      <td>359992.497656</td>\n",
       "      <td>3.968219e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>2.269280</td>\n",
       "      <td>1071.178079</td>\n",
       "      <td>-76.424879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-106.550955</td>\n",
       "      <td>35.848282</td>\n",
       "      <td>3200</td>\n",
       "      <td>POINT (-106.55096 35.84828)</td>\n",
       "      <td>3443.839992</td>\n",
       "      <td>359943.437179</td>\n",
       "      <td>3.968231e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1072.658592</td>\n",
       "      <td>-76.424688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-106.551500</td>\n",
       "      <td>35.848382</td>\n",
       "      <td>3200</td>\n",
       "      <td>POINT (-106.55150 35.84838)</td>\n",
       "      <td>3443.839992</td>\n",
       "      <td>359894.376729</td>\n",
       "      <td>3.968243e+06</td>\n",
       "      <td>145.248001</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>746 rows × 11 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      longitude   latitude  surface_el                     geometry  \\\n",
       "745 -106.173019  35.701873        2183  POINT (-106.17302 35.70187)   \n",
       "744 -106.173184  35.701974        2183  POINT (-106.17318 35.70197)   \n",
       "743 -106.173656  35.702265        2183  POINT (-106.17366 35.70226)   \n",
       "742 -106.174128  35.702555        2183  POINT (-106.17413 35.70256)   \n",
       "741 -106.174600  35.702846        2183  POINT (-106.17460 35.70285)   \n",
       "..          ...        ...         ...                          ...   \n",
       "4   -106.549319  35.847982        3196  POINT (-106.54932 35.84798)   \n",
       "3   -106.549865  35.848082        3198  POINT (-106.54986 35.84808)   \n",
       "2   -106.550410  35.848182        3198  POINT (-106.55041 35.84818)   \n",
       "1   -106.550955  35.848282        3200  POINT (-106.55096 35.84828)   \n",
       "0   -106.551500  35.848382        3200  POINT (-106.55150 35.84838)   \n",
       "\n",
       "        elev_AGL       long_UTM       lat_UTM       knots  climb_angle  \\\n",
       "745  2426.839992  393879.684606  3.951517e+06  145.248001     0.000000   \n",
       "744  2426.839992  393864.953629  3.951528e+06  145.248001     0.000000   \n",
       "743  2426.839992  393822.638210  3.951561e+06  145.248001     0.000000   \n",
       "742  2426.839992  393780.323011  3.951593e+06  145.248001     0.000000   \n",
       "741  2426.839992  393738.008215  3.951626e+06  145.248001     0.000000   \n",
       "..           ...            ...           ...         ...          ...   \n",
       "4    3439.839992  360090.619154  3.968196e+06  145.248001     2.269272   \n",
       "3    3441.839992  360041.558346  3.968208e+06  145.248001     0.000000   \n",
       "2    3441.839992  359992.497656  3.968219e+06  145.248001     2.269280   \n",
       "1    3443.839992  359943.437179  3.968231e+06  145.248001     0.000000   \n",
       "0    3443.839992  359894.376729  3.968243e+06  145.248001          NaN   \n",
       "\n",
       "     time_elapsed    heading  \n",
       "745      4.012023 -52.273734  \n",
       "744      5.408696 -52.273506  \n",
       "743      6.805371 -52.273182  \n",
       "742      8.202049 -52.272648  \n",
       "741      9.598731 -52.272361  \n",
       "..            ...        ...  \n",
       "4     1068.218217 -76.425576  \n",
       "3     1069.698726 -76.425362  \n",
       "2     1071.178079 -76.424879  \n",
       "1     1072.658592 -76.424688  \n",
       "0             NaN        NaN  \n",
       "\n",
       "[746 rows x 11 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# find the UTM zone from the elevation (.flt) file generated by `NMSIM Create Base Layers`\n",
    "zone = get_utm_zone(project_dir)\n",
    "\n",
    "# epsg codes for UTM zones in the United States (haven't tested Hawaii yet!)\n",
    "epsg_lookup = {1:'epsg:26901', 2:'epsg:26902', 3:'epsg:26903', 4:'epsg:26904', 5:'epsg:26905', \n",
    "               6:'epsg:26906', 7:'epsg:26907', 8:'epsg:26908', 9:'epsg:26909', 10:'epsg:26910', \n",
    "               11:'epsg:26911', 12:'epsg:26912', 13:'epsg:26913', 14:'epsg:26914', 15:'epsg:26915',\n",
    "               16:'epsg:26916', 17:'epsg:26917'}\n",
    "\n",
    "# convert from D.d (WGS84) to meters (NAD83)\n",
    "projector = pyproj.Transformer.from_crs('epsg:4326', epsg_lookup[zone])\n",
    "\n",
    "# convert into the NMSIM project's coordinate system\n",
    "long, lat = projector.transform(track[\"latitude\"].values, track[\"longitude\"].values)\n",
    "\n",
    "# write to the geopandas dataframe\n",
    "track[\"long_UTM\"] = long\n",
    "track[\"lat_UTM\"] = lat\n",
    "\n",
    "# convert the p(x, y, z) data into an array of vectors\n",
    "V = np.diff(np.array([[x, y, z] for x, y, z in zip(track['long_UTM'], \n",
    "                                                   track['lat_UTM'], \n",
    "                                                   track['elev_AGL'])]), axis=0) \n",
    "\n",
    "# headings using the x and y components of each vector; add np.nan to the end\n",
    "headings = np.degrees(np.arctan2(V.T[0], V.T[1]))\n",
    "headings = np.append(headings, np.nan)\n",
    "\n",
    "# a unit normal vector perpendicular to the xy plane\n",
    "n = np.array([0,0,1])\n",
    "\n",
    "# climb angle using the vector and the unit normal in the xy plane; add np.nan to the end\n",
    "climb_angles = np.array([climb_angle(v) for v in V])\n",
    "climb_angles = np.append(climb_angles, np.nan)\n",
    "\n",
    "# compute the time elapsed given cruising velocity of the aircraft in question; add np.nan to the end\n",
    "time_elapsed = np.cumsum(velocity/np.array([np.linalg.norm(v) for v in V]))\n",
    "time_elapsed = np.append(time_elapsed, np.nan)\n",
    "\n",
    "# convert airspeed to knots and write it to the geodataframe\n",
    "track[\"knots\"] = np.full(track[\"longitude\"].shape, 1.94384*velocity, dtype=None, order='C')\n",
    "\n",
    "# write the remaining data to the geodataframe\n",
    "track[\"climb_angle\"] = climb_angles\n",
    "track[\"time_elapsed\"] = time_elapsed\n",
    "track[\"heading\"] = headings\n",
    "\n",
    "track"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 5] Input an meaningful identifer for this modelling scenario\n",
    "For the **file name suffix** *and* a row in **the metadata header.** <br><br>\n",
    "For example:\n",
    "\n",
    ">\"N74PS beginning 20190610 173326 UTC\" <br><br>\n",
    ">\"C-185 ER_N 800ftAGL\" <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ER_N_GPSsim.shp 800ftAGL\n",
      "C-185_ER_N_800ft\n"
     ]
    }
   ],
   "source": [
    "print(os.path.basename(track_path), str(AGL)+\"ftAGL\")\n",
    "\n",
    "id_string = input()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 6] Save the trajectory file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing C:\\Users\\DBetchkal\\Desktop\\NMSIM_2014_local\\Data\\BAND002\\Input_Data\\03_TRAJECTORY\\C-185 WR_S 800ftAGL.trj\n",
      "\t\t\t...finished writing .trj \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# path to the specific .trj file to be written\n",
    "trj_path = project_dir + os.sep + r\"Input_Data\\03_TRAJECTORY\" + os.sep + id_string + \".trj\"\n",
    "\n",
    "print(\"writing\", trj_path)\n",
    "\n",
    "with open(trj_path, 'w') as trajectory:\n",
    "\n",
    "    # write the header information [exactly as NMSIM expects!]\n",
    "    trajectory.write(\"Flight track trajectory variable description:\\n\")\n",
    "    trajectory.write(\" time - time in seconds from the reference time\\n\")\n",
    "    trajectory.write(\" Xpos - x coordinate (UTM)\\n\")\n",
    "    trajectory.write(\" Ypos - y coordinate (UTM)\\n\")\n",
    "    trajectory.write(\" UTM Zone  \"+str(zone)+\"\\n\")\n",
    "    trajectory.write(\" Zpos - z coordinate in meters MSL\\n\")\n",
    "    trajectory.write(\" heading - aircraft compass bearing in degrees\\n\")\n",
    "    trajectory.write(\" climbANG - aircraft climb angle in degrees\\n\")\n",
    "    trajectory.write(\" vel - aircraft velocity in knots\\n\")\n",
    "    trajectory.write(\" power - % engine power\\n\")\n",
    "    trajectory.write(\" roll - bank angle (right wing down), degrees\\n\")\n",
    "    trajectory.write(\"FLIGHT \" + id_string + \"\\n\")\n",
    "    trajectory.write(\"TEMP.  59.0\\n\")\n",
    "    trajectory.write(\"Humid.  70.0\\n\")\n",
    "    trajectory.write(\"\\n\")\n",
    "    trajectory.write(\"         time(s)        Xpos           Ypos           Zpos         heading        climbANG       Vel            power          rol\\n\")\n",
    "\n",
    "    # now write the data section row by row \n",
    "    # (except for the last row - it's null due to the differencing)\n",
    "    for ind, point in track.iloc[:-1].iterrows():\n",
    "\n",
    "        # write the line\n",
    "        trajectory.write(\"{0:15.3f}\".format(point[\"time_elapsed\"]) + \\\n",
    "                         \"{0:15.3f}\".format(point[\"long_UTM\"]) + \\\n",
    "                         \"{0:15.3f}\".format(point[\"lat_UTM\"]) + \\\n",
    "                         \"{0:15.3f}\".format(0.3048*point[\"elev_AGL\"]) + \\\n",
    "                         \"{0:15.3f}\".format(point[\"heading\"]) + \\\n",
    "                         \"{0:15.3f}\".format(point[\"climb_angle\"]) + \\\n",
    "                         \"{0:15.3f}\".format(point[\"knots\"]) + \\\n",
    "                         \"{0:15.3f}\".format(95) + \\\n",
    "                         \"{0:15.3f}\".format(0) + \"\\n\")\n",
    "\n",
    "    print(\"\\t\\t\\t...finished writing .trj\", \"\\n\")\n"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}