EM1000-Numpy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "import parmap\n",
    "import numpy as np\n",
    "from pyproj import Proj, transform\n",
    "from IPython.core.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mb_reader import readEM1000, gen_input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "st_time = dt.datetime.now()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 129/129 [00:03<00:00, 35.49it/s]\n"
     ]
    }
   ],
   "source": [
    "em_files = list(gen_input(directory='/home/epinux/GreatSouthChannel/', \n",
    "                          output_type='numpy.array',\n",
    "                          reproject=True,\n",
    "                          in_srs='EPSG:4326',\n",
    "                          out_srs='EPSG:32619'))\n",
    "\n",
    "per_beam = parmap.map(readEM1000, \n",
    "                      em_files, \n",
    "                      pm_processes=48, \n",
    "                      pm_pbar=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_data = np.concatenate(([_ for _ in per_beam]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# em_data.view('f8,f8,f8,f8,f4').sort(order=['f3'], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "* **Transform the Longitude and Latitude to WGS85 UTM 19 N** `EPSG:32619`\n",
    "\n",
    "No longer needed, transformation is now performed by `PDAL` during data reading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#wgs84 = Proj(\"+init=EPSG:4326\")\n",
    "#UTM19N = Proj(\"+init=EPSG:32619\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#XUTM, YUTM = transform(wgs84, UTM19N, em_data['X'], em_data['Y'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Compute the index for the nadir beam of each ping**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Beam Count | Nadir Index   | integer resulting from: |\n",
    "|------|------|------------------------------------------------|\n",
    "|   54  | 27| $\\frac{54}{2}$ |\n",
    "|   56  | 82| $54+\\frac{56}{2}$ |\n",
    "|   59  | 139| $54+56+\\frac{59}{2}$ |\n",
    "|   59  | 198| $54+56+59+\\frac{59}{2}$ |\n",
    "|   58  | 257| $54+56+59+59+\\frac{58}{2}$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Compute the euclidean distance between `each beam` from the beam at the `nadir`**\n",
    "\n",
    "* **Compute the incidence angle for each beam**\n",
    "\n",
    "The incidence angle $\\theta$ is compute by:\n",
    "\n",
    "$$\n",
    "\\theta = atan_2{\\frac{z}{d}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAGkCAYAAADaJunmAAAYWGlDQ1BJQ0MgUHJvZmlsZQAAWIWVeQc4lt8b/3med7+89t6y997Ze++Zymtvem2VhGRUkpGMFCpKtENCVEqyslKkEEqlgYT8H6N+39/3d13///U/13We5/Pe5z73OOs+9/MCwJlIDg8PhukACAmNpNgZ6/G5uLrx4SYADEiACmAAkewVEa5rY2MBkPLn/d9lcRBAG+8XUhuy/rf9/1rovX0ivACAbBDs6R3hFYLgWwCgU7zCKZEAYFURukBMZPgGdkcwEwUxEMHhG9hvC6dsYM8tXLDJ42Cnj+BqAPDUZDLFDwCaOoTOF+3lh8ihGUbaGEK9A0IR1jkEa3n5k70B4JREeCRDQsI2sAuCRT3/Icfvv2R6/pVJJvv9xVu+bBa8QUBEeDA57v9zOP7fJSQ46o8OYaRS+1NM7DZ8RsZtOCjMfANTI3gu1NPKGsEMCP4Z4L3Jj2CY6B9l4rjFD3N5RegjYwZYECzrTTYwRzAXgo1Cg60stumevgFGpghGVggcGxBp6rDdN80nwtB+W2YxJczO+g/2pejrbvetIVM29W7wP4wKctTdlj/s72P6R/6PeH8H5y2bUcToACcrBNMgmCUiyN58iwclGO+vb/WHhxJlt2G/IILVfUKN9bbko/b4UozstvkpIRF//EWl+QeYWm3jwkh/B5NtOdVe5E372RBc5xOq6/hHjk+Ei8UfX7x9DAy3fEf1+IQ6bvuLGg+P1LPb7vstPNhmmx9N9Ak23qDvQDBXRLT9dl+0ViSyILfko63CI20ctuxEewaSzWy27EHHAgugDwwAH4hCqicIA4EgoGvu7hzya6vFCJABBfgBHyC1TfnTw3mzJRR52oN48AlBPiDibz+9zVYfEI3Q1/5St55SwHezNXqzRxB4j+AQYA6Ckd9Rm71C/2pzApMIJeB/tHshtgYjdaPtf2m6CMVimxL1Ry4f7R9OrCHWAGuCNcKKoTnQWmgNtAXy1EGqPFoVrfbH2v/wY95jejHvMAOYcczLvQFJlH/5wwcswTiiwWjbZ89/+owWRqQqofXQmoh8RDaaBc0BpNCKiCZdtDaiWwmh6m9bvuH9v2X/lw//GPVtPoIsASawEnQIov/uSSNOo/RXysaY/nOEtmz1/Duu+n9b/q1f/x8j7Y28zf/NiUpD3US1ox6gnqIaUXcBH6oZVYfqRN3fwH9X0eTmKvqjzW7TniBETsD/6CNv69wYyQjZK7KzsqtbbZE+sZEbG0w/LDyOEuDnH8mni5z8PnymoV7SknzysnJqAGzEka1j6rvdZnyAWLr/QyMjcUJVHgCi3n9oYcjZUJOHbI0z/6EJI3uXHZF2w84rihK9RUNvPJDIBGiRHcUOeIAAEEX8kQfKQAPoAENgBqyBA3AFe5BR9kfWMwXEgAPgMEgFmeAkyAOFoBSUg0pwFdwAd0EjeAAeg2egBwyAV8jqmQIfwTxYBCsQBOEgEsQIsUO8kBAkAclDqpAWZAhZQHaQK+QB+UGhUBR0AEqGMqFTUCF0HqqCrkP10APoKdQLvYTeQrPQN+gXjIKpYSaYGxaGZWBVWBc2hx3g3bAfvA+Oh1PgE3ABXAZXw3fgB/AzeAAehz/CCyiAokKxoPhRUihVlD7KGuWG8kVRUAmoDFQ+qgxVg2pA5vkFahw1h1pGY9GMaD60FLKCTdCOaC/0PnQC+hi6EF2JvoN+iH6BfoueR//GkDBcGAmMOsYU44Lxw8RgUjH5mIuY25hHyG6awixisVgWrAhWBdmNrthA7H7sMWwJthbbgu3FTmAXcDgcO04Cp4mzxpFxkbhU3BlcNa4Z14ebwv3EU+F58fJ4I7wbPhSfhM/HX8Y34fvw0/gVAh1BiKBOsCZ4E+IIWYQKQgOhmzBFWCHSE0WImkQHYiDxMLGAWEN8RHxN/E5FRbWDSo3KliqAKpGqgOoa1ROqt1TL1AzU4tT61O7UUdQnqC9Rt1C/pP5OIpGESTokN1Ik6QSpitRGGiP9pGGkkaYxpfGmOURTRHOHpo/mMy2BVohWl3YPbTxtPu1N2m7aOToCnTCdPh2ZLoGuiK6ebohugZ6RXo7emj6E/hj9Zfqn9DMMOAZhBkMGb4YUhnKGNoYJRhSjAKM+oxdjMmMF4yPGKSYskwiTKVMgUybTVaYupnlmBmZFZifmWOYi5vvM4ywoFmEWU5ZgliyWGyyDLL9YuVl1WX1Y01lrWPtYl9g42XTYfNgy2GrZBth+sfOxG7IHsWez32Uf5UBziHPYcsRwnOV4xDHHycSpwenFmcF5g3OEC+YS57Lj2s9VztXJtcDNw23MHc59hruNe46HhUeHJ5Anl6eJZ5aXkVeLN4A3l7eZ9wMfM58uXzBfAd9Dvnl+Ln4T/ij+8/xd/Cs7RHY47kjaUbtjVIAooCrgK5Ar0CowL8graCl4QPCK4IgQQUhVyF/otFC70JKwiLCz8FHhu8IzImwipiLxIldEXouSRLVF94mWifaLYcVUxYLESsR6xGFxJXF/8SLxbglYQlkiQKJEolcSI6kmGSpZJjkkRS2lKxUtdUXqrTSLtIV0kvRd6c8ygjJuMtky7TK/ZZVkg2UrZF/JMciZySXJNch9kxeX95Ivku9XICkYKRxSqFP4qiih6KN4VnFYiVHJUumoUqvSmrKKMkW5RnlWRVDFQ6VYZUiVSdVG9ZjqEzWMmp7aIbVGtWV1ZfVI9RvqXzSkNII0LmvM7BTZ6bOzYueE5g5NsuZ5zXEtPi0PrXNa49r82mTtMu13OgI63joXdaZ1xXQDdat1P+vJ6lH0bust6avrH9RvMUAZGBtkGHQZMhg6GhYajhntMPIzumI0b6xkvN+4xQRjYm6SbTJkym3qZVplOm+mYnbQ7KE5tbm9eaH5OwtxC4pFgyVsaWaZY/naSsgq1OquNbA2tc6xHrURsdlnc88Wa2tjW2T73k7O7oBduz2j/V77y/aLDnoOWQ6vHEUdoxxbnWid3J2qnJacDZxPOY+7yLgcdHnmyuEa4FrnhnNzcrvotrDLcFferil3JfdU98HdIrtjdz/dw7EneM/9vbR7yXtvemA8nD0ue6ySrcll5AVPU89iz3kvfa/TXh+9dbxzvWd9NH1O+Uz7avqe8p3x0/TL8Zv11/bP958L0A8oDPgaaBJYGrgUZB10KWg92Dm4NgQf4hFSH8oQGhT6MIwnLDasN1wiPDV8fJ/6vrx98xRzysUIKGJ3RF0kE3Jh74wSjToS9TZaK7oo+meMU8zNWPrY0NjOOPG49LjpeKP4C/vR+732tx7gP3D4wNuDugfPJ0AJngmthwQOpRyaSjROrDxMPBx0+HmSbNKppB/JzskNKdwpiSkTR4yPXEmlSaWkDh3VOFqahk4LSOtKV0g/k/47wzujI1M2Mz9z9ZjXsY7jcscLjq+f8D3RlaWcdfYk9mToycFs7ezKU/Sn4k9N5Fjm3Mnly83I/ZG3N+9pvmJ+6Wni6ajT4wUWBXVnBM+cPLNa6F84UKRXVFvMVZxevFTiXdJ3VudsTSl3aWbpr3MB54bPG5+/UyZcll+OLY8uf1/hVNF+QfVC1UWOi5kX1y6FXhqvtKt8WKVSVXWZ63LWFfhK1JXZavfqnqsGV+tqpGrO17LUZl4D16KufbjucX3whvmN1puqN2tuCd0qvs14O+MOdCfuzvxd/7vjda51vfVm9a0NGg2370nfu9TI31h0n/l+VhOxKaVpvTm+eaElvGXugd+Dida9ra/aXNr6H9o+7Hpk/ujJY6PHbe267c1PNJ80PlV/Wt+h2nH3mfKzO51KnbefKz2/3aXcdadbpbuuR62noXdnb1Ofdt+DFwYvHveb9j8bsBroHXQcHB5yHxof9h6eeRn88utI9MjKq8TXmNcZo3Sj+WNcY2VvxN7UjiuP339r8Lbznf27VxNeEx8nIyZXp1Lek97nT/NOV83IzzTOGs32fNj1Yepj+MeVudRP9J+KP4t+vvVF50vnvMv81FfK1/Vvx76zf7/0Q/FH64LNwthiyOLKUsZP9p+Vy6rL7b+cf02vxKziVgvWxNYafpv/fr0esr4eTqaQN68CKKTCvr4AfLsEAMkVAMYe5E6xayvP2y4o5PIBI28nSBr6CD9EJaPtMTpYERwHno3AS9SksqIOIp2kqaedo5di8GEsZ5pgEWeNY2vmoOV05qrg/s67ky+F/7kAvaCd0HHhZ6JATEHcV+K0ZIfUkoyorK1covwVhQElWFlOZbdqhtod9bc7SZqqWh7a6TrXdV/r4w2UDb2MThrXmYyZQeaCFsaWgVZZ1rdshm1/2rM4KDhaO4U4H3epcX3m9nbXvPvS7pW9wINIZveU8tL1tvPZ6+vjR/a3D9gZyBcEBY0HN4ecC00O8w+32adK4YvAR3yJHIxqiq6MyYlNiAuOd91vekDzoEqC8iG1RN3D5knOyT4pkUeOpOYerUi7md6S0Zk5eOzN8ekTn7K+nVzIXjy1kLOQ+ysffZq5QPKMcaFX0aHigpKas82lz871nx8pGy+frfhxEXWJuVK8Su+y+5WY6tyrN2p6a79ep7+hcNP+VsTtk3eq7jbUPahva2i5d6/x9v3apqrm8paSB3mtGW0HHgY+sn+s3M7Wvvxk/Gl3x+NnbZ0Pnjd21XYX9ET06veR+l68KOr3HVAaxAwODVUOR7/UGcGOtCPrS+n19Gj2mMbYxJvj4xrjH9+WvrObQE3UTjpOLk/lvpd83zxtNz05c2RWZnbyQ+XH0DmFuYVPtZ+9vtB/uT1vM//+64FvrN8ef8/6EbpAXvRF1tHkr0dr0uvrm/MvAF2DA1HyqBn0dUwi1gWniZciiBBFqHZQy5LUaWxpvegS6EsZmhhnmelYVFnJbGnstzjGuKi4FXh28Sbynedv3vFKYEGISphXREnUVMxDPE4iR/K6VKf0jCxajl9+p4KbYqRSpnKFSr3qc7V36j92YjU5teS0LXWCdbP0run3GHwywhtzm8ibGpo5mntZhFrGWiVYJ9scsU21S7PPcDjmmOGU4hzn4u/q4GawS9vdaLfbnpi9eR7XyK2eHV6PvG/7FPvu93P2lw2gDpgL7AlqCK4KKQrNCksKp+xzp+hE8EasRA5EXY1OjfGMNYyTjRfcz32A/SBzAt0h7KHFxHeHO5KuJ+elxBzZnWp21CDNIp2ccTjzwrHHx8dOfM5aOLmUvXDqe8587qe8ufzPp3+eoStUKwotvljSdXaidPbc1Pk3ZS/LeyueXGi62Hipo/LTZf4ru6uLr76sZbpmdT0NOb2Wb0vf8b5bVNfXgLmn2Lj3/pGmi82NLU0PLreebDv4MOZR4uOs9pIn5U/Pdpx4FtVp/1yqC9010n2jJ7M3sM/2hWG/4YDtoOdQ1HDKy6MjB1/5vtYf5RidG6t/c3Tc5a3UO/y79xNtkyVT+97rTFNP98+Uzx76EPDRe87/U8jn8C/h8+FfKd+iv8f9iFkIWDReol26+dPw57Nlt+VPv3pWqddGNudfAjyEzKFh2AeFRWWhJdDdmHisDHYWdwHvT5AhLBM7qEqpY0h2NPK0NLSLdC/pWxiqGHOYDjL7sdixarKJsTOzr3LMcPZxNXHX8JTzFvHl8+fuyBJIFYwWIgsbivCJ/BTtFCsVj5AwkeSXgqVmpYdknsg2yF2WL1BIVPRQUlPGKner5Km6qLGrvVQv0fDeKa+J1RzTuqOdpeOva6AnrE9nAAy+G04bDRrfM8k39TETMhs3L7CwtsRZtlklW5vasNl8sG2yy7H3d9BwJDmOOV11PuBi5srs+satclcYEv+Xd9/fk7hX3wPv0Usu9gzy2ulN7T3ic8l3n5+q36p/c0BioE4QCGoJPhyiH4oOfRR2JFw3/Oe+KxRXJGZXRVpH/ogqiN4ZPRaTGMsdez/OI54lfmT/lQPJB10SRBMWD7Ul5hz2SzJIFk9hO0KVClJ/HJ1Ie55em3Esk3xM8Tju+MiJa1kZJ4OyjU8xnHqcsytnLjc+Tzdf73TaGXxhRtFkCftZ+VK1c2rnlcpkykUr+C+wX6S/RKwkVNEiK0mz2uPq0ZqrtS+urd4Qvel269Tt3rtMda71xQ1DjZj7Yk3GzZ4thx6cbW1qe/Nw/TF/u/4Tv6fHOq4/G+xc6xLr3tVzunfshXz/8YHPQ/bD9SP8r/JGZd7QvI2ZzJyJ+2T1bXHZdmP+t773bRSsMgA5SJ7pdBypswBk30XyzHsAsBIBsCEB4KAG4KM1ADauAVDQkb/xA0ISTzySc7IAXiAGFJFM0wK4IVlzLEhHMspq0AT6wHuwCjFAYpAOkh9GQMeRfPARNAFDMD+sB3vDR5Esrw/+hRJAWaLiUZWoITQerY4OQZejX2IYMOZIRtaGhbA62ERsKw6DM8OdxA3j+fHB+HoCjuBMqCT8IloSzxOXqKyoKqnR1J7UbSQhUjrpM40DTSOS6WTTAbp9dJP0rvTdDEYM9xlVGe8wqTO1MdsxT7BEsWJZ89mE2erYrdhnONI45TgnuEq5PXkkeH7yPubL4/feoSiAFXgleFMoSzhYxFxUQowkNi8+IHFP8qxUgrS7jJosk+y83HP5ywrpiv5KZsrSKswq66qf1MbU+zQ6dj7SfKjVrt2lM6I7o7doAAyxyDmHN8GbEsyozZks+C0VraysQ21ybRvtphxIjopOrs4HXc65PnSbdqfaLbvHae8Bjwpyl+dPb0Efe98jfo3+vwL1g84EL4d6hfXtM6I0RipG1cZIxV6P37m/52DYIa7EwaTcFIsji0dz0yUzHh3zOcGc9Sb7ec5o3noBX6FascXZvefiys5VjFySqjpXLVszfv38rT13qeprGnc3S7TyPjJ6UtZJ3S3au9ifPST6svf12Ten3/W995hd/sTwpfob+CG7qLa0vpzxq26lf/XeWvnv8HWVzfMD2vzmwAA4gTCQB9rAEriDEJAAskEFqAfdYAqsQSyQDGQG+ULJUBn0AHoHo2ER2AKmwIVwG/wFxYUyRx1A1aIm0RxoO3Qm+hEGwmhi9mPuYVax2thk7FMcHc4VdwH3Da+Lz8G/J2gQcghzRCNkzlepXKhuIZkwhbqfpEY6R0NFE0szTetK20VnRNdCr0XfzKDP0MFozziKZKa/mLNYxFmese5jY2G7w27L/p4jjpPEWcGlwzXJnc1jxkvDO8p3k//EjgABPUE2wY9C94VPiviK6okJiTNI4CUxUnhpGhkGWXo5vNyy/IzCkGKH0gPlByodqq/UvmnQ7JTVtNUK0I7Uoej667noGxuoGSoaqRobm+w1TTA7b95uMW/FaW1oE4TEtFz70w55jrlO55ybXb66Ke1KdH++h2dvpEe3p4CXr3eez23fLr9J/5VAliCFYIeQ6NDCsJbwDxTWCKPI6KhL0SOxdHGW8Vn7hw8KJxw8NHHYL5kupSM1Mg2bfjQTfSztBGdWW3ZSjkue/mmNMxpFGiVqpWLn0WWPK6Ivcl66X+V5hbl6tObRte4bC7fl7h6of9ZI22TQQmm9+HC2Xe/pjU65ruKe0b4f/V8Hp4cnRmZe/3gDvSVOME0JTpvM5s+pfMn4fnEpeLlrJWW1be3H7+XN+YeR3U8PeIAU0AK2wBccBPngGugEHyACJAFZQhSoAGqBPsAssAEcCV+ER1D0KFNUCqoFtYbWQMejG9CrGF1MBmYIK4Y9jB3FaeHK8Hh8GL6foEYoIcLEQOIAlQHVPWo16gckG9J7miRaftoWOne6RfqTDFIMzxlDmUhMlcx6zK9Z4lh5WLvYTrB7cuhxinMxca1wj/LU8Z7iC+G32CErwCaIFVwW+ir8ReS76Jo4jYSgpI6Uh3SiTIlsndwL+e+KHEqmykkqbWrU6u4a1zRxyF21SXeHXo4Bi2GNsZspvVmvRaFVmI2jnbz9iKObU6eLieuLXb7uP/cke0DkcM8BbxWfYj+C/+FAYlB5iGUYCL9LCYvkiWqLiYrz3v85oSIx7vBg0moKfASfSndUIS0ivT/T8djsibST0tkvc9LyNPK/FlQV7ikmllwqVTl3v0y7vOWCwcWOSpuq/isO1T01RrX110VvnL6Fv33wzmpdeoPwvZ77Sc3KLbOtxQ+tH6Pb7z2NeCbROdl1tselj+lF30DWkNnw+kj1a+vRmTdR42vvkiZRU0nT8EzyB/THQ3OfPxt9iZsv+Xr8W9R3g+9LPy4vWC28WvRfXFyKXpr96f6ze1l/+cov0q/wX30rSisFK19XTVfLVlfWHNau/kb9dvldvQ6tO65f3pj/CF8F+c3wAVHrAYAZW1//LgwA7hQAa9nr6ytl6+tr5Uiy8RqAluCt/5A2Yw0dAMVvNlCH+EDiv//L+T+h799W3D36ZQAAAZ1pVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDUuNC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQaXhlbFhEaW1lbnNpb24+MzkxPC9leGlmOlBpeGVsWERpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxZRGltZW5zaW9uPjQyMDwvZXhpZjpQaXhlbFlEaW1lbnNpb24+CiAgICAgIDwvcmRmOkRlc2NyaXB0aW9uPgogICA8L3JkZjpSREY+CjwveDp4bXBtZXRhPgqCdjvZAAA4aUlEQVR4Ae2dB3hVVbbHF70nlNAhoYeOgDqOOhQbKlVK6F2aYJv50EHmDbz3AUqxDlakShdpoYMCCQQEURyeBZx5jjNvxlemPKVDgDfrzCSGs09IbnLLKb/9fZp71t7n7L1+64Z/9tmtyLW/JyFBAAIQgAAEchAomuMzHyEAAQhAAAIWAcSBLwIEIAABCBgEEAcDCQYIQAACEEAc+A5AAAIQgIBBAHEwkGCAAAQgAAHEge8ABCAAAQgYBBAHAwkGCEAAAhBAHPgOQAACEICAQQBxMJBggAAEIAABxIHvAAQgAAEIGAQQBwMJBghAAAIQ8I04fPPNN0QTAhCAAATCRMA34tClSxf5y1/+EiYsPAYCEIBAsAn4RhxOnjwp3bp1k3PnzgU7ongPAQhAIAwEfCMOyuLw4cPSv39/uXLlShjQ8AgIQAACwSXgK3HQMG7ZskXGjRsX3IjiOQQgAIEwEPCdOCiThQsXypQpU8KAh0dAAAIQCCYBX4qDhvK5556TmTNnBjOqeA0BCECgkASK+OWY0CJFihgo1LZ69WpJSUkx8jBAAAIQgEDuBHwtDup2yZIlrXGIe++9N3cK5EAAAhCAwHUEfC8O6m1cXJykpaVJmzZtrnOeCwhAAAIQcCYQCHFQ12vWrCmHDh2SpKQkZxJYIQABCEAgm4BvB6SzPfznh2+//Vbuu+8+VlHbwXANAQhAwIFAYMRBfT916pR07dqVVdQOXwRMEIAABHIS8LU41KlTJ6ev1ucPP/xQevXqxSpqgwwGCEAAAj8Q8LU47Ny5UypWrPiDt//8tHv3bhk7dqxhxwABCEAAAv8g4GtxaN68uaSmpkrp0qWNeC9atEh+8YtfGHYMEIAABCAg4mtx0ADfeeedsmrVKila1HRVV1DPmTOH7wEEIAABCNgI+Hoq699Xf2e7+8Ybb8iECROyr7M+qGi8++670rt37ywTPyEAAQgEnoD557RPkYwfP96xl3D16lUZNGiQpKen+9Rz3IIABCAQOoHA9Byy0Dz88MPWrq1Z11k/deD6wIED0qJFiywTPyEAAQgElkDgxEEPAtKprHrugz3Vrl3bWkVdt25dexbXEIAABAJFIHDioNE9f/683HXXXdbJcfZoN2vWTA4ePCiVKlWyZ3ENAQhAIDAEAjPmkDOiZcqUka1bt0pycnJOs/X5iy++kO7du8uFCxeMPAwQgAAEgkIgkOKgwa1cubLs2rVLatWqZcRaew4DBgwQHawmQQACEAgigcCKgwY7MTFRduzYIfHx8UbsN23a5Dj11SiIAQIQgIAPCQRaHDSerVq1ko0bN1qHAtnj+9Zbb8nkyZPtZq4hAAEI+J5A4MVBI9ypUydZvny5OB01Om/ePFmwYIHvvwg4CAEIQCAnAcThnzT69esnr7zySk422Z91ZfXmzZuzr/kAAQhAwO8EEIccEZ40aZLMmDEjh+UfH3VthA5QZ2RkGHkYIAABCPiRQCDXOeQVyBEjRsjSpUuNYjrDSVdR61oIEgQgAAE/E0AcHKKbmZlprXXQmUz2pDOc9Cxqpymw9rJcQwACEPAqAcQhl8idPXtWOnfuLEePHjVKtGzZ0upBOE2BNQpjgAAEIOBBAow55BK0cuXKybZt26Rx48ZGiX//93+XHj16yMWLF408DBCAAAT8QABxuEEUExISRI8arV69ulEqLS1NBg8ezCpqgwwGCEDADwQQhzyiWL9+fdm+fbtUqFDBKPnee+/JY489ZtgxQAACEPA6AcQhHxFs27atbNiwQUqUKGGUfvXVV+WZZ54x7BggAAEIeJkAA9IhRE/PotZXSTmPH826ffHixaJTYEkQgAAE/ECAnkMIURw4cKA8//zzjneMGTPGGsB2zMQIAQhAwGMEEIcQA/bkk0/KtGnTjLt0bYRuwXHkyBEjDwMEIAABrxHgtVIBIqavlYYMGSIrV6407tYZTnoeRJMmTYw8DBCAAAS8QgBxKGCkLl++LF27dpXdu3cbT6hXr561irpGjRpGHgYIQAACXiCAOBQiSmfOnJGOHTvKxx9/bDzlpptuEl0L4TQF1iiMAQIQgIDLCDDmUIiAlC9f3loD0aBBA+Mpx48fl169esmlS5eMPAwQgAAE3E4AcShkhKpVq2atoq5atarxpA8++ECGDRvmOPXVKIwBAhCAgIsIIA5hCEajRo2saay6H5M9rVmzRn7605/azVxDAAIQcDUBxCFM4bn55ptFt9MoXry48cSXXnpJpk+fbtgxQAACEHArAQakwxyZd955x3qV5PRYzdMpsCQIQAACbidAzyHMERo6dKjMnj3b8amjRo2SXbt2OeZhhAAEIOAmAvQcIhSNqVOnyqxZs4yn6wynffv2Sfv27Y08DBCAAATcQgBxiFAkdBX1gAEDZO3atUYNOsMpIyNDGjZsaORhgAAEIOAGAohDBKOgaxzuv/9+2bt3r1GLCoMKhAoFCQIQgIDbCCAOEY7I999/Lx06dJBPP/3UqElfLekrJn3VRIIABCDgJgIMSEc4GnFxcdYq6qSkJKOmY8eOSZ8+fUT3aSJBAAIQcBMBxCEK0ahZs6a1irpKlSpGbTp7SWcxkSAAAQi4iQDiEKVoJCcny9atW6Vs2bJGjcuXL5fJkycbdgwQgAAEYkUAcYgi+R/96EfW7KVixYoZtc6bN09mzpxp2DFAAAIQiAUBBqRjQH3RokUyevRoo+YiRYpYBwjpFFgSBCAAgVgSoOcQA/o6xjBjxgyjZl0bMXz4cNHdXEkQgAAEYkmAnkMM6T/11FMyd+5cowV6QJAeFKQHBpEgAAEIxIIA4hAL6v+s8+rVq9KvXz9Zv3690Qo9YvTQoUOiR46SIAABCESbAOIQbeK2+i5evCj33nuvpKen23JEmjRpIgcPHpSEhAQjDwMEIACBSBJgzCGSdPPx7FKlSsnmzZulZcuWRulTp05J165d5dy5c0YeBghAAAKRJIA4RJJuPp9dsWJF2bFjh9StW9e448iRI9arp8zMTCMPAwQgAIFIEUAcIkU2xOfWrl3bWkVdqVIl485t27bJmDFjDDsGCEAAApEigDhEimwBntusWTNJTU2V0qVLG3cvWbJEnnnmGcOOAQIQgEAkCCAOkaBaiGfecccdsnr1aila1AzNs88+K3PmzCnE07kVAhCAQP4IMFspf5yiXuqtt96ScePGGfWqaOgBQrqbKwkCEIBApAiYf55GqiaeGxKBsWPHOp5FrWsjBg8eLPv37w/peRSGAAQgEAoBeg6h0IpBWRWJBQsWGDXHx8dbayNatWpl5GGAAAQgUFgCiENhCUb4/itXrkjv3r2ttRD2qmrVqmWtok5MTLRncQ0BCECgUAQQh0Lhi87N58+fl3vuucc6c9peY9OmTa1V1JUrV7ZncQ0BCECgwAQYcygwuujdWKZMGWuKq051tacvv/xSunfvLiogJAhAAALhIoA4hItkhJ+jPQNdRa2vkuwpIyND9AwIfQVFggAEIBAOAohDOChG6Rk6trBz507RwWh70v2ZJkyYYDdzDQEIQKBABBCHAmGL3U26QZ8KgW7YZ086q4mzqO1UuIYABApCAHEoCLUY39OhQwdZsWKF4ypqPYtaF9CRIAABCBSGAOJQGHoxvFdXSL/yyiuOLdDXS5s2bXLMwwgBCEAgPwQQh/xQcmmZiRMnOp5FrauodYBaDwoiQQACECgIAdY5FISay+4ZOXKk6K6t9qTbfx84cECaN29uz+IaAhCAwA0JIA43xOONTD0IqGfPnqLnPtiTHiCkZ1HreREkCEAAAvklgDjkl5TLy+lRop07dxY9Oc6eWrRoYfUg9MQ5EgQgAIH8EGDMIT+UPFCmbNmysnXrVmnSpInR2s8++0x69OghFy9eNPIwQAACEHAigDg4UfGoLSEhwVokV6NGDcOD9PR0GTRokOhgNQkCEIBAXgQQh7wIeSy/Xr16sn37dqlQoYLR8vXr18ukSZMMOwYIQAACdgKIg52ID65vuukm2bhxo5QsWdLw5vXXX+csaoMKBghAwE6AAWk7ER9dr1mzRgYOHCjXrl0zvFq4cKGMGjXKsGOAAAQgoAToOfj4e9C/f3954YUXHD3UE+Z0AJsEAQhAwIkA4uBExUe2J554QqZNm2Z4pNt7p6SkyOHDh408DBCAAAR4rRSQ78DQoUNl+fLlhrdVqlSxttlITk428jBAAALBJYA4BCT2ly9flm7dusmuXbsMj5OSkqxV1DVr1jTyMEAAAsEkgDgEKO5nzpyRTp06ybFjxwyv27RpI2lpaRIXF2fkYYAABIJHgDGHAMW8fPny1v5LDRs2NLz+9NNPpVevXnLp0iUjDwMEIBA8AohDwGJerVo1axW1/rSnvXv3io5NOE19tZflGgIQ8DcBxMHf8XX0TnsOuoOr9iTsae3ataIznEgQgECwCSAOAY1/+/bt5b333pMSJUoYBPSEuenTpxt2DBCAQHAIMCAdnFg7eqrTW4cNG+b4KmnZsmXWaybHGzFCAAK+JkDPwdfhzdu5IUOGyOzZsx0L6vYaO3fudMzDCAEI+JsAPQd/xzff3k2dOlVmzZpllC9Xrpzs27dPbr75ZiMPAwQg4F8CiIN/YxuSZzpDSc97WL16tXFf1apVJSMjQxo1amTkYYAABPxJAHHwZ1wL5JWucXjggQfkgw8+MO5v0KCBJRDVq1c38jBAAAL+I4A4+C+mhfLo9OnT0qFDBzl+/LjxnHbt2sn+/fsdp8AahTFAAAKeJsCAtKfDF/7G6wlyepKcnihnTx9//LE89NBDovs0kSAAAX8TQBz8Hd8CeadnUOssJT2T2p727NkjI0aMcJz6ai/LNQQg4F0CiIN3YxfRljdp0sQ6DKhs2bJGPStXrpTJkycbdgwQgIB/CCAO/oll2D259dZb5d1335XixYsbz37++edl5syZhh0DBCDgDwIMSPsjjhH1YsmSJTJy5EijjiJFisiKFSusc6qNTAwQgICnCdBz8HT4otN4HWNw6iXo2ojhw4eLjkOQIAABfxGg5+CveEbUGx1nmDdvnlGHznDSKa5t27Y18jBAAALeJIA4eDNuMWn11atXJSUlxdrN1d4AXRx36NAhqV+/vj2LawhAwIMEEAcPBi2WTb548aJ06dLF6inY29G4cWM5ePCg6HYbJAhAwNsEEAdvxy8mrf/uu+/kJz/5iZw4ccKo/5ZbbhE9UU437CNBAALeJcCAtHdjF7OWx8fHy44dOyQxMdFow9GjR6Vv376SmZlp5GGAAAS8QwBx8E6sXNXSWrVqWauoK1eubLRLhePhhx827BggAAHvEEAcvBMr17W0adOmkpqaKmXKlDHatnTpUpkyZYphxwABCHiDAOLgjTi5tpW33367dQZEsWLFjDY+99xzMmfOHMOOAQIQcD8BBqTdHyNPtHDBggUyduxYo626inrNmjXSr18/Iw8DBCDgXgL0HNwbG0+1bMyYMY5nUesq6qFDh1pHjXrKIRoLgYAToOcQ8C9AuN0fP368vPnmm8ZjdYZTWlqatG7d2sjDAAEIuI8A4uC+mHi6RbqKunfv3rJp0ybDD53hpGdRJyUlGXkYIAABdxFAHNwVD1+05sKFC3LPPfdYq6XtDiUnJ1v2KlWq2LO4hgAEXESAMQcXBcMvTSldurQ1xbV58+aGSydPnpRu3brJ+fPnjTwMEICAewggDu6Jha9aUqlSJWsVde3atQ2/Dh8+bG3gd+XKFSMPAwQg4A4CiIM74uDLVtStW9daRV2xYkXDvy1btsi4ceMMOwYIQMAdBBAHd8TBt61o0aKFbN68WfRVkz0tXLiQs6jtULiGgEsIIA4uCYSfm6E7uOpxokWLml83PTzojTfe8LP7+AYBTxIwf1s96QaNdjsBnd46f/58x2ZOnDhRNmzY4JiHEQIQiA0BxCE23ANZ64QJE2TGjBmG77o2YtCgQXLgwAEjDwMEIBAbAqxziA33QNc6evRoWbRokcFAB65VIHScggQBCMSWAOIQW/6BrF2nsPbs2VO2bt1q+F+nTh3rLGr9SYIABGJHAHGIHftA13zu3Dm5++67Rdc82JMuntMehK6VIEEAArEhwJhDbLgHvtayZcuKrnXQ7TTs6fPPP5cePXqIbsNBggAEYkMAcYgNd2r9OwHdX2nnzp1Ss2ZNg4f2HAYOHCg6WE2CAASiTwBxiD5zasxBQHdo3b59u8TFxeWw/uPjxo0bRae5kiAAgegTQByiz5wabQTatGkjKgQlS5a05Yi1QG7q1KmGHQMEIBBZAgxIR5YvTw+BwNq1a2XAgAGip8fZkx5D+vDDD9vNXEMAAhEiQM8hQmB5bOgEUlJS5KWXXnK8UU+YS01NdczDCAEIhJ8A4hB+pjyxEAQee+wxmTZtmvEEXRvRv39/aw2EkYkBAhAIOwFeK4UdKQ8MB4Hhw4fLsmXLjEdVrlzZOkmuadOmRh4GCEAgfAQQh/Cx5ElhJJCZmWmdGKdTXe0pMTHR6kHomdQkCEAgMgQQh8hw5alhIHD27Fnp1KmTfPTRR8bTWrVqJenp6RIfH2/kYYAABApPgDGHwjPkCREiUK5cOdm2bZs0atTIqOHEiRPW/kwXL1408jBAAAKFJ4A4FJ4hT4gggapVq1qrqKtXr27Usn//fhkyZAirqA0yGCBQeAKIQ+EZ8oQIE2jQoIHVgyhfvrxR07p16+Txxx837BggAIHCEUAcCsePu6NEoF27dtZpcSVKlDBq1BPmpk+fbtgxQAACBSfAgHTB2XFnDAisXLnSepXktIp6yZIlolNgSRCAQOEJ0HMoPEOeEEUCepzo3LlzHWvU7TV0Ez8SBCBQeAL0HArPkCfEgIBuxjdr1iyjZp3htHfvXrnllluMPAwQgED+CSAO+WdFSRcR0NdKgwcPllWrVhmtSkhIkIyMDGncuLGRhwECEMgfAcQhf5wo5UICly9flgcffFD27NljtK5+/fqWQNSoUcPIwwABCORNAHHImxElXEzg9OnT0rFjR/nkk0+MVrZt21Z0LUSFChWMPAwQgMCNCTAgfWM+5LqcgP7Dr4PQ2lOwJxWMhx56SC5dumTP4hoCEMiDAOKQByCy3U9AV0/rBn26mtqe3n//fRkxYoTjAUL2slxDAAI/EEAcfmDBJw8T0MHnrVu3is5WsicdtP7Zz35mN3MNAQjcgADicAM4ZHmLgE5f1e00ihcvbjT8xRdflJkzZxp2DBCAgDMBBqSduWD1MIGlS5dar5LsLhQpUkTeeecdawqsPY9rCEDgegL0HK7nwZUPCOgWGs8++6zhia6NGDlypOzevdvIwwABCFxPgJ7D9Ty48hGByZMny7x58wyPdHdXneKqm/mRIAABZwKIgzMXrD4gcPXqVRkwYIC8++67hjc6w0lXUet24CQIQMAkgDiYTLD4iICeFHf//ffLvn37DK/0hDkVCKcpsEZhDBAIGAHEIWABD6K73333nXTo0EF+/etfG+7ffPPNlnA4TYE1CmOAQIAIMCAdoGAH1dX4+HhrFXVSUpKB4KOPPpI+ffpIZmamkYcBAkEmgDgEOfoB8r1WrVqyY8cOqVKliuG1rq4eNWqUYccAgSATQByCHP2A+d60aVPZsmWLlClTxvBc1z88/fTThh0DBIJKAHEIauQD6vdtt90ma9eulWLFihkE5syZI/ofCQIQEGFAmm9BIAksXLhQ9FhRe9JV1KtXr5aUlBR7FtcQCBQBeg6BCjfOZhEYPXq0zJ49O+sy+6euoh46dKh11Gi2kQ8QCCABeg4BDDou/0BgwoQJ8sYbb/xg+OenuLg4SUtLkzZt2hh5GCAQBAKIQxCijI+5EtBV1H379pUNGzYYZWrWrGktkqtXr56RhwECfieAOPg9wviXJ4ELFy7IvffeKwcOHDDKJicny8GDBx2nwBqFMUDARwQYc/BRMHGlYARKly4tqamp0qJFC+MBJ0+elK5du8q5c+eMPAwQ8DMBxMHP0cW3fBOoWLGitUiuTp06xj0ffvihNXuJVdQGGgw+JoA4+Di4uBYaARUGXS1dqVIl40Y9gnTcuHGGHQME/EoAcfBrZPGrQASaN28umzdvFn3VZE+LFi0SPSOCBIEgEEAcghBlfAyJwJ133imrVq2SokXNXw89POi1114L6XkUhoAXCZjffi96QZshEGYCvXr1kldffdXxqY8++qisX7/eMQ8jBPxCAHHwSyTxI+wExo8fLzNmzDCeq2sjBg0aZC2SMzIxQMAnBFjn4JNA4kbkCIwZM0befvttowKd4ZSeni4tW7Y08jBAwOsEEAevR5D2R5zAlStX5KGHHrLWQtgrq127thw6dEjq1q1rz+IaAp4mgDh4Onw0PloEzp8/L3fffbclBPY6mzVrZq2idpoCay/LNQS8QoAxB69EinbGlIAeEKQHBemBQfb0xRdfSPfu3UUFhAQBvxBAHPwSSfyIOIHKlStbi+T0yFF70v2XBg4cKPoKigQBPxBAHPwQRXyIGoHExETZvn27xMfHG3Vu2rRJHnnkEcOOAQJeJIA4eDFqtDmmBFq3bi0bN26UUqVKGe146623ZOrUqYYdAwS8RoABaa9FjPa6hsC6deukf//+ouse7ElFQqfAkiDgVQL0HLwaOdodcwJ6SNDLL7/s2A49YU73aCJBwKsEEAevRo52u4LApEmTZNq0aUZbdGB6wIAB1klyRiYGCHiAAK+VPBAkmuh+AiNGjJClS5caDdUZTnrCnK6FIEHASwQQBy9Fi7a6loAeBNSjRw9rJpO9kTrDKSMjQ3Q1NQkCXiEQFXH4j//4D8dfmvxC6tKlizRq1OiGxYsUKWLkX7t2zbBhgECkCJw9e1Y6d+4sR48eNarQ/Zd0Hybdj4kEAS8QiIo4bNiwQXr37l1gHqtXr7ZmhdzoAYjDjeiQFy0Cf/7zn+X222+Xr776yqiyQ4cOsmvXLscpsEZhDBCIMYGoDEjrP9wlSpTI93/FixePMRaqh0DBCCQkJFirqGvUqGE8IC0tTQYPHuw49dUojAECMSYQlZ5DqD6+9NJL8uSTT1q36esk7abn1R2n5xAqZcpHksAnn3wiHTt2lNOnTxvVTJw4UebPn2/YMUDATQRcJw579+6V++67T3SAr3z58nL48GFp0aJFnswQhzwRUSDKBN5//3158MEH5dKlS0bNOv11+vTphh0DBNxCICqvlfLr7DfffCMpKSmWMOg9S5YsyZcw5Pf5lINANAnoFt/6HXb6w+Vf//VfZfHixdFsDnVBICQCrhEH3e5YD1TRAT1NU6ZMkT59+oTkDIUh4DYCulPr888/79issWPHyrZt2xzzMEIg1gRc81ppyJAhsmLFCovH/fffL1u3bpWiRfOvXU5/nTGVNdZfL+rPIvDMM8/Is88+m3WZ/bNs2bLywQcfyI9+9KNsGx8g4AYCrhCHF198UX76059aPBo2bCgfffRRngPQdniIg50I124ioH+oDB06NPsPoJxt0xlOeh5EkyZNcpr5DIGYEoi5OOignS5y071oypUrZw1AF+TAdsQhpt8jKs8HgcuXL0vXrl1l9+7dRul69epZR5A6TYE1CmOAQBQI5P+9TQQa87vf/c5a3JZ1epYO3hVEGCLQNB4JgbAT0LU+69evl3bt2hnP1t+FBx54QL7//nsjDwMEYkEgZuJw7tw5awD6L3/5i+W3DkDrFsgkCPiZgE7P1kHoBg0aGG4eP37c+p1wmvpqFMYAgQgTiJk46EEo+sugSQegZ8yYEWFXeTwE3EGgevXq1irqqlWrGg3Swelhw4YJkykMNBiiTCAm4qBT+1auXGm5qgPQ+jmUmUlRZkR1EAg7AV35rz0I7UnY05o1a7J3CLDncQ2BaBGI+oD0nj17rJ5CYQeg7YAYkLYT4doLBHbu3Cndu3cXHay2J+1Ncx61nQrX0SIQ1Z7D119/bZ2OlTUArStEGYCOVqipx40EdKbewoULHVdR/8u//IssX77cjc2mTQEgEDVxsA9A//znP5d+/foFADEuQuDGBHT9w3PPPWcU0nGHUaNGWdt8G5kYIBBhAlF7raTbCOi5DJr0ryV93xrOcQZeK0X4m8LjI05g8uTJMm/ePKMeHZfYt2+ftG/f3sjDAIFIEYiKOMydO1eeeuopywcdgNYtuCtVqhRWnxCHsOLkYTEgoD0F/SNKB6TtqVq1atZRo/r7Q4JANAhEXBxUCH784x9bK6DVoUmTJkmrVq1C8i0pKcnqbdzoJsThRnTI8woBXeOgU7t163p7UmHQs6hVKEgQiDSBiIuD/hU0YMCAQvnRrVs3SU1NveEzEIcb4iHTQwR0lbQeKfrpp58ardZXS/qKyWkKrFEYAwQKQSBqA9KFaCO3QiBQBOLi4mT79u2i+y3Z07Fjx6yt7J2mvtrLcg2BwhCIeM+hMI0L5V56DqHQoqwXCJw8eVLuuOMOydpiJmebdYv7ZcuWOU6BzVmOzxAoKAF6DgUlx30QiDCB5ORk61wTPfPBnnT9Q9YkD3se1xAIBwHEIRwUeQYEIkRADwFau3atFC9e3KhBp73OmTPHsGOAQDgI8FopHBR5BgQiTGDRokUyevRooxZ9nap7kxV20ofxYAyBJ0DPIfBfAQB4gYCulJ49e7bRVF0bMXz4cNFDs0gQCCcBeg7hpMmzIBBhAhMnTpTXXnvNqKVChQqSlpYmN910k5GHAQIFIYA4FIQa90AgRgSuXr1q7UmmJ8rZkx4xqovk6tevb8/iGgIhE0AcQkbGDRCILYGLFy/KfffdZ/UU7C1p0qSJHDx4UBISEuxZXEMgJAKMOYSEi8IQiD2BUqVKyaZNmxy3uz916pR07dpVzp49G/uG0gJPE0AcPB0+Gh9UAhUrVpQdO3ZI3bp1DQRHjhyRlJQUyczMNPIwQCC/BBCH/JKiHARcRqB27drWWdROOxzrlvh6TjsJAgUlgDgUlBz3QcAFBJo1a2ZtSlmmTBmjNUuWLGEVtUEFQ34JIA75JUU5CLiUgO6/tGrVKilWrJjRQj1LZf78+YYdAwTyIoA45EWIfAh4gEDPnj0d1z9o0x9//HFZt26dB7ygiW4igDi4KRq0BQKFIDB27FiZMWOG8QRdG6G7uO7fv9/IwwCB3AiwziE3Mtgh4FECKhILFiwwWh8fHy/p6ekhn8RoPAhDIAggDoEIM04GicCVK1ekd+/esnnzZsPtWrVqyaFDhyQxMdHIwwCBnAQQh5w0+AwBnxA4f/683HPPPdZ2GnaXmjZtaq2irly5sj2LawhkE2DMIRsFHyDgHwI6tVXPXdeprvb05Zdfip7LrgJCgkBuBBCH3Mhgh4DHCWjPQFdR62I5e9JXS3oGhL6CIkHAiQDi4EQFGwR8QkDHFlQgdLsNe9IxifHjx9vNXEPAIoA48EWAgM8JtGzZ0tqoTzfss6e3335bpk6dajdzDQFBHPgSQCAABDp06CArVqyQokXNX/lZs2bJm2++GQAKuBgKAfObEsrdlIUABDxDoE+fPvLKK684tldPmNu4caNjHsZgEkAcghl3vA4oARWBadOmGd7rwPTAgQOtKa5GJoZAEmCdQyDDjtNBJzBq1ChZvHixgUG3/z5w4IA0b97cyMMQLAKIQ7DijbcQsAjoQUC6WZ+e+2BPeoCQnkVdp04dexbXASKAOAQo2LgKgZwEzp07J3fddZd8+OGHOc3W5xYtWlg9CKcpsEZhDL4kwJiDL8OKUxDIm0DZsmVly5Yt0qRJE6PwZ599Jj169JALFy4YeRiCQQBxCEac8RICjgQSEhKso0Zr1Khh5OsOroMHDxbd8psUPAKIQ/BijscQuI5AvXr1ZPv27RIXF3edXS/Wr18vkyZNMuwY/E8AcfB/jPEQAnkSuOmmm2TDhg1SsmRJo+zrr78u06dPN+wY/E2AAWl/xxfvIBASgTVr1ljrHa5du2bct3DhQtEpsKRgEKDnEIw44yUE8kWgf//+8sILLziW1RPmdACbFAwCiEMw4oyXEMg3gSeeeEKmTJlilNdV1Coehw8fNvIw+I8Ar5X8F1M8gkChCehrpWHDhsny5cuNZ1WpUsXaZiM5OdnIw+AfAoiDf2KJJxAIK4HLly9bJ8bt2rXLeG5SUpJ1FnXNmjWNPAz+IIA4+COOeAGBiBA4c+aMdOrUSY4dO2Y8v3Xr1qJrIZymwBqFMXiOAGMOngsZDYZA9AiUL1/e2n+pYcOGRqW//vWvpVevXnLp0iUjD4P3CSAO3o8hHkAgogSqVatmraLWn/a0d+9eGTp0KKuo7WB8cI04+CCIuACBSBPQnoPu4Ko9CXtau3atPPnkk3Yz1x4ngDh4PIA0HwLRItC+fXt57733pESJEkaVesLczJkzDTsG7xJgQNq7saPlEIgJAZ3eqtNcnVZRL1261MqLScOoNKwE6DmEFScPg4D/CQwZMkTmzJnj6Ojo0aNlx44djnkYvUWAnoO34kVrIeAaApMnT5Z58+YZ7SlXrpzoQPUtt9xi5GHwDgHEwTuxoqUQcBUBfa00aNAgWb16tdGuqlWrWkeNNmrUyMjD4A0CiIM34kQrIeBKArrG4cEHH5T333/faF+DBg0sgahevbqRh8H9BBAH98eIFkLA1QROnz4tHTp0kOPHjxvtbNeunezbt08qVKhg5GFwNwEGpN0dH1oHAdcT0H/49SS5+vXrG239+OOPpXfv3qL7NJG8RQBx8Fa8aC0EXElAz6DWWUp6JrU97dmzR0aMGOE49dVelmv3EEAc3BMLWgIBTxNo0qSJbN26VXS2kj2tXLlSdHYTyTsEEAfvxIqWQsD1BG699VbR7TSKFy9utPX555+XuXPnGnYM7iTAgLQ740KrIOBpAkuWLJGRI0caPhQpUsQ6QEinwJLcTYCeg7vjQ+sg4EkCOsYwe/Zso+26NkLzdByC5G4C9BzcHR9aBwFPE3j00Udl/vz5hg86w2n//v3Stm1bIw+DOwggDu6IA62AgC8JXL16Vfr37y/r1q0z/NPFcRkZGaKL5UjuI4A4uC8mtAgCviJw8eJF6dKli9VTsDvWuHFjOXjwoOh2GyR3EUAc3BUPWgMBXxL47rvv5Cc/+YmcOHHC8E836NON+pymwBqFMUSNAAPSUUNNRRAILoH4+HhrkVxiYqIB4ejRo9K3b1/JzMw08jDEjgDiEDv21AyBQBGoVauWdRZ15cqVDb91dbWeBUFyDwHEwT2xoCUQ8D2Bpk2bSmpqqpQpU8bwddmyZfLUU08ZdgyxIYA4xIY7tUIgsARuv/126wyIYsWKGQx0BfWvfvUrw44h+gQQh+gzp0YIBJ5Ajx495I033nDk8MQTT8i7777rmIcxegQQh+ixpiYIQCAHgYcfflhmzJiRw/KPj7o2YujQodY5EEYmhqgRYCpr1FBTEQQg4ERg/Pjx8uabbxpZOsMpLS1NWrdubeRhiDwBxCHyjKkBAhC4AYErV65YU1k3btxolNIZTrqKOikpycjDEFkCiENk+fJ0CEAgHwQuXLgg99xzj7Va2l48OTnZslepUsWexXUECTDmEEG4PBoCEMgfgdKlS1tTXJs3b27ccPLkSenWrZucO3fOyMMQOQKIQ+TY8mQIQCAEApUqVbJWUdepU8e46/Dhw9YGfvoKihQdAohDdDhTCwQgkA8CdevWtQSiYsWKRuktW7bIuHHjDDuGyBBAHCLDladCAAIFJNCiRQvZvHmz6Ksme1q4cKFMnTrVbuY6AgQQhwhA5ZEQgEDhCOgOritWrJCiRc1/ombNmiWvv/564Srg7jwJmOTzvIUCEIAABCJPoHfv3o6nyGnNkyZNkg0bNkS+EQGuAXEIcPBxHQJuJzBhwgT55S9/aTRTV1EPGjRI0tPTjTwM4SHAOofwcOQpEIBABAnodt6LFi0yatCB6wMHDoiOU5DCSwBxCC9PngYBCESAgE5h7dWrl+iMJXvSqa+6ilpnOpHCRwBxCB9LngQBCESQgC6Cu/vuu0XXPNiTLp7THoSulSCFhwBjDuHhyFMgAIEIEyhbtqzVc9DtNOzp888/F90GXLfhIIWHAOIQHo48BQIQiAIB3V9p586dUrNmTaM27TkMHDhQWEVtoCmQAXEoEDZuggAEYkVAd2jVM6fj4uKMJujOrhMnTjTsGEIngDiEzow7IACBGBPQMx5UCEqVKmW0RM+GmD59umHHEBoBBqRD40VpCEDARQTWrl1rvUrSdQ/2tGDBAtHT5kgFI0DPoWDcuAsCEHABgZSUFHnxxRcdW6InzKWmpjrmYcybAOKQNyNKQAACLibw2GOPyZQpU4wW6sB0//795dChQ0YehrwJ8Fopb0aUgAAEPEBg+PDhsmzZMqOllStXtk6Sa9q0qZGHIXcCiEPubMiBAAQ8RCAzM1O6d+9uzWSyNzsxMdHqQeiZ1KT8EUAc8seJUhCAgAcInD17Vjp37ixHjx41WtuqVStro774+HgjD4NJgDEHkwkWCEDAowTKlSsnW7dulUaNGhkenDhxQnr27CkXL1408jCYBBAHkwkWCEDAwwSqVq1qraKuXr264cX+/ftlyJAh4jT11SgccAPiEPAvAO5DwI8EGjRoINu2bZMKFSoY7q1bt04ef/xxw47hegKIw/U8uIIABHxCoF27drJ+/XopUaKE4dH8+fNl5syZhh3DDwQYkP6BBZ8gAAEfEli5cqX1KunatWuGd4sXL5YRI0YYdgwi9Bz4FkAAAr4moMeJzp0719HHMWPGyPbt2x3zgm6k5xD0bwD+QyAgBCZPnizz5s0zvNUZTh988IHceuutRl6QDYhDkKOP7xAIEAF9raQzlfQ1kz0lJCRYR402btzYnhXYa8QhsKHHcQgEj8Dly5flwQcflD179hjO169f3xKIGjVqGHlBNCAOQYw6PkMgwAROnz4tHTt2lE8++cSg0LZtW9G1EE5TYI3CPjcwIO3zAOMeBCBwPQH9h18HoXUthD2pYDz00ENy6dIle1bgrhGHwIUchyEAAV09rUeN6mpqe3r//fet6a1OU1/tZf18jTj4Obr4BgEI5EpAB591HyadrWRPq1atkp/97Gd2c6CuEYdAhRtnIQCBnARuueUW0e00ihcvntNsfdYT5nJbH2EU9qGBAWkfBhWXIACB0AjoIUF6WJA9FSlSRN555x0ZPHiwPcv31/QcfB9iHIQABPIiMGzYMJk9e7ZRTMcdRo4cKbt37zby/G6g5+D3COMfBCCQbwJ6HvWvfvUro3z58uWtKa66mV9QEj0HW6S1G7lv3z6blUsIQCAIBF566SXp16+f4eqZM2fkgQcekN/+9rdGnl8NiINfI4tfEAgQgXBNOy1atKg1xtCpUyeD3v/8z//I/fffL/ozCAlxCEKU8RECPiTwn//5nzJ69GjRs6F1y4vbb79djh07VmhPS5UqJRs3bpTWrVsbz/rNb34jXbt2Fe1J+D0hDn6PMP5BwIcEdB1CixYt5Ouvv7ZeA3/zzTfyt7/9TfQv/i+//LLQHsfHx1urqJOSkoxnffTRR9K3b1/RfZr8nBAHP0cX3yDgQwKpqakydOhQa/8jXaNQpUoVKV26tAwYMMD6i/6FF14Ii9e1atWyzqLW59vTzp07rV5LuF5n2Z/vhmvEwQ1RoA0QgEC+CGivQEXg6tWr1thA5cqVs+/L2m579erV2bbCfkhOTpYtW7ZI2bJljUfp+oenn37asPvFgDj4JZL4AYEAEJg0aZKcO3dOevbsKZ07d77O46x/wHXXVX3FFK502223yZo1a6RYsWLGI3UF9csvv2zY/WBAHPwQRXyAQAAIbNiwQXRTPE16qps9/d///V+2SQerw5m6desmb775puMjn3zySUs8HDM9bEQcPBw8mg6BIBF46623LHd//OMfWzOT7L5/++232Sanv/KzMwv4QWdGzZgxw7hbxx10hbUeNeqnhDj4KZr4AgGfEtB/+LO2sBg4cKCjlzlnKVWrVs2xTGGNU6dOlQkTJhiP0fMf9ByITz/91MjzqgFx8GrkaDcEAkRg+fLlcuXKFcvj7t27O3qedbKbTkN1mmHkeFMBjPPnz7eEwH7r999/b62i/t3vfmfP8uS1uU+tJ92g0RCAgJ8JLF261HJPewQqAllCkOWzzl764osvrEsdQNZtcCKVdBX1ypUr5b777pP09PTrqtEeTpcuXeTgwYOSkJBwXZ7XLhAHr0WM9kIgYAQ+++wz0f806QZ4es6CPemWFpmZmZa5Q4cO9uywX+u6is2bN8udd96Z3basSk6dOiU6gK1jEFkzqLLyvPST10peihZthUAACRw5ciTbaxWGtLQ047+szfL00B6ncxmyHxDGDxUrVrSOGq1bt67x1A8//FBSUlKyBcso4AED4uCBINFECASZwMcff5ztvs5Uckp6HrSmXr16Se3atZ2KRMRWp04dSyAqVapkPF+PIB07dqxh94oBcfBKpGgnBAJK4A9/+IPlec2aNaVq1aoGha+++kp0vyNNEydONPIjbWjevLn1iklfNdnT4sWLRWc4eTEhDl6MGm2GQIAI/PnPf7a8zW16qm5joUk34tON92KRdOxBNwN0Wl8xa9Ysee2112LRrELViTgUCh83QwACkSZQsmRJqwqnVze6lUbWyuVf/OIXkW7KDZ+vr7ReffVVxzKPPvqovPfee455bjUiDm6NDO2CAAQsAro7qianQ3b0L3K1615LuiFfrNO4cePkl7/8pdEMnWo7ePBgayDdyHSpgTOkbYEJ9/xoP2/pa0PHJQQiQkA3t3vqqadED+HRnoKuM9CkZzm0b99etEeRkZEh1atXj0j9BXnomDFj5O233zZu1RlOujaiZcuWRp7bDIiDLSKIgw0IlxCIMYG//vWvotNFVRgOHDggd9xxh+gme7orq45H6NTW+vXrx7iV11evq7l1Ow09e8KedDbVoUOHLJ/seW665rWSLRr6l344/7M9nksIQCBEAnpmg56poAPSffr0kSFDhoie0KZ/fR8/ftx1wqDu6cC0bvPtNPX2j3/8o7WKWkXPzYmeg5ujQ9sgAIFsArq5nW6R8b//+7/Zr5OyM136QQVAezo5NwXMaqradTPBMmXKZJlc9RNxcFU4aAwEIOA3Ar///e+tHsSf/vQnwzUdSNdZTE5TYI3CUTbwWinKwKkOAhAIFoHExERrFbXuFmtPmzZtkkceecRudsU14uCKMNAICEDAzwRatWolKgQ648qe9BCj6dOn280xv+a1UsxDQAMgAIGgEFi3bp30799fdN2DPeliPjftxUTPwR4hriEAAQhEiEDfvn3l5Zdfdny6vl7S3oVbEuLglkjQDghAIBAEJk2aJFOmTDF81bURegSqLuhzQ+K1khuiQBsgAIHAERg5cqQsWbLE8FvXdehiv2bNmhl50TQgDtGkTV0QgAAE/klAT67r0aOHbN++3WCiM5y0BxHNsynsjUAc7ES4hgAEIBAlAmfPnpW77rpLcp52l1W1rgDXfZh0P6ZYJF+LQyyAUicEIACBcBHQ87B37drlOAU2XHXk9hwGpHMjgx0CEIBAjAnopoK61bfT1NdINw1xiDRhng8BCECgEAR0ew09LCjaCXGINnHqgwAEIBAiAT3UaObMmSHeVbjijDkUjh93QwACEIgagUWLFolOgY1GKh6NSqJRByeuRYMydUAAAkEhwGuloEQaPyEAAQiEQABxCAEWRSEAAQgEhQDiEJRI4ycEIACBEAggDiHAoigEIACBoBBAHIISafyEAAQgEAIBxCEEWBSFAAQgEBQCiENQIo2fEIAABEIggDiEAIuiEIAABIJCAHEISqTxEwIQgEAIBBCHEGBRFAIQgEBQCCAOQYk0fkIAAhAIgQDiEAIsikIAAhAICgHEISiRxk8IQAACIRBAHEKARVEIQAACQSGAOAQl0vgJAQhAIAQCiEMIsCgKAQhAICgEEIegRBo/IQABCIRAwDcnwYXgM0UhkCeB9evXy7fffis1atSQPn365FmeAhDwGwHfnCHtt8DgT2wJ3HHHHZKRkSG33XabHDp0KLaNoXYIxIAAr5ViAJ0qIQABCLidAOLg9gjRPghAAAIxIIA4xAA6VUIAAhBwOwHEwe0Ron0QgAAEYkCA2UoxgE6VsSVw+fJl2bJlixw9etT6709/+pMkJiZKly5dZPTo0VKhQoXYNpDaIeACAoiDC4JAE6JHQKen9uvXTw4ePHhdpZ9//rns2LFDnnvuOVm7du11eVxAIIgEEIcgRj2gPp86dUo6duwo//Vf/2URaNq0qXWtP0+cOCHp6eny1VdfSY8ePeg9BPQ7gts/EEAcfmDBJ58TeOaZZ7KF4ZFHHpH58+dLkSJFsr2+cuWKPP744/Lqq6/Kd999l23nAwSCSIAB6SBGPYA+f/LJJ6KrnjU98MADlgDkFAa1FytWzBKMvn376iUJAoEmgDgEOvzBcX7lypVy7do1y+FHH330ho4//fTTN8wnEwJBIIA4BCHK+Ci//e1vLQolSpSwZiXdCMnNN98stWrVulER8iDgewKIg+9DjINK4De/+Y0FQv/RL1o07699vXr1rPL8DwJBJZD3b0lQyeC3rwicOXPG8qdixYr58isuLi5f5SgEAb8SQBz8Gln8uo5Aw4YNrevf//7319lzu/jDH/6QWxZ2CASCAOIQiDDjZOPGjS0If/vb3+SPf/zjDYFoL+Prr7++YRkyIeB3AoiD3yOMfxYBPZ8hK73yyitZHx1/LlmyRM6dO+eYhxECQSGAOAQl0gH3c8CAAZKcnGxReP311+XIkSOORHSF9LPPPuuYhxECQSKAOAQp2gH2VRe4/du//ZtF4PTp09KpUyd57bXXsl8f6eZ7uhZCexj6OT8zmgKME9cDQIDtMwIQZFz8B4GUlBT561//am2Rcf78eZk4caKVobuwqmBkpZ///OeSlpZmHROaZeMnBIJGgJ5D0CIecH/Hjx8ve/fulQ4dOki5cuUsGlnCULt2bZk7dy6vlQL+HcH9fxAo8vctBf6xpwBEIBAwArrRnm7V/d///d/Spk0bqVq1asAI4C4EcieAOOTOhhwIQAACgSXAa6XAhh7HIQABCOROAHHInQ05EIAABAJLAHEIbOhxHAIQgEDuBBCH3NmQAwEIQCCwBBCHwIYexyEAAQjkTgBxyJ0NORCAAAQCSwBxCGzocRwCEIBA7gQQh9zZkAMBCEAgsAQQh8CGHschAAEI5E4AccidDTkQgAAEAksAcQhs6HEcAhCAQO4EEIfc2ZADAQhAILAEEIfAhh7HIQABCOROAHHInQ05EIAABAJLAHEIbOhxHAIQgEDuBP4fq4WZZhU60VQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 9,
     "metadata": {
      "image/png": {
       "width": 300
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(\"Image1.png\", width=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique, counts = np.unique(em_data['GpsTime'], return_counts=True)\n",
    "beam_indexes = np.array([[i,v] for i,v in enumerate(list(counts))])\n",
    "\n",
    "ping_indexes = np.cumsum(beam_indexes[:,1])\n",
    "input_indexes = list(range(len(ping_indexes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#x, y, z = em_data['X'], em_data['Y'], em_data['Z']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit, prange\n",
    "\n",
    "\n",
    "@jit(nopython=True, nogil=True)\n",
    "def getAngle(ping_index, x, y, z, ping_indexes, beam_indexes):\n",
    "    if ping_index == 0:\n",
    "        start=0\n",
    "        stop=ping_indexes[ping_index]\n",
    "        shape=ping_indexes[ping_index]\n",
    "        nadir_index = int( beam_indexes[ping_index][1]/2)\n",
    "    else:\n",
    "        start=ping_indexes[ping_index-1]\n",
    "        stop=ping_indexes[ping_index]\n",
    "        shape=ping_indexes[ping_index] - ping_indexes[ping_index-1]\n",
    "        nadir_index = int(np.sum(beam_indexes[0: beam_indexes[ping_index][0]][:,1]) + \n",
    "                          ( beam_indexes[ping_index][1]/2))\n",
    "    XY = np.array([(x[nadir_index], y[nadir_index])])\n",
    "    pings = np.column_stack((x[start:stop], y[start:stop])) \n",
    "    Zi = z[start:stop] \n",
    "    a_min_b = XY - pings\n",
    "    d = np.sqrt(np.sum((a_min_b)**2, axis=1))\n",
    "    angle = np.rad2deg(np.arctan2(np.abs(Zi),d)).reshape(shape,)\n",
    "    #return [start, stop, angle]\n",
    "    return angle"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# usage\n",
    "getAngle(57, XUTM, YUTM, em_data['Z'], ping_indexes, beam_indexes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4320"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import multiprocessing\n",
    "ncpu = multiprocessing.cpu_count()\n",
    "chunksize = unique.shape[0]/ncpu\n",
    "chunksize = int(chunksize - (chunksize%48))\n",
    "chunksize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "211680it [00:39, 6880.09it/s]                            \n"
     ]
    }
   ],
   "source": [
    "angles = parmap.map(getAngle, \n",
    "                    input_indexes, \n",
    "                    em_data['X'], em_data['Y'], em_data['Z'], \n",
    "                    ping_indexes, \n",
    "                    beam_indexes, \n",
    "                    pm_processes=48, \n",
    "                    pm_pbar=True, \n",
    "                    pm_parallel=True, \n",
    "                    pm_chunksize=chunksize)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_values = [beam for ping in angles for beam in ping]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pings = np.column_stack((em_data['X'], \n",
    "                         em_data['Y'], \n",
    "                         em_data['Z'], \n",
    "                         em_data['Amplitude'], \n",
    "                         angle_values))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "end_time = dt.datetime.now()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "datetime.timedelta(0, 48, 303721)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "end_time-st_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://gist.github.com/59fd3dbaa7e10fa449bdf11ff4366bcb\n"
     ]
    }
   ],
   "source": [
    "!gist EM1000.ipynb mb_reader.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

mb_reader.py

import glob
from string import Template
import pandas as pd
import pdal
import sys

def gen_input(directory, 
              file_extension='.ALL', 
              reader_driver='mbio', 
              file_format='MBF_EMOLDRAW', 
              output_type='pandas.DataFrame',
              reproject=False,
              in_srs=None,
              out_srs=None,
              verbose=False):
    assert output_type in ['numpy.array', 'pandas.DataFrame'], "Wrong output type"
    dirlist= "{directory}/*{file_extension}".format(directory=directory, 
                                                    file_extension=file_extension)
        
    file_list = glob.glob(dirlist)
    for result in file_list:
        if reproject and in_srs and out_srs:
            yield {"file_name": result, 
                   "reader_driver": reader_driver, 
                   "file_format": file_format, 
                   "output_type": output_type, 
                   "verbose": verbose,
                   "reproject": reproject,
                   "in_srs": in_srs,
                   "out_srs": out_srs}
        else:
            yield {"file_name":result, 
                   "reader_driver":reader_driver, 
                   "file_format":file_format, 
                   "output_type":output_type, 
                   "verbose":verbose}

def readEM1000_1(args):
    file_name = args['file_name'] 
    if 'reader_driver' in args:
        reader_driver=args['reader_driver']
    else:
        reader_driver='mbio'        
    if 'file_format' in args:
        file_format=args['file_format']
    else:
        file_format='MBF_EMOLDRAW'
    if 'output_type' in args:
        output_type=args['output_type']
    else:
        output_type='numpy.array'
    assert output_type in ['numpy.array', 'pandas.DataFrame'], "rong output type"
    if 'verbose' in args:
        verbose=args['verbose']
    else:
        verbose=False
    t = Template('{"pipeline":[{"filename": "${file_name}","type":"readers.${reader_driver}","format" : "${file_format}"}]}')
    json = t.substitute(file_name=file_name, reader_driver=reader_driver, file_format=file_format)
    p = pdal.Pipeline(json)
    p.validate() # check if our JSON and options were good
    p.loglevel = 8 #really noisy
    count = p.execute()
    data = p.arrays[0]
    if verbose:
        if verbose == 1:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
        if verbose == 2:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
            print('Dimension names are', data.dtype.names)
        if verbose == 3:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
            print('Dimension names are', data.dtype.names)
            print('Metadata: ', p.metadata)
            print('Log: ', p.log)
    if output_type == 'numpy.array':
        return data
    if output_type == 'pandas.DataFrame':
        return pd.DataFrame(data)
    if output_type == 'count':
        return count

def readEM1000(args):
    file_name = args['file_name'] 
    if 'reader_driver' in args:
        reader_driver=args['reader_driver']
    else:
        reader_driver='mbio'        
    if 'file_format' in args:
        file_format=args['file_format']
    else:
        file_format='MBF_EMOLDRAW'
    if 'output_type' in args:
        output_type=args['output_type']
    else:
        output_type='numpy.array'
    assert output_type in ['numpy.array', 'pandas.DataFrame'], "wrong output type"
    if 'verbose' in args:
        verbose=args['verbose']
    else:
        verbose=False
    if all(opt in args for opt in ['reproject', 'in_srs', 'out_srs']):
    #if args['reproject']:
        in_srs=args['in_srs']
        out_srs=args['out_srs']
        t = Template('{"pipeline":[{"filename": "${file_name}","type":"readers.${reader_driver}","format" : "${file_format}"}, {"type":"filters.reprojection", "in_srs":"${in_srs}", "out_srs":"${out_srs}"}]}')
        json = t.substitute(file_name=file_name, 
                            reader_driver=reader_driver, 
                            file_format=file_format, 
                            in_srs=in_srs, 
                            out_srs=out_srs)
        

    else:
        t = Template('{"pipeline":[{"filename": "${file_name}","type":"readers.${reader_driver}","format" : "${file_format}"}]}')
        json = t.substitute(file_name=file_name, reader_driver=reader_driver, file_format=file_format) 
    #print(json)
    p = pdal.Pipeline(json)
    p.validate() # check if our JSON and options were good
    p.loglevel = 8 #really noisy
    count = p.execute()
    data = p.arrays[0]
    if verbose:
        if verbose == 1:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
        if verbose == 2:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
            print('Dimension names are', data.dtype.names)
        if verbose == 3:
            print('Read', count, 'points with', len(data.dtype), 'dimensions')
            print('Dimension names are', data.dtype.names)
            print('Metadata: ', p.metadata)
            print('Log: ', p.log)
    if output_type == 'numpy.array':
        return data
    if output_type == 'pandas.DataFrame':
        return pd.DataFrame(data)
    if output_type == 'count':
        return count