src/notebooks/save_nhd_extract.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c8e9b643",
   "metadata": {},
   "source": [
    "# Save a file for each HUC with reach geoms and ids\n",
    "\n",
    "To make things more cloud-friendly, we want to save a portion of the NHD database in a set of GeoJSON files which can be stored on S3. Each file will have the reach geometries and associated ids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "45c95b9c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR 1: PROJ: proj_create_from_database: Open of /srv/conda/envs/notebook/share/proj failed\n"
     ]
    }
   ],
   "source": [
    "import json\n",
    "import os\n",
    "import collections\n",
    "\n",
    "import psycopg2\n",
    "import shapely\n",
    "import shapely.wkt\n",
    "import geopandas as gpd\n",
    "import xarray as xr\n",
    "import fsspec\n",
    "import numpy as np\n",
    "import pyproj\n",
    "import dask.bag as db\n",
    "from dask.distributed import Client\n",
    "import logging\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ddf90be",
   "metadata": {},
   "source": [
    "To be able to generate a separate geojson for each of the huc8 units (around 2300) we need the following tables: `wbdhu8` (in database `nhdplushr`) and `nhdflowline` (in database `nhdplusv2`). This is not _the only_ way to do this...it's just that this is how we went about it. To generate each of these tables download all the necessary files pertaining to NHDPlus V2 and NHDPlus HR and follow the instructions of creating the database listed  in the [project's README.md](https://github.com/azavea/noaa-hydro-data#local-database-setup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "320dbbf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cursor(database):\n",
    "    connection = psycopg2.connect(host='database', database=database,user='postgres', password='password')\n",
    "    cursor = connection.cursor()\n",
    "    return cursor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc3e3138",
   "metadata": {},
   "source": [
    "# Get all HUC8s and write a GeoJSON file for each one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c2506429",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_huc8s(database):\n",
    "    cursor = get_cursor(database)\n",
    "    query = 'select huc8 from wbdhu8'\n",
    "    cursor.execute(query)\n",
    "    return [c[0] for c in cursor]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6b69f0a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_huc8_extract(huc8, out_dir):\n",
    "    # get huc8 boundary\n",
    "    cursor = get_cursor('nhdplushr')\n",
    "    query = 'SELECT wkb_geometry from wbdhu8 WHERE huc8=%s'\n",
    "    cursor.execute(query, [huc8])\n",
    "    ## huc_geom = shapely.wkb.loads(cursor.fetchone()[0].tobytes())\n",
    "    huc_geom = shapely.wkb.loads(cursor.fetchone()[0], hex=True)\n",
    "\n",
    "    # get reaches intersecting with huc boundary\n",
    "    cursor = get_cursor('nhdplusv2')\n",
    "    query = f'''\n",
    "        SELECT comid, ST_Force2D(wkb_geometry) from nhdflowline WHERE ST_Intersects(\n",
    "            ST_Transform(ST_GeomFromEWKB(wkb_geometry), 4326), ST_GeomFromGeoJSON(%s))\n",
    "        '''\n",
    "    huc_geom_str = json.dumps(shapely.geometry.mapping(huc_geom))\n",
    "    cursor.execute(query, [huc_geom_str])\n",
    "    reach_geoms = []\n",
    "    reach_ids = []\n",
    "    for reach_id, reach_geom in cursor:\n",
    "        reach_ids.append(int(reach_id))\n",
    "        reach_geoms.append(shapely.wkb.loads(reach_geom, hex=True))\n",
    "\n",
    "    # make dataframe with comid and geometries and save to GeoJSON\n",
    "    df = gpd.GeoDataFrame({'comid': reach_ids + [0], 'geometry': reach_geoms + [huc_geom]})\n",
    "    out_path = os.path.join(out_dir, f'{huc8}.json')\n",
    "    df.to_file(out_path, driver='GeoJSON', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26dd46f9",
   "metadata": {},
   "source": [
    "Having defined the necessary extract utility functions let's try to generate the json files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d95c104c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR 1: PROJ: proj_create_from_database: Open of /srv/conda/envs/notebook/share/proj failed\n"
     ]
    }
   ],
   "source": [
    "limit = npartitions = 10\n",
    "client = Client()\n",
    "\n",
    "huc8s = get_huc8s('nhdplushr')\n",
    "huc8_out_dir = f'/opt/data/noaa/huc8-extracts/'\n",
    "os.makedirs(huc8_out_dir, exist_ok=True)\n",
    "huc8s_processed = [f.replace('.json','') for f in os.listdir(huc8_out_dir)]\n",
    "if len(huc8s_processed) < len(huc8s):\n",
    "    huc8s_missing = list(set(huc8s) - set(huc8s_processed))\n",
    "    huc8_bag = db.from_sequence(huc8s_missing, npartitions=npartitions)\n",
    "    out8 = huc8_bag.map(save_huc8_extract, huc8_out_dir).compute()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "655e0bcd",
   "metadata": {},
   "source": [
    "The first time we ran the notebook we generated a geojson for each of the huc8 which contained not only the `comid` of each reach but also it's associated geometry (linestring). However, for our experiments we need only the boundary of huc8 and the associated `comid`s. Therefore, we transform the geojson to extract out the portions that we need for our experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1348c3bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform_geojsons(json_dir='/opt/data/noaa/huc8-extracts', out_dir='/opt/data/noaa/huc8-extracts/transformed'):\n",
    "    \n",
    "    logger = logging.getLogger('distributed.utils_perf')\n",
    "    logger.setLevel(logging.ERROR)\n",
    "    \n",
    "    jsons = [f for f in os.listdir(json_dir) if f.endswith('.json')]\n",
    "\n",
    "    for jsonfile in jsons:\n",
    "        huc8 = jsonfile.replace('.json', '')\n",
    "        out_path = f'{out_dir}/{huc8}.json'\n",
    "        if os.path.exists(out_path):\n",
    "            continue\n",
    "        dat = gpd.GeoDataFrame.from_file(f'{json_dir}/{jsonfile}', driver='GeoJSON')\n",
    "        comids = [c for c in dat.loc[:,'comid'][:-1]]\n",
    "        huc_geom = dat.iloc[-1,:]['geometry']\n",
    "        gdf = gpd.GeoDataFrame(data = {'huc8': [huc8], 'comids': [comids]}, geometry=[huc_geom], crs='EPSG:4326')\n",
    "        with open(out_path, 'w') as fw:\n",
    "            fw.write(gdf.to_json())\n",
    "\n",
    "\n",
    "transform_geojsons()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "488794e1",
   "metadata": {},
   "source": [
    "**NOTE:** We have to use `gdf.to_json()` instead of `gdf.to_file('02040202.json', driver='GeoJSON')` because neither geopandas nor shapely can read lists from geojsons. Geojsons can contain any json object ([geojson ref](https://datatracker.ietf.org/doc/html/rfc7946#section-3.2)) as properties but geopandas ([see here](https://github.com/geopandas/geopandas/issues/2113)) and fiona ([here](https://github.com/Toblerity/Fiona/issues/738)) are unable to handle list (or array in json)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2f39c4d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index(['id', 'huc8', 'geometry'], dtype='object')\n",
      "# of features => 1\n",
      "huc8 => 02040202\n",
      "first 5 comids => [4491350, 4495894, 4496294, 4491376, 4496504]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADoCAYAAADheS9sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABD+ElEQVR4nO3deVhV1f7H8fcCcQJRFAUScR5QwQFS1JzKIeeprjmWNsnVNMtS62bZtTQzTXEoLa+aU+acWmqp4ZSKiChOOOSIgqKAiIzf3x+SP01AhHPYh8N6Pc95OMPee32Wcr5nsc/aeysRQdM0TbNeNkYH0DRN08xLF3pN0zQrpwu9pmmaldOFXtM0zcrpQq9pmmbldKHXNE2zcoWMDpARZ2dnqVSpktExNE3T8o2DBw9eF5GyGb1mkYW+UqVKBAUFGR1D0zQt31BKnc/sNb3rRtM0zcrpQq9pmmbldKHXNE2zcrrQa5qmWTld6DVN06ycLvSapmlWLtuFXillq5Q6pJTakP64tFJqq1IqPP2nUwbrVFBKbVdKHVdKhSmlRpgyvFZwJCUlcffuXQBOnjzJjh07GDt2LM888wyzZs0iISHB4ISaZrmeZEQ/Ajj+wOMxwO8iUh34Pf3xP6UA74qIJ+AHDFVK1c5pWK3g2bVrFw0bNqRYsWI4OztTtWpVvL29GT58OMnJyZQuXZrFixfTvXt3UlNTjY6raRYpWwdMKaXcgU7AZ8A76U93A1ql318I7ABGP7ieiEQAEen345RSx4HywLFc5tYKiAULFtCqVSuCgoKIiooiJCSEli1bUrRo0fvLpKSk0LJlS6ZOnUqxYsVo1aoVdevWvf+6iHD16lVCQ0MJDw/Hz88PX19fI7qjaYZQ2bnClFJqJTARKAGMEpHOSqlbIlLqgWVuisgju28eeL0SEAjUFZHYDF5/A3gDwMPDw+f8+UwP8tIKiKSkJGrVqsXSpUvx8/PLctkjR47w7LPP0qhRI/bv34+/vz8xMTEcOnSIY8fujSvq1q1LuXLl+Omnn4iJicHR0TEvuqFpeUIpdVBEMhzBPHZEr5TqDESKyEGlVKscBnAAVgFvZ1TkAURkLjAXwNfXV1/fsIBLTk5m5MiR1KhRg8aNGz92eS8vL6KiogDYsmULa9eupXLlynTp0gUvLy/KlSuHUorAwEBCQkJ0kdcKlOzsumkGdFVKdQSKAo5KqcXANaWUm4hEKKXcgMiMVlZK2XGvyC8RkdWmCq5Zl6SkJPz9/Zk/fz4uLi4kJSXRqFEjvv/+e5RST7Stdu3a0a5duwxf8/b2Jjw8nMTERIoUKWKK6Jpm8R5b6EVkLDAWIH1EP0pE+iulvgReBial/1z3z3XVvXfo98BxEZlqutiatWnXrh0ODg5ER0dz584dChUqdH8UbkpFixalUKFC3L17Vxd6rcDIzTz6SUBbpVQ40Db9MUqpp5RSm9KXaQYMAJ5VSoWk3zrmKrFmVUSEKVOmcPXqVdatW4eTkxPly5fHxcXF5EUe7hX6Dh06sHHjRpNvW9Ms1ROdplhEdnBvdg0icgN4LoNlrgAd0+/vAkz/btXytdTUVLZu3cqWLVvYsmULdnZ2bN26FVtb2zxp39XVla1bt9K3b988aU/TjKaPjNXy3Mcff8zbb7+Ns7Mz8+bN4+DBg1SoUCHP2h81ahQ//vjj/dk4mmbtLPLCI5r1evfdd5k7dy579uzBy8vLkAw1atTglVdeoU6dOly8eBF3d3dDcmhaXtEjei3PpKWlMXXqVG7fvs2uXbtISUkxLEurVq1wcnIiNjbD2b6aZlV0odfyjI2NDRcuXOB///sf06ZNY/jw4YZluXbtGjdv3mTp0qWGZdC0vJKtI2Pzmq+vr+hrxlq3nTt30qdPH6pWrcqKFStwcXHJ0/YdHR35448/qF+/vllm92haXsvqyFg9otcM0bx5c0JDQ6lSpQr9+vVj4sSJHDhwIE/avnPnDg4ODoiILvJagaALvWaY0qVLM2vWLNq3b8+ZM2fo1q0bbm5u/PLLL2Ztd8WKFZQuXZp69eo90Xrr1q1j6NChWOJfwZqWFT3rRjNU8eLFee+994B757eZNm0a/v7+NGvWjJkzZ+LklOl58nIkOTmZhQsX0r179yeat797924GDBhAXFwcV65coX///vTq1cuk2TTNXHSh1yyGnZ0dgwYN4tatWwQGBjJixAgWLlxo0t0rhw8f5tSpU2zZsiXTZdLS0li+fDmLFi2ifv36HDlyhH379rFkyRL8/PzYsGEDw4cPp0KFCjRq1Mhk2TTNXPSuG82ilC1bls8//5xhw4bxww8/mHwKZqlSpbhy5QqtW7dm4sSJ9w+aSktLY/PmzQwaNAhvb28mTZpEu3btSE1NpVevXly8eJEuXbpQtmxZChUqREJCAiVLljRpNk0zFz2i1yzS9evXGTJkCHZ2dibdbrVq1UhKSmL58uXs27ePZs2aMWPGDHbu3MnGjRtp0qQJ8+bN4+mnn6ZQoYzfHvPmzWPmzJnUrFnTpNk0zVx0odcs0oULF3BwcDDLtu3s7BgwYAADBgygb9++NGvWjPr16xMaGkqZMmWyXHfnzp2cPn2a7t27myWbppmD3nWjWaSIiAh27NiBl5cXLi4uLFiwwCztNG3alBMnThAcHPzYIh8dHU2nTp2YO3cuxYsXN0seTTMHPaLXLFJAQADbt2+nQoUKfPDBB1y/ft1sbWV3F0yRIkVISUmhTZs2ZsuiaeagC71mkUqVKkWPHj2Ae5cJnDx5MkFBQSQnJ5OYmEhkZCTDhw+nRYsWeHh45Ekme3t7PD09OXjwIM2aNcuTNjXNFPQpEDSLFxcXx7p160hOTsbe3p7bt29z+/ZtvvnmG44fP56nV4t68803iYqKYsmSJRQrVoyUlBSWLFnC0qVLSUpK4tdff9VXrtIMoU+BoOVrJUqUoH///gwaNIgXXniByMhIFi5cyM2bN/nqq6/ytLBOmTKFIkWK4OLiwtixY6levTpTp07l3//+NzExMaxcuTLPsmhadulCr+Ub58+fp0WLFqxcuZJPP/2UCxcu8M477+RphhIlSrBs2TKOHDlCREQEH3/8MYcPH6Zbt25UrFjR0FMva1pmdKHXLF5iYiKff/45vr6+dOzYkf3799OpUyeTz7F/EhUrVmTBggW88sor95/r2rUrX375JampqYbl0rSM6C9jNYsWGRlJjx49cHJyYu/evVSrVs3oSJnq168f/v7+xMfH4+joaHQcTbtPj+g1ixQVFcWrr75KjRo1aN26NevXr7foIg9w4MAB3NzcsLe3NzqKpj1EF3rNYmzYsIGhQ4fSo0cPPD09KVmyJCdPnmTChAnY2Fj+r6qvry8VK1Zk2LBhel+9ZlEs/92jWb2bN28yZswY3nzzTWrUqEG/fv04ePAgU6dOzfMrT+VGkSJFWLt2LX/99Re1atXi8OHDRkfSNOAJCr1SylYpdUgptSH9cWml1FalVHj6zwxPHK6Umq+UilRKHTVVaM06iAirVq2ibt26XL9+nb179zJixAheeOEFKlasaHS8HClVqhS//PIL//3vf2nXrh2zZs3i6tWrRsfSCrgnGdGPAI4/8HgM8LuIVAd+T3+ckQXA8zlKp1mthIQE+vXrx0cffcScOXP47rvv8uwI17zQp08f1q9fz65du6hduzY//fST0ZG0AixbhV4p5Q50Ar574OluwML0+wuB7hmtKyKBQHTOI2rWaODAgVy/fp39+/fTtWtXo+OYRePGjVm2bBktWrTg8uXLRsfRCrDsTq/8GngfKPHAcy4iEgEgIhFKqXK5CaKUegN4A7CqkZ32sGvXrjFlyhQCAwPZu3ev2U5FbClWr17Nrl27WLJkidFRtALssSN6pVRnIFJEDpoziIjMFRFfEfEtW7asOZvSDLJ582bc3NyIjY1lz549VKlSxehIZjVjxgyGDBnCzz//rKdcaobKzoi+GdBVKdURKAo4KqUWA9eUUm7po3k3INKcQbX8Lzo6mtq1axMQEEDhwoWNjmN2hw8fpmvXrjRp0sToKFoB99gRvYiMFRF3EakEvARsE5H+wHrg5fTFXgbWmS2lZhVefPFFUlNTWb9+vdFR8sTQoUP57bffjI6habmaRz8JaKuUCgfapj9GKfWUUmrT3wsppZYBe4GaSqlLSqlXcxNYy78WLFjAuXPnaNiwodFRzC4tLY1Ro0bRqVMno6No2pOd60ZEdgA70u/fAJ7LYJkrQMcHHvfJVULNapQrV47ExESaNGnCtWvXjI5jVrt27eLMmTNs3brV6Ciapo+M1fJO165dCQwMxMkpw2PrrEqDBg1ISUkhJCTE6Ciapgu9lrdu3rxJ5cqVjY5hdiVKlODDDz/kww8/NDqKpulCr+WdxMREFi9eXGAutffaa69x6tQp/vjjD6OjaAWcLvRanlm+fDkXL15k8eLFRkfJE4ULF2b8+PF8+OGHWOK1mbWCQxd6LU/Ex8ezceNGateubfVHwz6ob9++3Lx5ky1bthgdRSvAdKHX8sTYsWOJi4tjwoQJRkfJU7a2tlSvXp0DBw4YHUUrwPSlBDWzO3ToEEuXLmXjxo24ubkZHSdP/f2XzPfff09ycjK2trakpaVRqJB+62l5R4/oNbMLDQ0lPj6+QO6+uHTpEikpKVSuXJmSJUvi4OCAm5sb169fNzqaVoDoQq+Z3YABA9i0aRPTpk3j7NmzRsfJU9WrV2fDhg2cPHmSqKgoIiMj8fDw4PXXX2fDhg2kpKRw+PBh/WWtZlbKEn/BfH19JSgoyOgYJhcZGclvv/3GtWvXuHv3LmlpaezevRsR4fLly6SkpNC7d2+6du1KgwYNjI5rcv3796dOnTqMHTvW6CiGio2NZcmSJUyYMIEyZcoQHh5O06ZNWbRoEeXLl8/WNi5evMiRI0do3rw5JUqUePwKmtVTSh0UEd8MX9OF3vzS0tIICQlh1KhRREdHU79+fUqXLs3du3dp27YtNjY2ODs7c/r0aYKDg1m8eDHPPfcckydPplKlStluw9IvoL1w4UKmTZvGzp07dXECrl+/TlBQEM2bN2fMmDHExsaycOHCh5b55ZdfKFq0KJUrV6ZSpUrcuHGDLl26EB4eToUKFahevTo//vijQT3QLIku9AZKSEhgxIgR/Pjjj3Tp0oUpU6bg6uqa5Tp/z06ZM2cOJUqUoHfv3iilKFWqFGfPnqVy5cr4+vpy8OBBfvrpJ+zt7QkODmbo0KG0bNmSrVu38vrrr1OvXr086mX27Ny5kxYtWhAVFYWzs7PRcSzKjRs3aNSoEa1ataJx48bUrl2bLVu2MH/+fJydnblw4QJPP/00ERERtG3bli+//JJbt27h5eVFmTJlaNSoEfXr12fAgAGULFnS6O5oBtCF3kxCQkLYtWsXy5Yto0yZMpQqVYrWrVvTqVMntm3bxg8//MCBAwdo3LgxCxcupHTp0tnetohw7NgxAFasWIG9vT2RkZG4u7tz7tw5du7cSbNmzWjZsiUpKSns2rWL2bNnU61aNWJiYujbty/Tpk0zV9dzZOXKlcyaNYvt27cbHcUiRUREsHz5coKDgzl79izlypUjICAAd3d3EhMTWb16NQ4ODnTu3BmlFHBvVs/Jkyf5448/+Pnnn/Hy8mL69OkG90Qzgi70JpKcnMySJUuIjY3l8OHD/Prrr/j5+TF48GCSk5OJiopixYoV7Nu3Dy8vL95++20aNGhAtWrV8iRfSkoKdnZ2AEyePJn33nsvT9rNrt27dzNixAgs8f/WGly5cgVvb29+//13UlJSiI2Nxc/Pj2LFihkdTcsDWRV6PZk3m0SEIUOGcPz4cerUqUNSUhInTpx4ZF/z66+/blBCKFSoEA0bNqRWrVoMGzbMsByZuX37NmFhYRw5cgQvLy+j41idp556ipkzZ9K8eXPKlSuHi4sLhw4dolOnTkyaNImqVasaHVEziB7Rpzt+/DhTpkzhzJkzdOjQgZ49e1KmTBni4uLYvHkzs2fPpkiRIvz6668WfZrd8PBwhg0bxvXr1wkMDLSoa5XeuHGD9957jw0bNrB582arnFlkCVJSUrC1tUUpxe3bt/H29qZ///58+umnRkfTzCirEb1lT9Mwg5SUFAIDA5k9ezatW7fGxsaGbdu28fbbbxMdHc0HH3zAoUOHaNGiBWXKlKFOnTrs2LGD8ePH8+eff1p0kYd787bXrVtHmTJlaNeuHRMnTiQtLc3oWACUKVOG+fPnM2jQIHx8fNi5c6fRkaxSoUKF7u/DnzBhAoUKFeLdd981OJVmpAK16+bAgQO0b9+eKlWq4OHhwTvvvEPdunXp1asXAwYMYOLEidjb29OuXTsAUlNTSU5OpmjRogYnfzJFixZl5syZrF27liVLlhATE8OkSZOMjnXf+PHjOXr0KOvWraN58+ZGx7FqZ8+e5T//+Y+eiVPQiYjF3Xx8fMTUIiMjpVOnTvL++++bfNuWbP/+/dKgQQOjYzxi8uTJUq9ePUlKSjI6ilVzcXGRM2fOGB3jIcuXL5eXXnpJfvvtN1myZInUq1dP9uzZk+PtpaWlyd69e+XmzZumC5kPAUGSSU0tELtu5syZg7e3N46Ojrz22mtGx8lTlStX5ty5c3z11Vf88MMP7N69m9WrV/PLL78QGRlpWK727dsjIowfP96wDNYuIiKC+Ph4izqR3Pz58xk3bhyurq689dZbTJw4kZdeeolevXoRFRWV7e188sknlChRgqZNm9K/f3969eqFl5eXvkZvZjL7BDDyZsoR/dWrVwWQ9evXm2yb+c2iRYvktddek+bNm0vNmjWlbNmyUq1aNalRo4akpaUZluvkyZPi5uYmgYGBhmWwZp988okMGzbM6Bj3paWlSdWqVTMcvY8ePVpat24tMTExj93O3r17xdnZWc6ePSvPP/+8ODg4yI0bN+TLL78UJycnGT16tFy4cMEcXXjIwYMHpWnTptK5c2fp1auX+Pv75+ovk9wiixG91e2jv337Nps3byY1NZXVq1ezfv16BgwYQOfOnY2OZpgBAwYwYMCAh55LS0ujbNmyLF68+JHXAH766Se2bdtGp06dzPZvV6NGDSpVqsStW7fMsn2NXJ9qIiIigoCAAJYsWXL/EpClSpXCx8cHLy8vWrZsSe3ate9/+ZuVn3/+mZiYGPz8/B55bcKECQwbNow6deowc+ZMunXr9tDrBw8eJCwsjPPnzzNz5kwWLFhA5cqVWbp0KSJC6dKlGTRoEABnzpyhT58+7Nq1K1d9z0pSUhJDhw6lTp06dO3alYSEBE6ePEnXrl1ZsWIFcXFxFClShPbt25stwxPJ7BPAyFtuRvTt27eXRo0aScOGDWXkyJFy+fJlQ0etliwgIEAqVKggp0+fvv9cSEiItG/fXhwcHGT06NFiY2Nj1lGKs7OzhISEmG37BdnUqVNlyJAhOVr3yJEj0r9/f3FycpI333xTjh07JseOHZN9+/bJ77//LgEBATJo0CDx8PCQTp06Zfkeu3TpkgQGBsozzzwjPXv2zLLdb7/9Vpo0afLQc5GRkVKqVCnp37+/jBo1Svbt25flNpKTk8XV1VWOHz+e/Q5n0507d+TQoUPi7+8v9erVk+Tk5Ide37hxo5QtW1YAAeTo0aMmz5AZshjRZ7v4ArbAIWBD+uPSwFYgPP2nUybrPQ+cBE4DY7LTVm4KfdOmTWXRokU5Xr+g+frrr6Vu3boSGBgoa9euFXd3d/H29r7/p+/06dPF0dFRvv76a5O3HRcXJ7Vr15atW7eafNuaSHh4uJQpU0b2798vU6dOlREjRsiGDRseKsopKSny008/ybZt22TUqFFSpEgRadGihZQpU0a++uoriYiIyLKNlJQUqVmzprz//vty7NgxWbNmjWzdulUiIiIkLi5Opk6dKoULF5YKFSrIW2+9JWfPns1ye5MmTZJnn332oS/pAwICpH///k/U9yFDhshXX331ROv8LSIiQg4fPnz/cVxcnOzevVvee+89cXJyktq1a0u/fv0kODg4w/VTU1Pl6NGjUqhQocd+KJmSqQr9O8DSBwr95L8LNzAG+CKDdWyBM0AVoDBwGKj9uLZyU+jXrl0rrVq1yvH6BU1aWpp8+umn4unpKUWLFpW33npLDh069NAyo0aNkldffdWk7SYkJMjgwYOlXr16kpqaatJta/9vxYoV4ujoKM7OzvLll1+Kt7e3dOzYUW7fvi2rV68WT09P8fPzE19fX3n++efl9OnTsmDBAomMjMx2G4GBgTJgwAABpGnTplK+fHlxdHSU4sWLS506deTcuXPZ3lZUVJQ0btxYFi1aJJcvX5bXXntNXF1d5Ycffniifs+fP1/atGmTrd+tuLg42b59u2zYsEH69+8v9vb24urqKlWrVpUyZcpI0aJFxcfHR3r27JntGUzz5s174g+n3Mp1oQfcgd+BZx8o9CcBt/T7bsDJDNZrAmx+4PFYYOzj2stpoQ8KChJPT0/p0qVLjtbXMjZv3jwpXLiw3L171yTbi4qKkrZt20qrVq3k2LFjJtmmlrmwsLD7u8dSUlKkcOHCUqxYMWndurWsWrXKZLs2L168KKmpqRITEyOxsbESHR2do9+Zjz76SBo1aiQeHh4yZswYCQ4OfuKMN2/elGbNmmVabGNjYyUsLExWrlwp7u7uUq9ePWnfvr2MHz9e4uLiJD4+Xk6cOCFXr17N0b/P559/LgMHDnzi9XLDFIV+JeADtHqg0N/6xzI3M1jvBeC7Bx4PAGY+rr2cFvr+/ftL1apVZd26dTlaX8vY7t27pXjx4iYZeYeHh4urq6u89dZbj+zf1PLGiRMnHrtLxkiRkZEyYsQI+e2333K1nfj4eHF2dpY6deqIv7+/pKSkiIjId999Jw4ODlKzZk3p2LGjLF++3BSxH7Jz505xdHSUuLg4k287M7kq9EBnYHb6/Sct9C9mUOgDMmnnDSAICPLw8MhRR8PDw+WFF16Qp556KkfraxmLiooSOzs7kxzctGPHDqldu7beXaPliRs3bsi2bdukQYMG0rJlS/nqq6/Ezc1NTp06ZdZ2U1JSpHfv3jJ8+PA8mwyS20I/EbgE/AVcBe4Aiy1x142IyF9//SUODg4SEBCgZ9uYSHJystja2kp8fHyutxUbGyuA3mWj5amUlBTx9PQUwCyzcTJy6dIlcXV1lXHjxuVJe1kV+sceGSsiY0XEXUQqAS8B20SkP7AeeDl9sZeBdRmsfgCorpSqrJQqnL7++se1mRsVK1Zk/fr1zJkzh6VLl5qzqQKjUKFCdO7cmcmTJ+d6W8HBwdSrV49atWqZIJmmZY+trS3Lly9n4cKFefa7V758eSZPnsyff/6ZJ+1lJTenQJgEtFVKhQNt0x+jlHpKKbUJQERSgGHAZuA4sEJEwnIX+fFat25NQEAAn3zyCUlJSeZurkD44IMPTPLBeejQIWrVqpWtA2w0zZS8vb0ZOHBgnrW3fPlyBg4cSJ8+ffKszcw8UaEXkR0i0jn9/g0ReU5Eqqf/jE5//oqIdHxgnU0iUkNEqorIZ6aNn7nWrVtTuXJl5s2bl1dNWjUfHx/Cw8NZs2ZNrrbTtWtXdu/eTUBAgImSaZplOn36NMOGDeOVV14xOor1ntRMKcWXX37Jxx9/zPr1Zt1bVCDY2tqyZcsWBg8eTP/+/fntt99ytJ0qVaoQGBjIzJkzmThxoolTaprlKFKkCAkJCUbHAKy40AN4enpSv359QkNDjY5iFdq2bcv3339PbGwsEyZMyPF2KleuzLZt25g9ezZ79uwxYUJNsxwrV67kxRdfNDoGYKWF/n//+x/u7u5UqlSJ4OBgypcvb3Qkq9GzZ08mTJiQ61Mcly9fnmHDhrF8+XITJdM0y3L06FFq165tdAzACq8wdfHiRQYPHsyyZcsoXLgwPXr00F/8mUFSUhIikqt/2woVKljEjARNM4cXX3yRtWvX8tZbbxkdxfpG9H9fDNvLy4uePXvqIm8GXl5eFCtWjE2bNuV6O2FhYX8fY6FpVqVs2bLEx8cbHQOwwkK/detWqlSpgqenp9FRrJZSismTJ/PSSy9RqVIlxowZk6NiXadOHQACAwNNHVHTDFenTh32799vdAzACgv9+vXr6d+/PzY2Vtc1i9KhQwdiYmLYvHkz27dv5+OPP37ibdjY2PCf//yHd955h9TUVDOk1DTjeHl5sXHjRvr27UtYmNkPH8qSVVXDtLQ0li5dipOTE8nJyUbHsXo2NjbUrFmTDRs2MHXq1Bx9QTtgwADs7OxYty6jA6s1Lf/y8fHhyJEjVKpUifbt27Nq1apMl71x44ZZs1hVobexseHjjz9mwYIFDB8+3Og4BUZaWho2NjaULl36iddVSvHSSy/xyy+/mCGZphmrRo0afP755yxcuJCxY8fy7rvvcurUKWbPnk2DBg1o1aoVtWvXxtnZOdffeWVFWeIXYb6+vhIUFJTj9WNjY6lSpQr169dnxYoVOSpAWvZ98MEHREVFPfFRyCkpKaxbt45Ro0YxfPhwRo4caaaEmma8a9eu8f7777Nz506cnZ356KOPcHBwwMHBgd9//53IyEimTp2a4+0rpQ6KiG9Gr1nd9EoAR0dHzp8/z9NPP01wcDBt2rQxOpJVu3nzJgcOHCAkJIT69etna50DBw7QvXt3PDw8mDVrFh07dnz8SpqWj7m4uLBw4cIMX9u2bRu3bt0yW9tWtevmQfb29pQqVQo7Ozujo1i9WbNmUatWLVasWJHtdd5//30mTJjA3r17dZHXCryQkBDc3d3Ntn2rHNHDvQOn9u7dS4UKFYyOYvX+nuFUs2bNbC0fGRlJSEgIv/76qzljaVq+kZycTOHChc22fasd0bu6utK5c2eGDh1KXFyc0XGs2pEjR9iyZQtPP/10tpY/ceIEVapUoUiRImZOpmn5Q6VKlYiJiSExMdEs27faQm9nZ8eiRYu4dOkS8+fPNzqOVZs3bx7dunXL9kFqTz/9NPHx8fpUxZqWzs/PjylTpvDee++ZZftWu+sGwMnJicTERMqVK2d0FKt14MABfvrpJxYuXJjt000UK1aMX3/9ldatW2NnZ8eQIUPMnFLTLNe5c+f45ptv8PT05L///a9Z2rDaEf3fzp07x44dO4yOYbVmzpxJgwYNaNeu3ROt9/epE3bu3GmmZJqWP0yaNImyZcuybds2SpYsaZY2rL7QX758mZUrV3LkyBGjo1gVEWHu3Lls2rQpx5dKO3bsWLanY2qatTp48CA+Pj64urqarQ2rL/TlypWjadOmFnNyIWvxxx9/8Nlnn7F9+3YGDBjwxOtHR0fz008/0apVK9OH07R84syZM5w/f97s15W1+kIP0Lt3b4YPH877779vdBSrsWTJEvz9/albt26O1v/tt9/w9fXN9kwdTbNGBw8epH79+ma/OJJVfxn7t/79+wP3TqA1adIkfWbLXBIR1qxZQ25OU3H69Glq1aplwlSalv8sXbqUzp07m72dAlPxfH19qVOnDo0aNWLOnDmcP3/e6Ej51s6dOylVqhQVK1bM8TbKly/PuXPnTJhK0/Kfjh078sEHH1ClShVefvlls+1iLjCFvlatWoSGhjJ06FD+/e9/8+677xodKV/asWMHvXv3ZsqUKbm6etfdu3dxdHQ0YTJNy3/eeOMNrly5wqeffsqiRYvMNkPwsYVeKVVUKbVfKXVYKRWmlBqf/nw9pdRepdQRpdTPSqkM37VKqRFKqaPp675t4vxPxMbGhkGDBuHp6UmZMmWMjJJvDR06lICAALp3756r7YSGhuLl5WWaUJqWj5UsWZK9e/cyevRos32PmJ0RfSLwrIjUA+oDzyul/IDvgDEi4gWsAR45pEspVRd4HWgE1AM6K6Wqmyh7jo0bN47g4GCjY+RL8fHxNGjQINfbOXPmDFWrVjVBIk3Lf6Kiojh+/Djjxo3jiy++YOPGjcae60buuZ3+0C79JkBN4O+LfW4FemWwuifwp4jcEZEU4A+gR65T51JgYKC+IHUOzJ8/n+TkZJMcaWxvb8/p06dNkErT8pf4+HgqVapE586diYyMZN++fbzyyit89NFHZmszW7NulFK2wEGgGjBLRPYppY4CXYF1wItARqeJPAp8ppQqAyQAHYGcT9UwkbNnz+Lv7290jHwlNTWVcePGsWrVKkqUKJGrbYkIN2/exNbW1kTpNC1/SE5O5oUXXqB37955eg6ubH0ZKyKpIlIfcAcape+SGQwMVUodBEoASRmsdxz4gnsj/l+Bw0BKRm0opd5QSgUppYKioqJy0pdsS0pKYvbs2XpU/wR27dpFYmIifn5+ud7WtWvXCAkJ0R+2WoFz6NAhQkND+fbbb/O03SeadSMit4AdwPMickJE2omID7AMOJPJOt+LSEMRaQFEA+GZLDdXRHxFxLds2bJPEuuJbd68mYSEBH0+9CdQs2ZNbty4YZKr4Dg5OelTR2sFzscff8yyZcu4cuUKc+bMydO2H7vrRilVFkgWkVtKqWJAG+ALpVQ5EYlUStkA/wG+yWT9v5fzAHoCTUyYP0fs7OyoWbMma9asoVy5cvj4+BgdyeKVLl0aEaFQoXu/MgkJCQQEBNCmTRsaNmz4yPLR0dHcunWLKlWqcOrUKebMmcPdu3exs7MjISEBe3t7kpOT9RXAtAIhNTWVCRMm0LdvX5o0aZLnp/7IzojeDdiulAoFDgBbRWQD0EcpdQo4AVwB/geglHpKKfXg5cxXKaWOAT8DQ0Xkpkl7kEMVK1Zkw4YNPPvss3z55ZdGx7F4sbGxlChRgtjYWA4fPoy3tzfffvst77zzziPLigidOnWievXqzJgxgyVLlrB+/Xpq1apFpUqVqFmzJoGBgRQvXtyAnmha3ktMTMTDwwN/f3/27NmDt7d33gYQEYu7+fj4SF45c+aM1K5dW3r37i1paWl51m5+ExYWJo6OjhIcHCwuLi6yePFiSUxMFFdXVwkLC3to2Y0bN0rt2rXl9OnT0qhRIwFk7NixBiXXNONNnz5dmjRpIgkJCWZrAwiSTGqq4UU9o1teFnoRkcuXL4udnZ3cvXs3T9vNb5ydnQWQxYsX33/uo48+kueff/7+v11qaqrUrVtX1q5dKyIiaWlpEh8frz9EtQKtW7duMmXKFLO2kVWhLzCnQMhKYGAgycnJHD582OgoFq169ep06dKFfv363X/uww8/xN7envbt27Njxw4mTJiAg4MDXbt2BUApRfHixXN1ugRNy8+2b9/O7t27efPNNw3LUCDOXvk4vXv3JigoiClTprBixQqj41isnTt3PnLmzyJFirB8+XLmzJnDqFGjcHV1ZenSpbqwa1q6mJgYKleujIODg2EZ9Iiee6PONWvWUKNGDaOjWDRbW9sMC3ihQoV46623CAoKYsOGDVSuXNmAdJr25BITE83exp49eyhVqhTR0dFmbyszutCn++ijj9i2bZvRMTRNM5OoqCiGDx/O888/T5MmTWjRogVFixbl+++/N2u79vb2/Pnnn3Ts2NGs7WRFF/p0qampODk5GR1D0zQT+9///sfrr79O7dq1SUpKYtiwYXz11Vd88sknrFq1ioCAgPvLJiYmsnHjRhYsWGCyI+c//vhjzp07x40bN/juu+9Mss0nltm3tEbe8nrWjYjIqlWrBJCvv/46z9vWNM18HB0dZcSIEfL7778/8tqdO3ekdOnS8tlnn8mSJUukUqVK0rx5c3F1dZXPP//8oWVXr14tQ4YMkbNnzz6ynePHj0tMTEyWOYKCgqRKlSryr3/9SxITE3PXqQygp1dmz4IFC6Rly5aGtK1pmnn4+fnJmjVrMn394MGDMmLECOnRo4fMmjVLRERCQ0NFKSXDhw+XY8eOiYhIhQoVxMPDQzp06CAi9wr3yy+/LG3atBEHBwdp3Lix3Lp1K8ssCQkJUr9+fdmwYYNpOveArAq93nXzgNq1a3P69GlCQkKMjqJpmokUKlQoyzOuNmzYkK+//prVq1fz73//GwAvLy/Onz+Po6MjLVu2xN3dnZiYGD777DOCg4M5evQoY8aMwdbWlgoVKhAeHk7lypXp1asXaWlpGbYTGxvL+++/z4ULF/L+WgyZfQIYeTNqRC8iMm7cOPH29jasfU2zVsePH5fRo0dnuAvFnHx9fWXlypU5Xj8xMVGOHTsmFy5ckCNHjkizZs3Ezs5OevbsKbGxsfeXS0lJkcqVK8uBAwckKSlJkpKSHtrO22+/LSVLlpTQ0NAcZ8kKWYzo9Tz6f3BwcDDJFZQ0Lb9LS0tDKUVqaio2NjaPHEORkfj4eOzt7TN8bejQoXh4ePDqq69Sv359unbtyp07d3B3d6dq1aoUKVKEKlWqmPw6BePHj2fUqFH06NEjW334p8KFC+Pp6QlAhQoV2LVrF2lpaY9sy9bWljFjxtC+fXsAypUrx+LFi/n666+JiIjg6NGjhIWFUb58+dx36gnpXTf/UKVKFS5fvkxqaiqXL182Oo6mGSIiIgJbW1vs7OwoVqwYnp6eHD58mNu3b/PXX3+xfft2GjRowPjx49myZQvdu3enevXqODg4kJyc/NC2RIQrV64QEhLC119/TVhYGK1atWL+/PkcOXKEefPm8a9//Qtvb29cXV0ZM2YML774IgsWLDBJXzp06EBUVBSbNm16/MLZlNkHxhtvvEFoaCiHDx8mNTWVRo0a0bBhQ4YMGcK+ffsMKfKA3nXzT4cOHZLChQtL+fLlBZCPPvpIn6dFK3ASExOlWLFicuPGDUlMTJQpU6ZI1apVpVixYsK9S4mKi4uLNGvWTABp1aqV/PjjjwLIe++9JxcvXpS4uDh58803pXjx4lK6dGl57bXXsmwzISFBBg8eLIMGDZIZM2YIILt27Xri7GlpaTJ8+HDp3r27DBkyRNq1a3c/84ULF3L6T/LETp8+LUFBQXnWHlnsulH3Xrcsvr6+EhRk3BUHz507R1RUFCVKlKBVq1bs3r2batWqGZZH04zg5+dH7969GTly5P3nbt++TVJSEnfu3MHNze2R3Szffvst06dP5/jx4yilaN68OXPnzqV48eK4u7tn+9QYMTExlCpVivPnz+Ph4fFEuU+cOIGnpydjx46lXLly1KxZk2vXrrFkyRJmzJhxfzeMtVFKHRQR3wxfzOwTwMibkSP6f+rQoYN89913RsfQNJNJS0uTrVu3yqJFix76MvFBd+7ckYCAAKlTp06O2pg7d658//33uYkplSpVuj+1MTvu3r0rU6dOlTJlysi7774rqampuWo/v0HPo8+ZPXv2SNmyZWXfvn1GR9E0k5k6daq4u7tLyZIlpX///o+8HhkZKb6+vgKIp6enAQnv6dixo3h4eMg333zz2GX3798vnp6e0qVLF/njjz/yIJ3l0YU+B9atWyfA/fOqa1p+dezYMenQoYO88cYbMnLkSHF2dpajR49KeHi42NraPjL1cP/+/QJIdHS0xMXFGZT63rUN9u7dK4C0b99etm3b9sgyaWlpsn79evHx8ZGPP/64QH+fllWh17NuMlGkSBEAmjVrZnASTcudZcuWERsbS/369SlRogSrVq2iTp06VKtWje+++45ly5bdXzYyMpLFixfTpEkTnJycjD21ro0Nfn5+7Nixg759+9K3b18WLlz40DJhYWEMHDiQZs2aMXToUH167EzoefSZaN++PS1btuTgwYP358VqWn7UpEkT/vjjD/z9/R95rXr16pw7dw6AF154gVWrVvHqq6+yZMmSvI6ZqZYtWwLg6+tLmzZteOqpp2jbti0Af/31Fz4+PkyfPt3IiBZPF/osFC5cONPDmTXN0l26dIkTJ06wadMmwsLCMlzG3d2d8PBw9u3bx61bt5g9e3aGHwiWoHbt2ixatIiBAweyaNEi6tWrx7Jly+jcubPR0SyeLvRZcHJyIjw8nA4dOhgdRdMytH37dqZMmUKHDh3466+/uHLlCgkJCVy7do2TJ09Sr149UlNTGTVqVIbrV6xYkc8++4w+ffpQuHBhevXqlcc9eDJt2rRh7ty590f0AK+//rqBifIHPY8+C+vWrWPq1Kn88ccfDz1//vx5ihQpgqurq0HJNEt3584dfv75Z1JSUqhSpQoNGza8/72PqezZs4cuXbrw7rvvcuTIEVxdXfHx8aFYsWK4uLjg6elJmTJlsrWtv+tAftnHnZaWxsWLF7l48SKNGzfGzs7O6EiGy2oevS70WYiPj8fLy4vWrVszfPhw6tWrR2xsLFWrViUlJYVq1arh4+PDiBEjrPYgDC1n+vXrx+XLl3Fzc+PkyZMcO3YMR0dHqlWrxq5du3J0zpW/iQg///wzQ4YM4bvvvjP0ykWa5ciq0D/2t00pVVQptV8pdVgpFaaUGp/+fD2l1F6l1BGl1M9KKcdM1h+Zvt5RpdQypVTR3HUn79jb27Nz506cnZ3p0KED/v7+dOnSha5du7Jv3z4CAgJwc3OjefPm/Pjjj0bH1SzIzp07+f7771m2bBnBwcHcunWL8ePHc+zYMXIzuLpx4wZDhgxh7NixzJkzRxd5LXsym3f59w1QgEP6fTtgH+AHHABapj8/GPhvBuuWB84BxdIfrwBeeVybljCP/p/OnDkjgwYNksmTJz8ytzgkJERcXV1lx44dDz2/bds28ff3L7AHcBRkbdq0kdWrVz/0XGhoqDg4OMjWrVufeHsXL16UBQsWSLly5WTAgAFy/fp1U0XVrASmOmAKKA4EA42BWP5/108F4FgGy5cHLgKluffF7wag3ePascRC/zhr1qwRDw8PuXr1qoiIhIWFibOzs3z66afi4uIiu3fvNjihlldOnDghFStWzPAUAGPHjhVfX99sbefUqVPSvXt3qVu3rjg7O0uHDh1kz549po6rWYmsCn22Zt0opWyBg0A1YJaI7FNKHQW6AuuAF9OL/T//WrislJoCXAASgC0isiU7beY33bt3JzQ0lCpVquDs7MydO3f49NNP8ff3p0SJEnzyySd88cUX+lz3BcC0adPo0aMHgwYNeuS1woULc+7cOY4fP06tWrUy/PJTRPjhhx94++236dixI3PnzqVx48a52q+vFWzZ+s0RkVQRqQ+4A42UUnW5t7tmqFLqIFACSPrnekopJ6AbUBl4CrBXSvXPqA2l1BtKqSClVFBUVFSOOmO0cePGERERwfbt2/n9998ZMmQIAC+//DK2trY888wzD11xXrNOHh4e3LlzJ8Mi3qhRI1xdXXn++edp1aoVZ8+eBeDMmTOEh4fz5ptvUrFiRSZNmsTmzZvvH6Wqi7yWK5kN9TO7AR8Do/7xXA1gfwbLvgh8/8DjgcDsx7WRH3fdZMe5c+ekVq1a8s477xS4M+sVJI0bN5bNmzdnuUxqaqpMnDhRypQpIxUrVpSyZcuKg4ODDBs2TE6cOFGgz9mi5Qy5OdeNUqqsUqpU+v1iQBvghFKqXPpzNsB/gG8yWP0C4KeUKq7uDW+eA47n/GMpf6tUqRJ79uwhKCiI3r17k5KSYnQk7QmICG+88QYzZszIcuaMUirTy+n9zcbGhjFjxnDx4kUCAgIIDw8nLi6OgIAAatasmW/ms2v5Q3b+HnQDtiulQrk302ariGwA+iilTgEngCvA/wCUUk8ppTYBiMg+YCX3vsA9kt7eXJP3Ih9xcnJiy5YtREdH88UXXxgdR3sCiYmJzJs3j88++4yWLVuye/duUlNTH1muRIkSxMTEZGubxYoVo0uXLpQsWdLUcTXtPn3AlEEuXrxIw4YN2bp1K/Xr1zc6jpYNO3bs4LXXXiM0NJRFixYxffp0bt68SceOHfnss89wc3MjISEBBwcHLl68yFNPPWV0ZK0AydUBU5p5VKhQgddff/2hU8RqlunUqVNcvHiRlStX4u3tTfHixRkyZAjHjx9nz549JCUlMXXqVADs7OwQEdzc3AxOrWn/T5/UzEBNmjRhxowZRsfQsnD16lW8vLwoXbo0xYoVY9euXQ+9XqVKFTw9PYmNjQXuHblaqlQpvY9dsyh6RG8gNzc3rl+/bnQMLRNxcXF89tlndOzYkYiICM6ePZvh7ph69erx448/EhYWRlRUFC4uLgak1bTM6RG9gSpWrMj58+eNjqFlYsuWLcycOZM5c+ZkuVznzp25fv06zzzzDMWKFaNEiRJ5lFDTskeP6A1UvHhxbt68yfLly7ly5YrRcbR/6NWrF0uXLmXEiBGPXfaVV14hIiICLy+v+/vrNc1S6EJvIHt7e77++muWLl2Kn59frs5qqJnHwYMHqVGjRraWLVq0KJs3b6ZTp05mTqVpT0YXeoONGDGCdevWYWNjk+nl3jTjdO/enejoaI4fL7DH+WlWQBd6C6CUomfPno9c4V4zXt26dSlSpAjz5883Ooqm5Zgu9BbirbfeYv78+XpUb2Hmzp1LtWrVmDBhgtFRNC3HdKG3ACLChAkTiI6OZt++fUbH0R5QvXp1rl69avLrvWpaXtLTKw0mIqxYsYJffvmFc+fOUalSJaMjaQ8oXbo0V69eZc+ePTRt2tToOJqWI7rQG2T48OH8/PPPFCpUiLS0NKZNm6aLvAVq2bIlX3zxBQMGDOD06dP6iFctX9InNTPAnTt3sLe3Z+bMmbRq1YrKlStTvHhxo2Npmbh9+zalSpXi/PnzlC9f3ug4mpYhfVIzC2NjY0O9evU4efIkderU0UXeQl26dIlBgwbh4+PDwIED9dkotXxLF3oDFC1alO3bt7NhwwY9pdJCpaWl0alTJ1xdXRkzZgzffvut3m2j5Vt6H71BnJyc+Pnnn3n++edJSkri9ddfNzqS9oBvvvkGR0dHPv/8c13gtXxPj+gNVKdOHbZs2cLo0aOJjo42Ok6BsmfPHiZPnpzhFaIA1q5dy+jRo3WR16yCLvQG8/T0pGPHjixYsMDoKAXKv//9byZOnMjIkSOJiYnh7t27908ZHRQUxKFDh/D29jY4paaZhi70FsDf358PPviA5ORko6NYNRHhypUrrF27lsOHDxMcHMytW7fw9PTkmWeewcXFhUGDBtGhQwe+/fZbPDw8jI6saSah99FbgJIlS1K4cGFsbPTnrjm9/fbbzJo1i5IlS9K1a1cqV67MokWL2LNnD6GhoTRp0oRdu3bRvXt3unXrZnRcTTMZXegtwIULF6hfvz62trZZLhcXF0dgYCCtW7fWUzKzKSwsjGXLllG6dGmWLVvGtWvXKFOmzEPLNG3a9P5Rr/Xq1TMipqaZlS70FuDSpUtERUWRmJhI4cKFiYqKYtWqVbi4uHDs2DGuXLmCra0tixYtwtXVlWefffaxVz3S4Mcff+Tdd9/F29sbR0dHPvnkk0eKvKYVBLrQW4DXXnuNNWvWULZsWRISErCxsaFu3boULVqURo0a4eHhwZ07dwgNDaVIkSLUrFkTDw8PGjZsSLt27fTMkAykpKTQt29fJk2axDvvvPPYv5Y0zZrpQm8BbGxs2LRpE5GRkZQsWZKiRYtmufzq1auZPn06kydPZtGiRXTp0iWPkhrjwIEDXLp0iR49ejx2WRFh2bJlbNmyhbS0NP766y9d5LUC77Hf/imliiql9iulDiulwpRS49Ofr6eU2quUOqKU+lkp5ZjBujWVUiEP3GKVUm+boR/5nlIKFxeXxxZ5gOeee47Vq1cTHx9v1Yfliwi7du2iY8eO9OzZk1OnTj12ndmzZzNp0iT8/PzYvn07U6ZMyYOkmmbZsjOiTwSeFZHbSik7YJdS6hcgABglIn8opQYD7wEfPbiiiJwE6gMopWyBy8AaE+YvsAoVKkTPnj1ZsGABDRo0sKoZOyLCpEmTWLZsGVFRUcyZM4dbt27RvHlzRo4cSUREBDVq1KBHjx489dRTxMfHU7x4ce7evcv48eP57bff9Bx4TXvAYwu93Du95e30h3bpNwFqAoHpz28FNvOPQv8PzwFnROR8jtNqD5k+fTovvvgiTZo04ccff7Sa0xxfu3aNcePGsXDhQl566aX7H2IpKSns3buXunXrsmfPHv7zn/+QmJhIQkICtra22Nvb89xzz+kir2n/kK3TFKePxg8C1YBZIjJaKbUH+EJE1iml3gHGi0iJLLYxHwgWkZmZvP4G8AaAh4eHz/nz+vMgO1JTUxk5ciRnzpxh/vz5uLi4GB0p10SEF154gerVqzNp0qRMl0tJSeHMmTPY29vj4uJCdHQ0zs7Oep+8ViDl+jTFIpIqIvUBd6CRUqouMBgYqpQ6CJQAkrIIUBjoCvyURRtzRcRXRHzLli2bnVgaYGtry6RJk3B1dcXV1ZWTJ08aHSnXlFJ06NCBtWvXZrlcoUKFqFmzJu7u7tjZ2eHi4qKLvKZl4Il27IrILWAH8LyInBCRdiLiAywDzmSxagfujeav5TSolrnixYszc+ZM3NzcuHPnjtFxcu3ChQtMmTKFTz75xOgommYVsjPrpqxSqlT6/WJAG+CEUqpc+nM2wH+Ab7LYTB/ufRhoZrJx40Zq1aqFk5MTx48fNzpOrgQEBNC6dWt69+5tdBRNswrZGdG7AduVUqHAAWCriGwA+iilTgEngCvA/wCUUk8ppTb9vbJSqjjQFlht6vDa//vzzz+JiIigcePGNG/eHF9fX4YOHUqdOnXw9/fn3LlzpKSkGB3zsc6ePcuaNWt48cUX9YFgmmYi2Zl1Ewo0yOD56cD0DJ6/AnR84PEdQB93bmbPPfccN2/epFmzZvTu3Zs///yTP//8k9mzZ7Np0yZ8fX25c+cOly5dsujTALz33nt07NiR1q1bGx1F06yGvjh4AVKzZk1++ukni55++NJLL9GtWzf69OljdBRNy1f0xcE1AJo3b84XX3xhdIwsOTg4cPv27ccvqGlatulz3RQgM2bM4KmnnuLq1au4uroaHecRqamp3LhxI9PL+2maljO60BcgxYsXp3nz5qxZswZ/f3+j43D79m3Wr1/PkSNHaNKkCYcOHeLMmTP6oh+aZmK60Bcw77//Pt27dyc2NpaqVavi7e1NjRo18jRDdHQ0gwcPZvv27dStW5d27doxevRoIiMjCQ4Oxs3NLU/zaJq101/GFkCHDx/m888/59atWwQHBzN06FD8/PwoVKgQzz33nFmmNcbExODt7Y2XlxdhYWF07doVPz8/unbtir29PSkpKSQlJekrZ2laDmX1Zawu9AXc0aNHmTdvHseOHePUqVN8/vnn9OvXz2TbFxHOnj3L4sWL2bx5M/369aNEiRIMHDjQZG1omqYLvZZNGzduZNy4cezatYtixYoBkJaWxtKlSylatCiNGzemQoUKT7TNsWPH8v3331OuXDmmTZtG27ZtzRFd0wq8rAq93kev3de2bVt++OEHnJycqF27Ni1atGDLli2UKFECe3t7Bg8ezMCBAxk+fDhVq1YlKSmJO3fuULp06Ux392zZsoV169bRpEmTPO6Npml/0/PotfsKFy7M8uXLiYmJ4b333sPJyYkZM2bw559/sm3bNkJCQoiPj6dly5b4+PjQoEEDypYtS8WKFYmPj89wm05OTuzcuTOPe6Jp2oP0rhvtiaWlpfHLL78QERHBq6++iru7O+vXr8fHx+eRZbdv386zzz5LWlqaPneNppmRPjJWMykbGxs6derEa6+9hlKKN998kz59+hASEnJ/mePHjzNu3DieffZZhg0bpou8phlIF3ot1z744AOGDx9O+/btmT17NvHx8TRp0oSEhARWrlypL9CtaQbTu240kwkODsbf35+goCBatmzJtm3bjI6kaQWGnnWj5YmGDRsSGBhIdHS0PrpV0yyI3nWjmVSRIkV0kdc0C6MLvaZpmpXThV7TNM3K6UKvaZpm5XSh1zRNs3K60Guaplk5Xeg1TdOsnC70mqZpVs4ij4xVSkUB53O5GWfgugniGEn3wXJYQz90HyyDufpQUUTKZvSCRRZ6U1BKBWV2OHB+oftgOayhH7oPlsGIPuhdN5qmaVZOF3pN0zQrZ82Ffq7RAUxA98FyWEM/dB8sQ573wWr30Wuapmn3WPOIXtM0TSOfF3ql1I9KqZD0219KqZD05ysppRIeeO2bLLbxllLqpFIqTCk1Oc/C/3/7ue5D+vKjlFKilHLOk+APt52rPiilvlRKnVBKhSql1iilSuVl/vQMue1DaaXUVqVUePpPpzztAJn34YHXPZRSt5VSozJZv75S6s/09YOUUo3yJPijOXLVj/RlLPJ9/cDrj+1D+nImeV/n6wuPiEjvv+8rpb4CYh54+YyI1M9qfaVUa6Ab4C0iiUqpcmYJmoXc9iF9vQpAW+CCyQNmgwn6sBUYKyIpSqkvgLHAaJMHzYIJ+jAG+F1EJimlxqQ/tqQ+AEwDfsliE5OB8SLyi1KqY/rjVqbO+Ti57Uc+eF/D4/8vTPq+ztcj+r+pe1ee/hew7AlX9QcmiUgigIhEmjpbduWiD3Dvl+Z9wNAvXHLaBxHZIiIp6Q//BNxNnS27cvH/0A1YmH5/IdDdhLGeSEZ9UEp1B84CYVmsKoBj+v2SwBUzRcyWXPTDot/X2ewDmPB9bRWFHmgOXBOR8Aeeq6yUOqSU+kMp1TyT9WoAzZVS+9KXe9r8UTOVoz4opboCl0XkcJ6kzFpO/x8eNJjHjHTMLKd9cBGRCID0n3k+inzAQ31QStlz76+L8Y9Z723gS6XURWAK9/6yMlJO+2Gx7+vs9sHU72uL33WjlPoNcM3gpQ9FZF36/T48PAKLADxE5IZSygdYq5SqIyKx/9hGIcAJ8AOeBlYopaqIiacimasPSqniwIdAO1PmzYiZ/x/+buNDIAVYYsLoD27f7H0wtxz2YTwwTURu3xtgZsofGCkiq5RS/wK+B9qYIPYjzNwPS35fP7YPZnlfi0i+vnHvP/Ua4J7FMjsA3wye/xVo9cDjM0DZ/NIHwAuIBP5Kv6Vwb3+ea37pwwOvvQzsBYrn09+lk4Bb+n034KSl9AHY+cDvyC0gGhiWwbox/P+UawXEWtL/xRP0w2Lf19npgzne1xY/os+GNsAJEbn09xNKqbJAtIikKqWqANW5t0/sn9YCzwI7lFI1gMIYc8KkHPVBRI7wwC4CpdRf3CtC+aYP6cs9z70/Z1uKyJ28CpyB3Pwurefeh9Wk9J/rMlgmLzzSBxG5v7tJKfUJcFtEZmaw7hWgJfc+zJ4FwjNYJq/kph9rsdD3dXb6YI73tTXso3+JR784awGEKqUOAyuBISISDaCU+k4p9fcJheYDVZRSR4HlwMuS/pGax3LTB0uRmz7MBEoAW1U2ppKaUW76MAloq5QK595MiUl5lPmfMupDpv7Rh9eBr9L7+jnwhhnyZVdu+mHJ7+tMmfN9rY+M1TRNs3LWMKLXNE3TsqALvaZpmpXThV7TNM3K6UKvaZpm5XSh1zRNs3K60Guaplk5Xeg1TdOsnC70mqZpVu7/AACMXZj/uC2LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "huc8 = '02040202'\n",
    "fname = f'/opt/data/noaa/huc8-extracts/transformed/{huc8}.json'\n",
    "d = gpd.read_file(fname,driver='GeoJSON') ## neither fiona nor geopandas can read lists from geojson\n",
    "d.plot(facecolor='none')\n",
    "print(d.columns) ## comids column is missing!\n",
    "\n",
    "with open(fname) as fd:\n",
    "  js = json.load(fd)\n",
    "print(f'# of features => {len(js[\"features\"])}')\n",
    "print(f'huc8 => {js[\"features\"][0][\"properties\"][\"huc8\"]}')\n",
    "print(f'first 5 comids => {js[\"features\"][0][\"properties\"][\"comids\"][:5]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f5a934",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}