detect_cont_shelf.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b9b94de9-a606-4bbf-a882-03faa0e423db",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.features import shapes\n",
    "from shapely.geometry import shape, MultiPolygon, Polygon, box\n",
    "import geopandas as gpd\n",
    "import ocsmesh\n",
    "from osgeo import gdal\n",
    "import rasterio\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8e8477d1-f13b-4183-9f4c-64b0197c0914",
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_interiors(poly):\n",
    "    \"\"\"\n",
    "    Close polygon holes by limitation to the exterior ring.\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    poly: shapely.geometry.Polygon\n",
    "        Input shapely Polygon\n",
    "\n",
    "    Returns\n",
    "    ---------\n",
    "    Polygon without any interior holes\n",
    "    \"\"\"\n",
    "    if poly.interiors:\n",
    "        return Polygon(list(poly.exterior.coords))\n",
    "    else:\n",
    "        return poly\n",
    "\n",
    "\n",
    "def pop_largest(gs):\n",
    "    \"\"\"\n",
    "    Pop the largest polygon off of a GeoSeries\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    gs: geopandas.GeoSeries\n",
    "        Geoseries of Polygon or MultiPolygon objects\n",
    "\n",
    "    Returns\n",
    "    ---------\n",
    "    Largest Polygon in a Geoseries\n",
    "    \"\"\"\n",
    "    geoms = [g.area for g in gs]\n",
    "    return geoms.pop(geoms.index(max(geoms)))\n",
    "\n",
    "\n",
    "def close_holes(geom):\n",
    "    \"\"\"\n",
    "    Remove holes in a polygon geometry\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    gseries: geopandas.GeoSeries\n",
    "        Geoseries of Polygon or MultiPolygon objects\n",
    "\n",
    "    Returns\n",
    "    ---------\n",
    "    Largest Polygon in a Geoseries\n",
    "    \"\"\"\n",
    "    if isinstance(geom, MultiPolygon):\n",
    "        ser = gpd.GeoSeries([remove_interiors(g) for g in geom])\n",
    "        big = pop_largest(ser)\n",
    "        outers = ser.loc[~ser.within(big)].tolist()\n",
    "        if outers:\n",
    "            return MultiPolygon([big] + outers)\n",
    "        return Polygon(big)\n",
    "    if isinstance(geom, Polygon):\n",
    "        return remove_interiors(geom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "49952b69-89b1-4e59-916b-cdc5b01a06f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eddy(path,hmax,outdir):\n",
    "    r = ocsmesh.Raster(path)\n",
    "\n",
    "    shps = shapes(r.get_values(band=True), transform=r.src.transform)\n",
    "    polys = []\n",
    "    sids = []  # For contours this will be elev levels\n",
    "    for shp, sid in shps:\n",
    "        if sid<=hmax:\n",
    "            polys.append(shape(shp))\n",
    "            sids.append(sid)\n",
    "\n",
    "    gdf = gpd.GeoDataFrame(\n",
    "        geometry=polys,\n",
    "        data={'sid': sids},\n",
    "        crs=r.crs\n",
    "    ).dropna().dissolve()\n",
    "\n",
    "    gdf.to_file(outdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8812112d-9f9c-4fe3-a2cc-108f11550c69",
   "metadata": {},
   "outputs": [],
   "source": [
    "def land_mask(path,hmin,outdir):\n",
    "    r = ocsmesh.Raster(path)\n",
    "\n",
    "    shps = shapes(r.get_values(band=True), transform=r.src.transform)\n",
    "    polys = []\n",
    "    sids = []  # For contours this will be elev levels\n",
    "    for shp, sid in shps:\n",
    "        if sid>=hmin:\n",
    "            polys.append(shape(shp))\n",
    "            sids.append(sid)\n",
    "\n",
    "    gdf = gpd.GeoDataFrame(\n",
    "        geometry=polys,\n",
    "        data={'sid': sids},\n",
    "        crs=r.crs\n",
    "    ).dropna().dissolve()\n",
    "\n",
    "    gdf.to_file(outdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "26060b1a-6929-4802-8193-700c9e96716b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_survey_shapes(path,hmax,hmin,outdir):\n",
    "    r = ocsmesh.Raster(path)\n",
    "\n",
    "    shps = shapes(r.get_values(band=True), transform=r.src.transform)\n",
    "    polys = []\n",
    "    sids = []  # For contours this will be elev levels\n",
    "    for shp, sid in shps:\n",
    "        if sid>=hmin and sid<=hmax:\n",
    "            polys.append(shape(shp))\n",
    "            sids.append(sid)\n",
    "\n",
    "    gdf = gpd.GeoDataFrame(\n",
    "        geometry=polys,\n",
    "        data={'sid': sids},\n",
    "        crs=r.crs\n",
    "    ).dissolve(by='sid').dropna()#.dropna()#.dissolve()\n",
    "\n",
    "    gdf.to_file(outdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b01a6988-868f-4c6c-8e23-f55ed6c666fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_slope(DEM):\n",
    "    gdal.DEMProcessing('slope.tif', DEM, 'slope')\n",
    "    with rasterio.open('slope.tif') as dataset:\n",
    "        slope=dataset.read(1)\n",
    "    return slope\n",
    "    gdal.DEMProcessing('slope.tif', DEM, 'slope',scale=111120)\n",
    "    with rasterio.open('slope.tif') as dataset:\n",
    "        slope=dataset.read(1)\n",
    "    return slope"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b78e8c-f9a8-4224-881b-cf8ef3f07739",
   "metadata": {},
   "source": [
    "## Find non-eddy zones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bfac6c1b-738d-4262-9341-5ef204fc48d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "hmax=0\n",
    "hmin=-1000\n",
    "interv=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0b64f288-dd03-411b-8986-f29036bb4aff",
   "metadata": {},
   "outputs": [],
   "source": [
    "#rast_path =r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\dem\\smaller_domain\\GEBCO_28_Oct_2024_cea7192fe05d/gebco_2024_n26.0_s22.0_w-82.0_e-78.0.tif\"\n",
    "rast_path =r\"C:/Users/Felicio.Cassalho/Work/Python_Development/Mesh/OCSMesh_Contributions/ocean_mesh/dem/small_domain/GEBCO_28_Oct_2024_994bac2e22a9/gebco_2024_n28.0_s15.0_w-90.0_e-73.0.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6002a884-7d76-4d06-85c6-0df747392fc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate_slope(rast_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e5ef15eb-1995-4ca9-95ab-0304b1e374e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_path=r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\slope.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "032f1511-9555-4b4d-9a97-7ca4357f3e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "50fc49fa-a13e-41a5-981b-c714d7581f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "get_survey_shapes(path=rast_path,hmax=hmax,hmin=hmin,outdir=outdir+r\"/greater_1000_ind.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "ff172180-5492-49c3-996c-77ad931117f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/greater_1000_ind.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "557e2ec1-431a-469f-8310-f198fadac010",
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in range(hmin,hmax,interv):\n",
    "    y=x+interv\n",
    "    elev.loc[(elev['sid'] >= x) & (elev['sid'] < y), 'range'] = f'{x}-{y}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "23dae2d2-075c-474d-8fa0-f561ffef3696",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = elev.dissolve(by='range').dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "e5655ba5-95c3-422f-9392-e2de75152975",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = elev.sort_values(by=['sid'],ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "69d6f7af-9877-4168-abb9-6f0d6123351d",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.reset_index(drop=False, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "3474fb3a-45ea-4b4c-a9f6-c0cb594ef07c",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.to_file(outdir+r\"/greater_1000_by10.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "7bf034ed-f78e-46b8-a19c-8efac2b08376",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterstats import zonal_stats\n",
    "\n",
    "# Read the raster\n",
    "with rasterio.open(slope_path) as raster:\n",
    "    affine = raster.transform\n",
    "    array = raster.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "2a37da1d-de9d-4a20-acdd-091a9e6baa4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = elev#.explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "0b0b50a5-5b83-45d3-9f68-d807821f3ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = zonal_stats(polygons, slope_path, affine=affine, stats=[\"median\",\"mean\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "d72eb2bf-5115-4661-b7f5-1c74be53822a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate percentage\n",
    "list_median,list_mean=[],[]\n",
    "for stat in stats:\n",
    "    list_median.append(stat['median'])\n",
    "    list_mean.append(stat['mean'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "453fd860-f10f-4199-bfe1-8deff32a319d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# polygons[\"median_slope\"]=list_median\n",
    "polygons[\"mean_slope\"]=list_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b03675ec-da81-4567-8bd2-edeb1d3913a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "1cf93ac0-9ee8-4c5b-a5f3-093bf24640b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='mean_slope', ylabel='sid'>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGxCAYAAACtEoj/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+g0lEQVR4nO3dfXTU5Zn/8c8QTYRIQmAgIRBILELBiEAiGiKkYAm1wErtUYEuEkrdQ3kqDXRPg60iFbEL+Os2LVS7FqWnC90uWrtVNKyiCKKEADaoK2JIgoYQCCFDgiYhzO8PnDGTzHPmed6vc3IOmbln8p2NZ3P1vq7v5zaYzWazAAAA4LYewb4AAACAcEMBBQAA4CEKKAAAAA9RQAEAAHiIAgoAAMBDFFAAAAAeooACAADwEAUUAACAh64J9gVEoitXrqimpka9e/eWwWAI9uUAAAA3mM1mXbx4UampqerRw/keEwWUH9TU1CgtLS3YlwEAALxw6tQpDR482OkaCig/6N27t6Srv4CEhIQgXw0AAHCHyWRSWlqa9e+4MxRQfmBp2yUkJFBAAQAQZtwZv2GIHAAAwEMUUAAAAB6igAIAAPAQBRQAAICHKKAAAAA8RAEFAADgIQooJzZv3qyMjAxdd911ysrK0ltvvRXsSwIAACGAAsqBP//5z1qxYoUeeughHTlyRBMnTtRdd92l6urqYF8aAAAIMoPZbDYH+yJC0W233aZx48Zpy5Yt1sdGjhypWbNmaf369U5fazKZlJiYqMbGRoI0AQAIE578/WYHyo7W1laVlZUpPz/f5vH8/Hy9/fbbXda3tLTIZDLZfAEAgMhFAWXHuXPn1N7eruTkZJvHk5OTVVtb22X9+vXrlZiYaP3y50HCFWebtOejOp081+y3nwEAAJzjLDwnOp+FYzab7Z6PU1RUpMLCQuv3lsMIfenCpVYt335Uez8+a31s0o39VTxnrBJ7XevTnwUAAJxjB8oOo9GomJiYLrtNdXV1XXalJCkuLs56cLC/DhBevv2o9p84Z/PY/hPntGz7EZ//LAAA4BwFlB2xsbHKysrS7t27bR7fvXu3JkyYEPDrqTjbpL0fn1V7p3n/drNZez8+63E7jzYgAADdQwvPgcLCQs2bN0/Z2dnKycnR008/rerqai1atCjg11J1/pLT5yvrm5VhjHf5PrQBAQDwDQooB+6//37V19dr7dq1On36tDIzM/Xyyy9r6NChAb+WoX17OX0+vZ/r4kly3gbctnC819cHAEC0oYXnxOLFi1VZWamWlhaVlZVp0qRJQbmOG/pfr0k39ldMpwH2GINBk27s79buk6/bgAAARDMKqDBRPGescocZbR7LHWZU8Zyxbr3+3ZPnnT5fWU8BBQCAu2jhhYnEXtdq28LxOnmuWZX1zUrvF+/13JM97rYBAQAABVTYyTC6VzhZ2Jt76ijGYFDuMKNH7wkAQLSjhRfBHM09dZQ7zKiV+cOJNQAAwAPsQEUwV/EHD88YqTc+Oqe7f7vf+hixBgAAuMYOVARzFX+w61gt6eYAAHiBAiqCOYs/yB6apNLKBmINAADwAgVUhHMUf7BgQrrT171f0+jHqwIAILwxAxXhHMUfVJxtcvq6596u1IzRqQG6SgAAwgs7UFEiwxivySMGWOMKbuh/vbKHJjlcX1rZQBsPAAAHKKCimKs23jsV9cQbAABgBy28KDYyNcHp80XPl1v/TbwBAABfYQcqijm6S88e4g0AAPgKBVSUs3eXnj2+ijeoONtEWxAAEPZo4UW5znfpnWn8Qj/t0LrrrLK+2atz8+wdakxbEAAQrtiBgqSv7tIbn9HX6br0ft4dOmzvUGPaggCAcEUBBRvO0ssn3djfq90nR4cak3oOAAhXFFDowlF6efGcsV69n6tDjSvrfVNAMV8FAAgUZqDQhaP0cm+5OtTY27agBfNVAIBAYwcKDnVOL/eWP9qCHTFfBQAINAooBISv24IWwZ6vom0IANGJFh4CwtdtQQt35qt88XM6o20IANGNHSgElK/aghb+nq9yhLYhAEQ3CiiENX/PV9kT7Lahv9GWBADXaOEh7BXPGatl24/YtNN8MV/lSLDahv5GWxIA3EcBhbDnr/kqR4LVNvQ3Z23JbQvHB+mqACA00cJDxPD1fJUjwWgb+luktyUBwNcooAAv+CuWIVgClRbfGfNWAMIVLTzAC4FuG/pboNuSzFsBCHfsQAHdEKi2ob8Fui1JDASAcEcBBUBS4NqS3sxb0eoDEGpo4QGQFLi2pCcxELT6AIQqdqAA2PB3W9KTeStafQBCFQUUgIByd96KaAUAoYwCCkDAWGaZVuUPdzlvFaxoBQBwBzNQAPzO0SzT35bmqr651e68VaQmvgOIDBRQAPzO0SyTJIfHxFhafftPnLNp48UYDModZgz76AgA4Y0CCoBfWWaZOus4y9S5GKo426Sq85e0Kn+4JAXsoGgAcBcFFAC/8kVsgbNWHwAEAwUUAL/yRWyB5LjVBwDBwF14APwqELEFJJUDCDR2oAD4XfGcsVq2/YjTWSZPWn0WJJUDCBYKKAB+584xMd7EFjhLKqflB8CfaOEBCBhnx8S42+qzIKkcQDBRQAEIGcVzxrpMKLfwJKmcGSkAvkYLD0DIcKfVZ+FOy48ZKQD+wg4UgJDjrNVn4U7Lz9mMFAB0BwUUgLDlrOXnakZq7/GztPUAeI0WHoCw5azld/hUg9PXPvCHg9Z/09YD4Cl2oACEPXstP1czUh3R1gPgKQooABHJ0YyUPcGKPuDuQCB80cIDELHsJaA7Yy/t3B+4OxAIfxRQACJW5xmpGIPBZvapM3tp5/5AgjoQ/iigAES8DONXw+WTbuyv/SfO2dydF2MwKHeYMSC7T5a7Azvr2EYMxHUA6B5moABEFU/Szv3BkwT1jpiXAkILO1AAooonaef+4OmhycxLAaGJHSgAUcmdtHN/8PTQZNLUgdAUMQVUZWWlFi5cqIyMDPXs2VNf+9rX9Mgjj6i1tdVmXXV1tWbOnKn4+HgZjUYtX768y5ry8nLl5eWpZ8+eGjRokNauXStzpzRjANHJF600d9uIrtLUaecBwRMxLbz/+7//05UrV/TUU09p2LBhOnbsmB588EE1Nzdr48aNkqT29nZNnz5d/fv31759+1RfX6/58+fLbDaruLhYkmQymTR16lRNnjxZpaWlOn78uAoKChQfH6+VK1cG8yMCCCJfttLcbSO6My/FwDkQHAZzBG+tbNiwQVu2bFFFRYUkadeuXZoxY4ZOnTql1NRUSdKOHTtUUFCguro6JSQkaMuWLSoqKtKZM2cUFxcnSXriiSdUXFysTz/9VAY3QvlMJpMSExPV2NiohIQE/31AAAHzwDMHHd6956/ogYqzTZqy6U2Hz+9Z9Q0KKMCHPPn7HTEtPHsaGxvVt29f6/cHDhxQZmamtXiSpGnTpqmlpUVlZWXWNXl5edbiybKmpqZGlZWVdn9OS0uLTCaTzReAyBGsVpqn81IAAidiC6hPPvlExcXFWrRokfWx2tpaJScn26xLSkpSbGysamtrHa6xfG9Z09n69euVmJho/UpLS/PlRwEQZO600vwVMxDs2AUA9oX8DNSaNWv06KOPOl1TWlqq7Oxs6/c1NTX61re+pXvvvVc/+MEPbNbaa8GZzWabxzuvsXQ5HbXvioqKVFhYaP3eZDJRRAERxFX0wOY9J1Ra2WD93pcxA8GOXQBgX8gXUEuXLtXs2bOdrklPT7f+u6amRpMnT1ZOTo6efvppm3UpKSl69913bR5raGhQW1ubdZcpJSWly05TXV2dJHXZmbKIi4uzafkBiCyWVpq9GaiEntfocNUFm/X+OJalY5o6gOAL+QLKaDTKaDS6Xijps88+0+TJk5WVlaWtW7eqRw/bDmVOTo7WrVun06dPa+DAgZKkkpISxcXFKSsry7pm9erVam1tVWxsrHVNamqqTaEGILrYO5h47JA+OlTV0GWtu8eyVJxtUtX5S+wqAWEoYu7Cq6mpUV5enoYMGaJt27YpJibG+lxKSoqkqzEGY8aMUXJysjZs2KDz58+roKBAs2bNssYYNDY2asSIEZoyZYpWr16tjz/+WAUFBXr44YfdjjHgLjwgcnVspVXWN2vB1lKHa7cuuFWTRwzo8jjp4kBo8uTvd8jvQLmrpKREJ06c0IkTJzR48GCb5yw1YkxMjF566SUtXrxYubm56tmzp+bOnWvNiZKkxMRE7d69W0uWLFF2draSkpJUWFhoM+MEIHp1bKW5+t+fnY9lsXCWLu6vSAQAvhUxO1ChhB0oIHp4mg/l62wn2oCA70TlDhQABIO92ShnMQO+ShenDQgEFwUUAHSDpzEDriIRHLX9OqMNCARXxAZpAkAgZRjjNXnEAJe7R75IF+eQYSD4KKAAIMC6my7uqg34TkW9X1LRAXyFFh4ABFh308VdtQGLni+3/pu5KMA/2IECgCBxt+3XmaM2oD2WuSgAvkUBBQBhyF4b0B7LXNRbHe7WA9B95ED5ATlQAALF0gY80/iFftqhdWcP7TzAOU/+frMDBQBhzNIGHJ/R1+Va2nmA71BAAUAEcGcuipgDwHcooAAgQrg7F1VZ/1UBVXG2icgDwAvEGABAhLDEI+w9flYP/OGgw3Xp/eI5CgboJnagACDCTBre32XaubOjYAC4RgEFABHIXjtv7JA+uj97sPYer+MoGKCbaOEBQATqmHZ+rKZR296uVGllgw5VNbh8bWV9s8fhnkC0YQcKACJYhjFefyn9VIerLrj9mvR+FE+AK+xAAUAEqzjbZDMo7kyMwaDcYUaXu08VZ5tUdf6Sx2f4AZGEAgoAIljV+Utur80dZlTxnLEOn+fOPeArFFAAEMGG9u3l9Pk/Lhyvy1fMbu0mObtzb9vC8d2+ViCcUEABQASzJJTvP3HO5q47S7tu4o393XofR63Ajnfu0c5DNGGIHAAinL1IA1ftus5ctQI7ppt7giR0hCt2oAAgwnWMNKisb/Zq+NtVK9DTO/eYp0K4YwcKAKJEhjFek0cM8KrV5uiw4o7p5p4gCR3hjgIKAOAWX7QCpa/mqXyVhE4bEMFACw8A4BZftAIl9+ap3Hlf2oAIJnagAAAe6U4rUPLdPBVtQAQTBRQAIKB8MU/l6zYg4CkKKABAwHV3nsrXsQrMUcFTzEABAAKuu/NUvmoDMkcFb7EDBQAIGm/nqXwVq8AcFbxFAQUACEvdbQN6OkdFmw8d0cIDAISl7rYB3Y1ToM0He9iBAgCENW/bgO7OUdHmgz0UUACAqOTOHJUv4hJo/UUmWngAgKhVPGeslm0/YtOe6zhH1Z3UdFp/kY0CCgAQtVzNUXUnLsFZ62/bwvHdu3AEHS08AEDUczRH5W1cAknpkY8CCgAAJ7yJS/B1UjpCDy08AACc8CYuwVXrr198rC8vEUHADhQAAG7wJC7B0vpzZOOrx315aQgCCigAAHygc1zByvzhDtcyBxX+aOEBANANjuIK7rt1sNPXOYtAQOijgAIAoBscxRV83nbZ6eucRSAg9NHCAwDAS87iCkorG3Tr0CSPIxAQHiigAADwkqu4goIJ6R5HICA80MIDAMALFy61avOeE07XjBqUqG23pHoUgYDwQAEFAIAXlm8/qsNVF+w+F2MwKHeY0VosZRgpnCINBRQAAG6qONukqvOXFGOQzV13nY0b0oc2XYSjgAIAwAV7UQXOLJ4yTIm9rvXzVSGYGCIHAMAFe1EFzhBREPnYgQIAwAlLVIE7Os8+WVp+DI9HHgooAACccBVV0JElosBROnnxnLG09iIEBRQAAE4M7dvL6fN/XDhel6+YbXaZHnjmoN108mXbj2jbwvF+u1YEDgUUAABO3ND/ek26sb/2nzhnkzhuaddNvLG/zXpHLb92s9l6iDDtvPDHEDkAAC4UzxnrdqK4q5ZfZX1zt66l4myT9nxUp5Pnuvc+6B52oAAAcCGx17XatnC8W4nirlp+3t6hx1xVaGEHCgAAN2UY4zV5xACnLThLy8/Xhwjbi1KwzFUh8CKygGppadGYMWNkMBh09OhRm+eqq6s1c+ZMxcfHy2g0avny5WptbbVZU15erry8PPXs2VODBg3S2rVrZe500jYAAI540vJzh2Wuqr3T36KOc1UIrIhs4f3rv/6rUlNT9d5779k83t7erunTp6t///7at2+f6uvrNX/+fJnNZhUXF0uSTCaTpk6dqsmTJ6u0tFTHjx9XQUGB4uPjtXLlymB8HABAmPGk5ecOV3NVy7Yf1p8W3k4rL4AiroDatWuXSkpKtHPnTu3atcvmuZKSEn3wwQc6deqUUlNTJUmbNm1SQUGB1q1bp4SEBP3pT3/SF198oWeffVZxcXHKzMzU8ePH9eSTT6qwsFCGTluyAAA44qtDhF3NVX1QYyIiIcAiqoV35swZPfjgg/rjH/+oXr26/sd24MABZWZmWosnSZo2bZpaWlpUVlZmXZOXl6e4uDibNTU1NaqsrLT7c1taWmQymWy+AADwFctclaM/2lfMopUXYBFTQJnNZhUUFGjRokXKzs62u6a2tlbJyck2jyUlJSk2Nla1tbUO11i+t6zpbP369UpMTLR+paWldffjAABgo3jOWI1KTXC65p2KeiIOAiTkC6g1a9bIYDA4/Tp06JCKi4tlMplUVFTk9P3steDMZrPN453XWAbIHbXvioqK1NjYaP06deqUpx8TAACnEntdq1+7GEIver5cC7aWavLGN/TAMwfVeKktQFcXfUJ+Bmrp0qWaPXu20zXp6el67LHH9M4779i03iQpOztb3/ve9/Tcc88pJSVF7777rs3zDQ0Namtrs+4ypaSkdNlpqqurk6QuO1MWcXFxXX4uAAC+5igV3R6OjvGvkC+gjEajjEajy3W//vWv9dhjj1m/r6mp0bRp0/TnP/9Zt912myQpJydH69at0+nTpzVw4EBJVwfL4+LilJWVZV2zevVqtba2KjY21romNTVV6enpPv50AAB4pnjOWC3bfsTucTEduXN0TMXZJlWdv9TtuwSjkcEcoQFHlZWVysjI0JEjRzRmzBhJV2MMxowZo+TkZG3YsEHnz59XQUGBZs2aZY0xaGxs1IgRIzRlyhStXr1aH3/8sQoKCvTwww+7HWNgMpmUmJioxsZGJSQ471cDAOANS0TCmcYv9NPnyx2u27rgVk0eMcDmMVLN7fPk73fIz0D5UkxMjF566SVdd911ys3N1X333adZs2Zp48aN1jWJiYnavXu3Pv30U2VnZ2vx4sUqLCxUYWFhEK8cAABbllT08Rl9na6zd3QMqebdF/ItPG+lp6fbTQ8fMmSI/v73vzt97c0336y9e/f669IAAPAZR3NRMQaDcocZu7TmLKnmnbnT8sNXomoHCgCASOTu0TEXLrVq+Q7nu0yV9UQguCNid6AAAIgW7h4ds3z7UX1Q4zzs2V7LD11RQAEAECGcHR3jqHVn0UPSHTf2p33nJgooAADsiLRb/F0dSDwqNaFLyw+OUUABANBBpN7i7+pA4uK548L68wUaQ+QAAHQQqbf4W+7Wi+l0LFmMwaBJtO48RgEFAMCXLHNCnY9J6XiLv6vXh/Jhvu7erQfXaOEBAPAlV3NClfX2M5LCpe3n7t16cI0dKAAAvuRqTsjRLf7h1vazpJhTPHmPAgoAgC95MyfU3bYfwhMFFAAAHXg6J+Sq7fe39z6jiIpABrO9A+PQLZ6c5gwACE3uzglVnG3SlE1vuny/UJyJgi1P/n6zAwUAgB3uzgk5avt1FsozUfAcBRQAAN1kr+3XmWUm6i0nx6kgfNDC8wNaeAAQnU6ea9bf3vtM/2/3x07X0c4LTbTwAAAIggxjvGaOTnW5jnZe+KOAAgDAh9yZifJFxEGop55HOpLIAQDwseI5Y7Vs+xGbZHJ7HCWbOxMuqeeRjh0oAAB8zHJkyrbvj3e6zlGyuTPhlnoeqSigAADwk0nD+3ucbO4MqeehgwIKAAA/8jTZvLOOs07uHHaMwGAGCgAAP7K089xNNrewN+t0a3qS09d40xKEd9wuoAoLC91+0yeffNKriwEAIFJlGN0rnCzszTodrrqgpF7XyvT5ZZs2XozBoNxhRo9bgvCe2wXUkSO2w2llZWVqb2/XiBEjJEnHjx9XTEyMsrKyfHuFAABEsIqzTao6f8lmZ8oy69RZu9mshkttunVokkqrGqyPe9IShG+4XUDt2bPH+u8nn3xSvXv31nPPPaekpKvbiQ0NDVqwYIEmTpzo+6sEACDCOIsjcDXrtHjKMKX3i/eoJQjf8uool0GDBqmkpEQ33XSTzePHjh1Tfn6+ampqfHaB4YijXAAArjzwzEHtP3HObituzT+N0pRNbzp87Z5V36Bo8gO/H+ViMpl05syZLo/X1dXp4sWL3rwlAABBE+hUb1dxBIYvYw58FX8A3/PqLrzvfOc7WrBggTZt2qTbb79dkvTOO+/oJz/5ie655x6fXiAAAP4SrFRvd+II7KWZM+sUOrwqoH73u99p1apV+ud//me1tbVdfaNrrtHChQu1YcMGn14gAAD+4izVe9tC5yni3TG0by+nz6f3i/c6/gCB4dUMlEVzc7M++eQTmc1mDRs2TPHx/GIlZqAAIBxUnG0K6pyRsxkofxZvcMzvM1AW8fHxGj16tG655RaKJwBAWHn35Hmnz/s71bu7CeUILrdbePfcc4+effZZJSQkuJxzev7557t9YQAA+IO9uSd7/J3qTYsuvLldQCUmJsrw5d0AiYmJfrsgAAD8yd7cU0eBTvX2NKEcocHtAmrr1q3Wf2/evFlXrlyxtu0qKyv117/+VSNHjtS0adN8f5UAAPiAo4TvjkYO7K1V+cMDdEUIV17NQN1999364x//KEm6cOGCbr/9dm3atEmzZs3Sli1bfHqBAAD4iqv4AEk6VmPSP/12vx545qAaL7UF4KoQjrwqoA4fPmw9suW///u/lZycrKqqKm3btk2//vWvfXqBAAD4iqv4gI4scQaAPV4VUJcuXVLv3r0lSSUlJbrnnnvUo0cP3X777aqqqvLpBQIA4Cs39L/ebsK3PZZU8EClkyO8eFVADRs2TH/961916tQpvfrqq8rPz5d09SgXco8AAIHg7fEr9uIDnPF3nAHCk1dJ5A8//LDmzp2rH//4x7rzzjuVk5Mj6epu1Nix5FcAAPynu8evdI4PiDEY9MAfDjpc7+84A4Qnr5PIa2trdfr0ad1yyy3q0ePqRtbBgweVkJCgr3/96z69yHBDEjkA+I8/ErxJBYfk2d9vr3agJCklJUUpKSk2j40fz39kAAD/cRRD0HFeydNMpYqzTbo/e7A+b72s0qoG6+OkgsMZrwsoAAACzVUMQWW9+wWUvVbgrelJmj8hXTelJhJuCae6dRYeAACB5CqGwJN5JXuJ5IerLui/Sj+leIJLFFAAgLDhKIYgxmDQpBv7u134WFqB7Z3GgLsbXeDtnYEIP7TwAABhpXjOWC3bfsSm9ebpvJIvW4FS9+8MRPihgAIAhJXOMQTp/Tw/jNeXrUDJfjvQkmTOXXyRiRYeACAsZRjjNXnEAK/mlXzVCpT81w5EaKOAAgBEJXuJ5N5EF7jTDnSEmanwRQsPABCVfNEKlLxrBzIzFf7YgQIARLXutAIl79qBzmamEB4ooAAA6AZLkvm4IX1sHnfUDmRmKjLQwgMAwAveJpn7OkIBwcEOFAAAXvA2ydzXEQoIDgooAAA81J02nC8jFBA8FFAAgKjUnQiB7kQXSL6LUEDwMAMFAIgqvogQ6G4bzlcRCggedqAAAFHFFxECvmrDdTdCAcETcQXUSy+9pNtuu009e/aU0WjUPffcY/N8dXW1Zs6cqfj4eBmNRi1fvlytra02a8rLy5WXl6eePXtq0KBBWrt2rcyd+twAgPDjywgB2nDRLaJaeDt37tSDDz6oxx9/XFOmTJHZbFZ5ebn1+fb2dk2fPl39+/fXvn37VF9fr/nz58tsNqu4uFiSZDKZNHXqVE2ePFmlpaU6fvy4CgoKFB8fr5UrVwbrowEAuunCpVYt3+F8l8mTCAHacNHNYI6QrZXLly8rPT1djz76qBYuXGh3za5duzRjxgydOnVKqampkqQdO3aooKBAdXV1SkhI0JYtW1RUVKQzZ84oLi5OkvTEE0+ouLhYn376qQydtmvtMZlMSkxMVGNjoxISEnz3IQEAXnvgmYPad+Ksrjj5q7dn1TcogqKYJ3+/I6aFd/jwYX322Wfq0aOHxo4dq4EDB+quu+7S+++/b11z4MABZWZmWosnSZo2bZpaWlpUVlZmXZOXl2ctnixrampqVFlZafdnt7S0yGQy2XwBAEKHpXXnqHjqIQUlQoDDhMNXxBRQFRUVkqQ1a9boZz/7mf7+978rKSlJeXl5On/+vCSptrZWycnJNq9LSkpSbGysamtrHa6xfG9Z09n69euVmJho/UpLS/PpZwMAdI+r2IFRqQkBnV26cKlVDzxzUFM2vakFW0s1eeMbeuCZg2q81Bawa0D3hHwBtWbNGhkMBqdfhw4d0pUrVyRJDz30kL773e8qKytLW7dulcFg0F/+8hfr+9lrwZnNZpvHO6+xdDkdte+KiorU2Nho/Tp16lS3PzcAwHdcxQ4Uzx3ndoSBL3CYcPgL+SHypUuXavbs2U7XpKen6+LFi5KkUaNGWR+Pi4vTDTfcoOrqaklSSkqK3n33XZvXNjQ0qK2tzbrLlJKS0mWnqa6uTpK67Ex1/DkdW34AgNBiiR3Yf+KczR14MQaDcocZA9q6s7QTO+t4JyBzWKEv5HegjEajvv71rzv9uu6665SVlaW4uDh99NFH1te2tbWpsrJSQ4cOlSTl5OTo2LFjOn36tHVNSUmJ4uLilJWVZV2zd+9em2iDkpISpaamKj09PTAfGgDgc6ESO9DdFHNvMW/lWyG/A+WuhIQELVq0SI888ojS0tI0dOhQbdiwQZJ07733SpLy8/M1atQozZs3Txs2bND58+e1atUqPfjgg9Zp+7lz5+rRRx9VQUGBVq9erY8//liPP/64Hn74YbfuwAMAhKZQiR0I9GHCvkheR1cRU0BJ0oYNG3TNNddo3rx5+vzzz3Xbbbfp9ddfV1JSkiQpJiZGL730khYvXqzc3Fz17NlTc+fO1caNG63vkZiYqN27d2vJkiXKzs5WUlKSCgsLVVhYGKyPBQDwoQxjcPOaAt1OdDZvtW3heJ/+rGgSMTlQoYQcKACAM42X2rRs+xG/7wpVnG3SlE1vOny+Y+5VxdkmVZ2/FNWBoJ78/Y6oHSgAAMJBoNqJ7sxbJfW6lhafF0J+iBwAgEjl78OE3Zm3IlLBOxRQAABEKMu8VUynm6BiDAZNurG/zF9GJ/jicOVoQwEFAECEsBdV4Cy+IViRCpGAGSgAAMKcq6gCR/NWgY5UiCTsQAEAEObcmWOyN2/lqsUXrXfjuYMCCgCAEOFNWrjlaBhv55hCJaE93NDCAwAgyLqTFu7OHJOznaRQSWgPN+xAAQAQZN2JEvDVHJO/IxUiDQUUAAA+4O1hvd1twTHHFBy08AAA6IbuHtbb3RacdHWOqfPRMMwx+RcFFAAA3dDdw3p90YJjjinwaOEBAOCl7rbfJN+24JhjChwKKAAAvPTuyfNOn3c3yZsogfBDCw8AAA/Zm3uyx9074GjBhR8KKAAAPGRv7qmjGINBucOMHhdBGUYKp3BBAQUAgAcsc0/OjBzYW6vyhwfoihAMzEABAOABV7EDknSsxqR/+u1+PfDMQTVeagvAVSHQKKAAAPCAq9iBjtxNE0f4oYACAMADjmIH7PEkzgDhhQIKAAAP2YsdcKZjnIG3R74gtDBEDgCAhzrHDsQYDHrgDwcdrk/vF9/tI18QWtiBAgDAS5bk70nD+7tME3d25AvCDwUUAAA+4CxN3N0jX2jvhQ9aeAAA+ICzNPHDpxqcvvZYTaMeefF92nthhB0oAAB8yN6Bvq6iD7a9XUl7L8xQQAEA4GeOog9iDAZlD01SaWWDw/bejoPVtPRCEAUUAAAB4GhGasGEdKev++nz5Zq88Q1SzUMMM1AAAASAoxmpirNNbr3e0tLbtnC8n68U7mAHCgCAAOo8I+Vusjmp5qGFAgoAgCDzJNmcVPPQQAsPAIAg69jee6eiXkXPlztcS6p5aGAHCgCAEJFhjNec8UNINQ8DFFAAAIQYey29cUP6eJRqDv+igAIAIMQk9rpWv54zRremJ1kfK61q0LLtR/TBaZPT13ackYL/MAMFAEAIWr79qA5XXbB5bP+Jc/q87bLT16X3i3f6PHyDHSgAAEKMszZdaWWDbh2a5HRGCv5HAQUAgJf8FSNQdf6S0+cLJqTbTTUvnjPWp9cBx2jhAQDgIX/HCLg6fHjUoERtuyW1S6o5AocdKAAAPOTvGAFnhw93bNN1TjVH4FBAAQDggTc/qgtIjICjw4dp04UGWngAALjBXtvOnsr6Zp/sCDk6fBihgQIKAAA32Gvb2ePrGIEMI4VTKKKAAgDABUusgDMxBoPGDuljDbKk6IlsFFAAALjgKlZAkhJ6XqNDVQ1asLVUEof7RjqGyAEAcMFVrMDNgxJk+tw2IZzDfSMbBRQAAC44ixXIHpqk8s9MHO4bZSigAABwg6NYgW9npjh9HYf7RiZmoAAAcEPnWIG+vWK1qeS41r70odPXcbhvZKKAAgDAA5ZYgQeeOeg01iDGYFDuMCN340UoCigAANxQcbZJVecvKb1fvMxfzjc5Q2p4ZKOAAgDACXsJ5JmDEpy+5ol7btbs8UP8fWkIIobIAQBwwl4C+Qc1Jqevue2Gfv68JIQAdqAAAHDAUQL5lS8TC3pIutLhceaeogc7UAAAOOAqgXxUqm0rL9BzTxVnm7TnozqypoKAHSgAABxwlUBePHecpKtZT+n9Anfor725LI6OCayI2oE6fvy47r77bhmNRiUkJCg3N1d79uyxWVNdXa2ZM2cqPj5eRqNRy5cvV2trq82a8vJy5eXlqWfPnho0aJDWrl0rc6eEWQBA5HOWQD7pxv7WSIPJIwYEtG1nby6Lo2MCK6IKqOnTp+vy5ct6/fXXVVZWpjFjxmjGjBmqra2VJLW3t2v69Olqbm7Wvn37tGPHDu3cuVMrV660vofJZNLUqVOVmpqq0tJSFRcXa+PGjXryySeD9bEAAEHkKIE8WBEFlrksfx0dQ1vQPQZzhGytnDt3Tv3799fevXs1ceJESdLFixeVkJCg//3f/9Wdd96pXbt2acaMGTp16pRSU1MlSTt27FBBQYHq6uqUkJCgLVu2qKioSGfOnFFcXJwk6YknnlBxcbE+/fRTGTr9rxB7TCaTEhMT1djYqIQE57e6AgDCgyWBPJCtOnv2fFSnBVtLHT6/dcGtmjxigMfvS1vQs7/fEbMD1a9fP40cOVLbtm1Tc3OzLl++rKeeekrJycnKysqSJB04cECZmZnW4kmSpk2bppaWFpWVlVnX5OXlWYsny5qamhpVVlYG9DMBAEJHMFp19riay/L26Bjagp6JmCFyg8Gg3bt36+6771bv3r3Vo0cPJScn65VXXlGfPn0kSbW1tUpOTrZ5XVJSkmJjY61tvtraWqWnp9ussbymtrZWGRkZXX52S0uLWlparN+bTM7zQQAA8JZlLmv/iXM2bTxXEQodk9Q7r3EU19CxLRjswjHUhPwO1Jo1a2QwGJx+HTp0SGazWYsXL9aAAQP01ltv6eDBg7r77rs1Y8YMnT592vp+9lpwZrPZ5vHOayxdTkftu/Xr1ysxMdH6lZaW5ouPDgCAXZ7MZV241KoHnjmoKZve1IKtpZq88Q098MxBNV5qs65xFddQWc88VGchvwO1dOlSzZ492+ma9PR0vf766/r73/+uhoYGa99y8+bN2r17t5577jn99Kc/VUpKit59912b1zY0NKitrc26y5SSkmLdjbKoq6uTpC67VxZFRUUqLCy0fm8ymSiiAAB+k9jrWm1bON6tuSxnrbltC8dL8l9bMJKFfAFlNBplNBpdrrt06Wr13KOH7aZajx49dOXK1ZzYnJwcrVu3TqdPn9bAgQMlSSUlJYqLi7POSeXk5Gj16tVqbW1VbGysdU1qamqX1p5FXFyczcwUAACBYIlRcMTd1py3bcFoFvItPHfl5OQoKSlJ8+fP13vvvafjx4/rJz/5iU6ePKnp06dLkvLz8zVq1CjNmzdPR44c0WuvvaZVq1bpwQcftO5azZ07V3FxcSooKNCxY8f0wgsv6PHHH1dhYaFbd+ABABAqPGnNhVpcQ6gL+R0odxmNRr3yyit66KGHNGXKFLW1temmm27Siy++qFtuuUWSFBMTo5deekmLFy9Wbm6uevbsqblz52rjxo3W90lMTNTu3bu1ZMkSZWdnKykpSYWFhTYtOgAAwoGr1ly/+Fjrvz1pCyKCcqBCCTlQAIBQ8cAzB+228aSrOU+WOShEaQ4UAADRzFGC+Mr84Q5f44vk8mgVMS08AACikasE8fOXWp28+uocFK06z7EDBQBAGHOVIE5EgX9QQAEA4Cf+PpjXnYOFLREFMZ3uJI8xGDTpxv7sPnmJFh4AAD4WqIN53YkpyDDGq3jOWC3bfsTmeogo6B4KKAAAfMyd9G9fcLc9R0SB79HCAwDAh9xpq/mKp+25DGO8Jo8YQPHkAxRQAAB8yRczS4E+mJcE8eCghQcAiHq+nFkK9F1vtOeCgx0oAEDUcxUF4Ilg3fVGey6wKKAAAAHh71v6veWPmSXaapGPFh4AwK8CdUu/t9yNAvAEbbXIxw4UAMCvfNke8wd/zizRVotcFFAAAL8J5C393iKpG96ggAIAhe58TrgL9C393mJmCZ5iBgpAVAv1+ZxwFy4H2TKzBE+xAwUgqoX6fE64C7f2GDNLcBcFFICoFQ7zOZGA9phjtI7DFy08AFHLH7evoyvaY13ROg5/7EABiFrhMp8TKWiPfYXWcfijgAIQtcJtPgeRgdZx94VC65MWHoCoVjxnrJZtP2LTSmE+B/5E69h7odT6pIACENWYz0Gg0Tr2nrPW57aF4wN6LbTwAEDM5yBwaB17J9RanxRQAAAEWDRHO3g7vxRqqfa08AAACLBobB13d34p1Fqf7EABABAk0dQ67m50Q6i1PimgAACAX7k7v+SqvRdKrU9aeAAAwK9czS8dq2nUIy++77K9F0qtT3agAACAX7maX9r2dqVH7b1QaH1SQAEAAL9yNr+UPTRJpZUNIRNP4C4KKAAA4HeO5pcWTEh3+rrK+uaQOLqlM2agAACA3zmaX6o42+T0dZv3nFBpZYP1+2Ad3dIZO1AAACBgOs8vOWvvJfW6VoerLtg87kn0gT9RQAEAgKCy1967oX+8Gi61hexsFAUUAAAIKkt778UlucoclCBJ+rjOeWsv0Ee3dEYBBQAAQsKmkuP6sOaiW2sDfXRLZwyRAwCAgKs426Sq85dshsk7Bmk6EmMwKHeYMejH31BAAQCAgHF0qPB9tw526/XBOrqlMwooAAAQMI4OFf687bLT162/52bdfkO/oO88WTADBQAAAsLZocKllQ26dWiS3TiDSTf215zxQ0KmeJIooAAAgAO+TgB3dahwwYR0u2nlodCy64wWHgAAsOFoTqm7CeCuDhUeNShR225J7ZJWHorYgQIAADYczSl1NwHcWer4pBv7W4ulzmnloYgCCgAASLrastt+sMrhnJIvEsAdHSocim06Z2jhAQAQ5ey17ByprG/u1s6Qo0OFww0FFAAAUc5ey84RXyWAZxjDs3CyoIACACCKhVsCeKiggAIAIIq5ihawsMwpdT6CJVpRQAEAEMVcRQtYEsCTel2rZduP+DzaIFxxFx4AAFHMVbSAJQHcX9EG4YoCCgCAKLcyf7hGDuxt81jHaAFnR7D4ItogHNHCAwAgStmLL8gclKDHv3OzRg/uY33M1ZxUd6MNwhE7UAAARCl7bbkPay5q46vHbR5zNSflq2iDcEIBBQBAFHK3LWe56+7WoUkuj2CJJrTwAACIQq7acsdqGvXIi+/btPeSel2rhktt1u/D8QgWX6GAAgAgCrlqy217u1KHqy7YPGb6/LJuTU/S4snDoj4HKmxaeOvWrdOECRPUq1cv9enTx+6a6upqzZw5U/Hx8TIajVq+fLlaW1tt1pSXlysvL089e/bUoEGDtHbtWpk7bV+++eabysrK0nXXXacbbrhBv/vd7/z1sQAACApn8QXZQ5NUWtlgt71XWtkQ9cWTFEYFVGtrq+6991798Ic/tPt8e3u7pk+frubmZu3bt087duzQzp07tXLlSusak8mkqVOnKjU1VaWlpSouLtbGjRv15JNPWtecPHlS3/72tzVx4kQdOXJEq1ev1vLly7Vz506/f0YAADxVcbZJez6q8ypKoHjOWOUOM9o8ljvMqAUT0p2+rrI++mILOgubFt6jjz4qSXr22WftPl9SUqIPPvhAp06dUmpqqiRp06ZNKigo0Lp165SQkKA//elP+uKLL/Tss88qLi5OmZmZOn78uJ588kkVFhbKYDDod7/7nYYMGaJf/epXkqSRI0fq0KFD2rhxo7773e8G4qMCAOCSvQgCT5PBE3tdq20Lx+vkuWZV1jdbd5YqzjY5fV003nXXWdjsQLly4MABZWZmWosnSZo2bZpaWlpUVlZmXZOXl6e4uDibNTU1NaqsrLSuyc/Pt3nvadOm6dChQ2praxMAAKHAl8ngGcZ4TR4xwNqWc5VOHu3tOymCCqja2lolJyfbPJaUlKTY2FjV1tY6XGP53tWay5cv69w52/9QLVpaWmQymWy+AADwl0Akgztq70XrXXedBbWFt2bNGmtrzpHS0lJlZ2e79X6GTpWyJJnNZpvHO6+xDJB7uqaj9evXu/wcAAD4SiCSwR2193BVUAuopUuXavbs2U7XpKenu/VeKSkpevfdd20ea2hoUFtbm3VHKSUlxbrTZFFXVydJLtdcc8016tevn92fXVRUpMLCQuv3JpNJaWlpbl03AACeCmQyeIaRwsmeoBZQRqNRRqPR9UI35OTkaN26dTp9+rQGDhwo6epgeVxcnLKysqxrVq9erdbWVsXGxlrXpKamWgu1nJwc/c///I/Ne5eUlCg7O1vXXmt/KC8uLs5mrgoAAH+yzCjtP3HOpo0XYzAod5iRgicAwmYGqrq6WkePHlV1dbXa29t19OhRHT16VE1NV+8UyM/P16hRozRv3jwdOXJEr732mlatWqUHH3xQCQkJkqS5c+cqLi5OBQUFOnbsmF544QU9/vjj1jvwJGnRokWqqqpSYWGhPvzwQ/3hD3/QM888o1WrVgXtswMA0BkzSkFmDhPz5883S+rytWfPHuuaqqoq8/Tp0809e/Y09+3b17x06VLzF198YfM+//jHP8wTJ040x8XFmVNSUsxr1qwxX7lyxWbNG2+8YR47dqw5NjbWnJ6ebt6yZYtH19rY2GiWZG5sbPT68wIA4I6Ks03m1//vjLnibFOwLyXsefL322A2dxrhR7eZTCYlJiaqsbHRuvsFAABCmyd/v8OmhQcAABAqKKAAAAA8RAEFAADgIQooAAAAD1FAAQAAeIgCCgAAwEMUUAAAAB6igAIAAPAQBRQAAICHKKAAAAA8dE2wLyASWU7HMZlMQb4SAADgLsvfbXdOuaOA8oOLFy9KktLS0oJ8JQAAwFMXL15UYmKi0zUcJuwHV65cUU1NjXr37i2DwWB3jclkUlpamk6dOsWBw2GA31d44fcVXvh9hZdI/n2ZzWZdvHhRqamp6tHD+ZQTO1B+0KNHDw0ePNittQkJCRH3H2Ak4/cVXvh9hRd+X+ElUn9frnaeLBgiBwAA8BAFFAAAgIcooIIkLi5OjzzyiOLi4oJ9KXADv6/wwu8rvPD7Ci/8vq5iiBwAAMBD7EABAAB4iAIKAADAQxRQAAAAHqKACoLNmzcrIyND1113nbKysvTWW28F+5LgwN69ezVz5kylpqbKYDDor3/9a7AvCU6sX79et956q3r37q0BAwZo1qxZ+uijj4J9WXBgy5YtGj16tDVPKCcnR7t27Qr2ZcFN69evl8Fg0IoVK4J9KUFBARVgf/7zn7VixQo99NBDOnLkiCZOnKi77rpL1dXVwb402NHc3KxbbrlFv/nNb4J9KXDDm2++qSVLluidd97R7t27dfnyZeXn56u5uTnYlwY7Bg8erCeeeEKHDh3SoUOHNGXKFN199916//33g31pcKG0tFRPP/20Ro8eHexLCRruwguw2267TePGjdOWLVusj40cOVKzZs3S+vXrg3hlcMVgMOiFF17QrFmzgn0pcNPZs2c1YMAAvfnmm5o0aVKwLwdu6Nu3rzZs2KCFCxcG+1LgQFNTk8aNG6fNmzfrscce05gxY/SrX/0q2JcVcOxABVBra6vKysqUn59v83h+fr7efvvtIF0VELkaGxslXf2jjNDW3t6uHTt2qLm5WTk5OcG+HDixZMkSTZ8+Xd/85jeDfSlBxVl4AXTu3Dm1t7crOTnZ5vHk5GTV1tYG6aqAyGQ2m1VYWKg77rhDmZmZwb4cOFBeXq6cnBx98cUXuv766/XCCy9o1KhRwb4sOLBjxw4dPnxYpaWlwb6UoKOACgKDwWDzvdls7vIYgO5ZunSp/vGPf2jfvn3BvhQ4MWLECB09elQXLlzQzp07NX/+fL355psUUSHo1KlT+tGPfqSSkhJdd911wb6coKOACiCj0aiYmJguu011dXVddqUAeG/ZsmX629/+pr1792rw4MHBvhw4ERsbq2HDhkmSsrOzVVpaqn//93/XU089FeQrQ2dlZWWqq6tTVlaW9bH29nbt3btXv/nNb9TS0qKYmJggXmFgMQMVQLGxscrKytLu3bttHt+9e7cmTJgQpKsCIofZbNbSpUv1/PPP6/XXX1dGRkawLwkeMpvNamlpCfZlwI4777xT5eXlOnr0qPUrOztb3/ve93T06NGoKp4kdqACrrCwUPPmzVN2drZycnL09NNPq7q6WosWLQr2pcGOpqYmnThxwvr9yZMndfToUfXt21dDhgwJ4pXBniVLlug///M/9eKLL6p3797W3d7ExET17NkzyFeHzlavXq277rpLaWlpunjxonbs2KE33nhDr7zySrAvDXb07t27yzxhfHy8+vXrF5VzhhRQAXb//fervr5ea9eu1enTp5WZmamXX35ZQ4cODfalwY5Dhw5p8uTJ1u8LCwslSfPnz9ezzz4bpKuCI5Z4kG984xs2j2/dulUFBQWBvyA4debMGc2bN0+nT59WYmKiRo8erVdeeUVTp04N9qUBLpEDBQAA4CFmoAAAADxEAQUAAOAhCigAAAAPUUABAAB4iAIKAADAQxRQAAAAHqKAAgAA8BAFFAAAgIcooADAz9544w0ZDAZduHAh2JcCwEcooAAAADxEAQUAAOAhCigAIecb3/iGli1bphUrVigpKUnJycl6+umn1dzcrAULFqh379762te+pl27dllf88EHH+jb3/62rr/+eiUnJ2vevHk6d+6c9flXXnlFd9xxh/r06aN+/fppxowZ+uSTT6zPV1ZWymAw6Pnnn9fkyZPVq1cv3XLLLTpw4IBb11xVVaWZM2cqKSlJ8fHxuummm/Tyyy87XL9z507ddNNNiouLU3p6ujZt2mTzfHp6un7xi19o7ty5uv7665Wamqri4mKbNY2NjfqXf/kXDRgwQAkJCZoyZYree+89t64XQPdQQAEISc8995yMRqMOHjyoZcuW6Yc//KHuvfdeTZgwQYcPH9a0adM0b948Xbp0SadPn1ZeXp7GjBmjQ4cO6ZVXXtGZM2d03333Wd+vublZhYWFKi0t1WuvvaYePXroO9/5jq5cuWLzcx966CGtWrVKR48e1fDhwzVnzhxdvnzZ5fUuWbJELS0t2rt3r8rLy/XLX/5S119/vd21ZWVluu+++zR79myVl5drzZo1+vnPf65nn33WZt2GDRs0evRoHT58WEVFRfrxj3+s3bt3S5LMZrOmT5+u2tpavfzyyyorK9O4ceN055136vz58x7+XxuAx8wAEGLy8vLMd9xxh/X7y5cvm+Pj483z5s2zPnb69GmzJPOBAwfMP//5z835+fk273Hq1CmzJPNHH31k92fU1dWZJZnLy8vNZrPZfPLkSbMk83/8x39Y17z//vtmSeYPP/zQ5TXffPPN5jVr1th9bs+ePWZJ5oaGBrPZbDbPnTvXPHXqVJs1P/nJT8yjRo2yfj906FDzt771LZs1999/v/muu+4ym81m82uvvWZOSEgwf/HFFzZrvva1r5mfeuopl9cLoHvYgQIQkkaPHm39d0xMjPr166ebb77Z+lhycrIkqa6uTmVlZdqzZ4+uv/5669fXv/51SbK26T755BPNnTtXN9xwgxISEpSRkSFJqq6udvhzBw4caP0ZrixfvlyPPfaYcnNz9cgjj+gf//iHw7UffvihcnNzbR7Lzc3Vxx9/rPb2dutjOTk5NmtycnL04YcfSrq6i9XU1KR+/frZfO6TJ0/atCYB+Mc1wb4AALDn2muvtfneYDDYPGYwGCRJV65c0ZUrVzRz5kz98pe/7PI+liJo5syZSktL0+9//3ulpqbqypUryszMVGtrq8Of2/FnuPKDH/xA06ZN00svvaSSkhKtX79emzZt0rJly7qsNZvN1vfu+Jg7Ol7TwIED9cYbb3RZ06dPH7feC4D3KKAAhL1x48Zp586dSk9P1zXXdP1/a/X19frwww/11FNPaeLEiZKkffv2+fw60tLStGjRIi1atEhFRUX6/e9/b7eAGjVqVJef//bbb2v48OGKiYmxPvbOO+/YrHnnnXesO2vjxo1TbW2trrnmGqWnp/v8swBwjhYegLC3ZMkSnT9/XnPmzNHBgwdVUVGhkpISff/731d7e7uSkpLUr18/Pf300zpx4oRef/11FRYW+vQaVqxYoVdffVUnT57U4cOH9frrr2vkyJF2165cuVKvvfaafvGLX+j48eN67rnn9Jvf/EarVq2yWbd//37927/9m44fP67f/va3+stf/qIf/ehHkqRvfvObysnJ0axZs/Tqq6+qsrJSb7/9tn72s5/p0KFDPv1sALqigAIQ9lJTU7V//361t7dr2rRpyszM1I9+9CMlJiaqR48e6tGjh3bs2KGysjJlZmbqxz/+sTZs2ODTa2hvb9eSJUs0cuRIfetb39KIESO0efNmu2vHjRun//qv/9KOHTuUmZmphx9+WGvXrlVBQYHNupUrV6qsrExjx47VL37xC23atEnTpk2TdLWV9/LLL2vSpEn6/ve/r+HDh2v27NmqrKy0zocB8B+D2d3GOwAgYNLT07VixQqtWLEi2JcCwA52oAAAADxEAQUAbrjrrrts4gI6fj3++OPBvjwAAUYLDwDc8Nlnn+nzzz+3+1zfvn3Vt2/fAF8RgGCigAIAAPAQLTwAAAAPUUABAAB4iAIKAADAQxRQAAAAHqKAAgAA8BAFFAAAgIcooAAAADxEAQUAAOCh/w8vW+oJJ8BFEAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "polygons.plot.scatter(x='mean_slope', y='sid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "85a5fbfc-957f-4c78-8334-65054bcb8ca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# polygons = polygons.sort_values(by=['sid'],ascending=True)\n",
    "# polygons.reset_index(drop=False, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "f1a49bc4-7094-45e9-bed3-953914fe16b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/polygons.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "79ae6b75-6d52-4911-976c-8a93a770b6b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/polygons.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "fa277295-cb30-46a9-88b7-a16f438811be",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\AppData\\Local\\Temp\\1\\ipykernel_5336\\2290710304.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
      "\n",
      "  polygons[\"area\"] = polygons['geometry'].area\n"
     ]
    }
   ],
   "source": [
    "polygons[\"area\"] = polygons['geometry'].area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "506e95e6-ed8f-4e2a-9d9f-ef8db32bdbb5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='area', ylabel='sid'>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGwCAYAAACaW3CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2rUlEQVR4nO3df3TU1Z3/8dcQkiEJZBxJTDIQCAgibFjFBDEiRO0xlAKF0qMiXUpcyqlFSjWwVXCt1CNClx/taVrQWgvYesStiNstVUkVReSHAaPFHyVoEhIaYhKMCZAlCeTz/cNvpg7MjwDz64bn45w5h8y9n0/uZ/iQvLj3zntslmVZAgAAQFD1iPQAAAAAuiNCFgAAQAgQsgAAAEKAkAUAABAChCwAAIAQIGQBAACEACELAAAgBHpGegDdUUdHh2pqatSnTx/ZbLZIDwcAAHSBZVk6fvy4XC6XevS4+HkoQlYI1NTUKCMjI9LDAAAAF6C6ulr9+/e/6PMQskKgT58+kr78S0pKSorwaAAAQFc0NzcrIyPD/Xv8YhGyQqBziTApKYmQBQCAYYK11YeN7wAAACFAyAIAAAgBQhYAAEAIELIAAABCgJAFAAAQAoQsAACAECBkAQAAhAAhy4+1a9dq0KBB6tWrl7Kzs/XWW29FekgAAMAQhCwfnn/+ed1333166KGHVFpaqnHjxmnixImqqqqK9NAAAIABbJZlWZEeRDQaM2aMrrvuOq1bt8793PDhwzVt2jQtX77c77HNzc1yOBxqamqi4jsAAIYI9u9vZrK8aGtr0/79+5Wfn+/xfH5+vnbt2nVO/9bWVjU3N3s8AABAaJTXn9D2g3WqaDgZ6aH4xWcXetHQ0KAzZ84oNTXV4/nU1FTV1tae03/58uX66U9/Gq7hAQBwSfqipU0LnntPOw7Vu58bPzRFRXeNkiMhNoIj846ZLD/O/oBIy7K8fmjk4sWL1dTU5H5UV1eHa4gAAFwyFjz3nt7+pMHjubc/adAPnyuN0Ij8YybLi+TkZMXExJwza1VXV3fO7JYk2e122e32cA0POC/l9Sd0+PMWZfZN1KDkxEgPBwAuSHn9CY8ZrE5nLEs7DtWrouFk1P2MI2R5ERcXp+zsbBUXF+tb3/qW+/ni4mJNnTo1giMDus60aXUA8Ofw5y1+2yuPRV/IYrnQh8LCQv32t7/V7373O3388ce6//77VVVVpXvuuSfSQwO6xLRpdQDwZ+DlCX7bM/tGV8CSmMny6c4779SxY8f06KOP6ujRo8rKytJf/vIXDRw4MNJDAwIycVodAPwZnNJb44em6O1PGnTmK9WnYmw2jR2SHJU/05jJ8mPevHmqrKxUa2ur9u/fr/Hjx0d6SDBUuN9u3JVpdQAwTdFdozR2SLLHc2OHJKvorlERGpF/zGQBIRSpfVEmTqsDQCCOhFg9M+d6VTScVOWxk1H/hh5msoAQitS+qM5p9ZizSo7E2GwaPzQlqn8oAUAgg5ITdcuwK6L+ZxkhC8aL1sq/nfuizpz1yVVf3RcVSqZNqwNAd8NyIYwV7SUKIv12Y9Om1QGgu2EmC8aK9hIF0bIvypRpdQDobpjJglE6q5fH2BT1JQpMfLsxACB4CFkwgrelQX+ipfJv0V2j9MPnSj3Gzb4oALg0ELJgBG9Lg/5ES4kC9kUBwKWLkIWo56t6uTfRuhQ3KJlwBQCXGkIWolbn/qvaplNdPoalOABAtCBkIeqc7/6r38+5Xqc7LJbiAABRhZCFqNPV/VedS4PjhqaEYVQAAJwfQhYuWOdyXjBnkM5n/xVLgwCAaEbIwnkLZaX1QFXSV0wfqVRHL5YGAQBRj4rvOG+hrLQeqEr6mMF9qV4OADACIQvnJdQfetxZJT3GZvN4PsZm0/ihKYQrAIAxCFlwK68/oe0H6/wGpa586PHFKrprlMYOSfZ4jv1XAADTsCcL57XHKhwfekyVdABAd8BMFs5rj1U4l/MGJSey/woAYCxC1iXE23LgheyxYjkPAIDAWC68BPhbDuzKHquzZ5JYzgMAIDBmsi4B/pYDL2aPFct5AAD4RsjqZs5eEgy0HGj7/3upKJkAAEBwsVzYTfhaErxjdH+/x1UeO6miu0bph8+VehzLHisAAC4OIaub8LUk+H/tp/0el9k3kT1WAACEACGrG/D1ocpnLEsllY0aPdCpd6u+8FgyjLHZNHZIskeYGpRMuAIAIFjYk9UNBHqHYMGNmZRcAAAgzJjJMtwXLW1au/0Tv31G9HPomWtcLAcCABBGhCzDLXjuPb17+AuvbWcvCbIcCABA+BCyDOZrL1an6wZcxpIgAAARwp4sgwXaizXv1iHnfMAzAAAID0KWwS6mWjsAAAgtQhYAAEAIELIM1pUPdwYAAJFByDJYoOXCvolxYRoJAAA4GyHLYINTemv80BSf7ateLQvjaAAAwFcRsgy3MP8qn207DtWrooElQwAAIoGQZbiqRvZlAQAQjQhZhntmV6Xfdso4AAAQGYQsg5XXn1BJZaPP9tEDnXyMDgAAEULIMligEg4FN2aGZyAAAOAchCyDBSrhMKKfI0wjAQAAZyNkAQAAhAAhy2BUfAcAIHoRsgzmjI/12847CwEAiBxClsHWFB/y2TZ+aArvLAQAIIIIWYYqrz+hHYfqfbYv8lMJHgAAhB4hy1CB9mMda2kL00gAAIA3hCxDBSrfwH4sAAAii5AFAAAQAoQsQ1G+AQCA6EbIMhTlGwAAiG6ELENRvgEAgOhGyDIQ5RsAAIh+hCwDfXi02W875RsAAIg8QpaBntlV6bed/VgAAEQeIcsw5fUnVFLZ6LN99EAn+7EAAIgC3SZkVVZWas6cORo0aJDi4+N15ZVX6pFHHlFbm+fSWVVVlaZMmaLExEQlJydrwYIF5/Q5cOCA8vLyFB8fr379+unRRx+VZVnhvByfApVuONF6Wk0t7WEaDQAA8KVnpAcQLH//+9/V0dGhJ598UkOGDNEHH3yguXPn6uTJk1q1apUk6cyZM5o0aZJSUlK0c+dOHTt2TLNnz5ZlWSoqKpIkNTc367bbbtMtt9yikpISlZWVqaCgQImJiVq4cGEkL1FS4ErvBz87rh8+V6pn5lwfphEBAABvbFa0TNGEwMqVK7Vu3TqVl5dLkl5++WVNnjxZ1dXVcrlckqRNmzapoKBAdXV1SkpK0rp167R48WJ99tlnstvtkqQVK1aoqKhIR44ckc1mO+f7tLa2qrW11f11c3OzMjIy1NTUpKSkpKBf13effkc7D9Wrw0+f7YtuZtkQAIDz0NzcLIfDEbTf391mudCbpqYmXX755e6vd+/eraysLHfAkqQJEyaotbVV+/fvd/fJy8tzB6zOPjU1NaqsrPT6fZYvXy6Hw+F+ZGRkhOaC/r+iu0ZphMv/Xz4V3wEAiKxuG7I+/fRTFRUV6Z577nE/V1tbq9TUVI9+TqdTcXFxqq2t9dmn8+vOPmdbvHixmpqa3I/q6upgXso5HAmxemxalt8+vMMQAIDIivqQtXTpUtlsNr+Pffv2eRxTU1Ojr3/967r99tv1ve99z6PN23KfZVkez5/dp3NF1duxkmS325WUlOTxCDUqvgMAEN2ifuP7/PnzNWPGDL99MjMz3X+uqanRLbfcotzcXP3mN7/x6JeWlqa9e/d6PNfY2Kj29nb3bFVaWto5M1Z1dXWSdM4MV6RQ8R0AgOgX9SErOTlZycnJXer7j3/8Q7fccouys7O1fv169ejhOVGXm5urZcuW6ejRo0pPT5ckbdu2TXa7XdnZ2e4+S5YsUVtbm+Li4tx9XC6XR5iLpEBlHKj4DgBA5EX9cmFX1dTU6Oabb1ZGRoZWrVql+vp61dbWesxK5efna8SIEZo1a5ZKS0v12muvadGiRZo7d657iW/mzJmy2+0qKCjQBx98oC1btujxxx9XYWGhz+XCcAtUxoH9WAAARF7Uz2R11bZt2/TJJ5/ok08+Uf/+/T3aOvdUxcTEaOvWrZo3b57Gjh2r+Ph4zZw5011HS5IcDoeKi4t17733KicnR06nU4WFhSosLAzr9QAAALN16zpZkRLsOhtn236wTnevL/HZvv7u0bpl2BVB/74AAHRn1MmCnPGxfttZLgQAIPIIWQaifAMAANGPkGUYyjcAAGAGQpZhKN8AAIAZCFmGoXwDAABmIGQBAACEACHLMIGWCyuPnQzTSAAAgD+ELMOwXAgAgBkIWQAAACFAyDIMy4UAAJiBkGUYlgsBADADIcswg1N6y5ng/WN1nAmxVHsHACBKELIMU15/Qo0t7V7bGlvaVdHAciEAANGAkGUY9mQBAGAGQpZh2JMFAIAZCFkAAAAhQMgyDMuFAACYgZBlGGe893cWdmK5EACA6EDIMsya4kM+28YPTaGEAwAAUYKQZZDy+hPacajeZ/ui/KvCOBoAAOAPIcsggfZjHWtpC9NIAABAIIQsg1C+AQAAcxCyAAAAQoCQZRDKNwAAYA5ClkFYLgQAwByELAAAgBAgZBmE5UIAAMxByDII1d4BADAHIcsgVHsHAMAchCxDUO0dAACzELIMQbV3AADMQsgyBOUbAAAwCyELAAAgBAhZhqB8AwAAZiFkGYLlQgAAzELIMsTglN5yJnivk+VMiKV8AwAAUYaQZYjy+hNqbGn32tbY0q6KBpYLAQCIJoQsQ7AnCwAAsxCyDMGeLAAAzELIAgAACAFCliFYLgQAwCyELEOwXAgAgFkIWQAAACFAyDIEy4UAAJiFkGUIlgsBADALIcsQVHwHAMAshCxDUPEdAACzELIMwZ4sAADMQsgyBHuyAAAwCyELAAAgBAhZhmC5EAAAsxCyDMFyIQAAZiFkAQAAhAAhyxAsFwIAYBZCliGc8d4LkXZiuRAAgOhCyDLEmuJDPtvGD02h4jsAAFGGkGWA8voT2nGo3mf7ovyrwjgaAADQFd0yZLW2turaa6+VzWbTe++959FWVVWlKVOmKDExUcnJyVqwYIHa2to8+hw4cEB5eXmKj49Xv3799Oijj8qyrDBegadA+7GOtbT5bQcAAOHXM9IDCIUf//jHcrlcev/99z2eP3PmjCZNmqSUlBTt3LlTx44d0+zZs2VZloqKiiRJzc3Nuu2223TLLbeopKREZWVlKigoUGJiohYuXBiJy6F8AwAABup2Ievll1/Wtm3btHnzZr388ssebdu2bdNHH32k6upquVwuSdLq1atVUFCgZcuWKSkpSc8++6xOnTqlDRs2yG63KysrS2VlZVqzZo0KCwtls9nO+Z6tra1qbW11f93c3BzaiwQAAFGvWy0XfvbZZ5o7d65+//vfKyHh3Nmf3bt3Kysryx2wJGnChAlqbW3V/v373X3y8vJkt9s9+tTU1KiystLr912+fLkcDof7kZGREdTronwDAADm6TYhy7IsFRQU6J577lFOTo7XPrW1tUpNTfV4zul0Ki4uTrW1tT77dH7d2edsixcvVlNTk/tRXV19sZfjgeVCAADME/Uha+nSpbLZbH4f+/btU1FRkZqbm7V48WK/5/O23GdZlsfzZ/fp3PTu7VhJstvtSkpK8ngE0+CU3nImeK+T5UyIpXwDAABRKOr3ZM2fP18zZszw2yczM1OPPfaY9uzZ47HMJ0k5OTn6zne+o40bNyotLU179+71aG9sbFR7e7t7tiotLe2cGau6ujpJOmeGK1zK60+osaXda1tjS7sqGk4StAAAiDJRH7KSk5OVnJwcsN8vf/lLPfbYY+6va2pqNGHCBD3//PMaM2aMJCk3N1fLli3T0aNHlZ6eLunLzfB2u13Z2dnuPkuWLFFbW5vi4uLcfVwulzIzM4N8dV3TlT1ZhCwAAKJL1C8XdtWAAQOUlZXlflx11ZcFOq+88kr1799fkpSfn68RI0Zo1qxZKi0t1WuvvaZFixZp7ty57iW+mTNnym63q6CgQB988IG2bNmixx9/3Oc7C8OBPVkAAJin24SsroiJidHWrVvVq1cvjR07VnfccYemTZumVatWufs4HA4VFxfryJEjysnJ0bx581RYWKjCwsIIjhwAAJgm6pcLL1RmZqbXKu0DBgzQn//8Z7/Hjhw5Ujt27AjV0M4by4UAAJjnkprJMhXLhQAAmIeQBQAAEAKELANQ8R0AAPMQsgwQaLmwb2JcmEYCAAC6ipBlgMEpvTV+aIrP9lWvloVxNAAAoCsIWYZYmH+Vz7Ydh+pV0cCSIQAA0YSQZYjPW9r8trMvCwCA6ELIMgRlHAAAMAshCwAAIAQIWYagjAMAAGYhZBmC5UIAAMxCyDLE4JTecibEem1zJsTy2YUAAEQZQpYhyutPqLGl3WtbY0s7JRwAAIgyhCxDsCcLAACzELIMwZ4sAADMQsgCAAAIAUKWIT482uy3neVCAACiCyHLEM/sqvTbznIhAADRhZBlgPL6EyqpbPTZPnqgkxIOAABEGUKWAQK9s7DgxszwDAQAAHQZIcsAgd5ZmNHXfzsAAAg/QpYBBqf01vihKT7bV71aFsbRAACArujZ1Y6FhYVdPumaNWsuaDDwbWH+VdpxqN5r245D9apoOMm+LAAAokiXQ1ZpaanH1/v379eZM2c0bNgwSVJZWZliYmKUnZ0d3BFCkvR5S5vf9spjhCwAAKJJl0PW9u3b3X9es2aN+vTpo40bN8rpdEqSGhsbdffdd2vcuHHBHyWo+A4AgGEuaE/W6tWrtXz5cnfAkiSn06nHHntMq1evDtrgAAAATHVBIau5uVmfffbZOc/X1dXp+PHjFz0onIsPiAYAwCwXFLK+9a1v6e6779YLL7ygI0eO6MiRI3rhhRc0Z84cTZ8+PdhjhFguBADANF3ek/VVTzzxhBYtWqR/+7d/U3t7+5cn6tlTc+bM0cqVK4M6QHxpcEpvORNi1djSfk6bMyGWTe8AAESZCwpZCQkJWrt2rVauXKlPP/1UlmVpyJAhSkzkF32olNef8BqwJKmxpZ0SDgAARJkLClmdEhMT9a//+q/BGgv8CLQn68OaJkIWAABRpMsha/r06dqwYYOSkpIC7rt68cUXL3pg8BRoT9bGXZWa/K+uMI0GAAAE0uWQ5XA4ZLPZ3H9GeA1O6a2cgU7tO9zotb2kspElQwAAokiXQ9b69evdf167dq06Ojrce7AqKyv10ksvafjw4ZowYULwRwlJ0h3Z/X2GLImq7wAARJMLKuEwdepU/f73v5ckffHFF7rhhhu0evVqTZs2TevWrQvqAPFPf3z3iN92yjgAABA9Lihkvfvuu+6Pz3nhhReUmpqqw4cP65lnntEvf/nLoA4QXyqvP6GSSt+zWKMHOpnFAgAgilxQyGppaVGfPn0kSdu2bdP06dPVo0cP3XDDDTp8+HBQB4gvBXp3YcGNmeEZCAAA6JILCllDhgzRSy+9pOrqar366qvKz8+X9OXH6iQlJQV1gPhSoHcXZvT13w4AAMLrgkLWT37yEy1atEiZmZkaM2aMcnNzJX05qzVq1KigDhBfGpzSW+OHpvhsX/VqWRhHAwAAArFZlmVdyIG1tbU6evSorrnmGvXo8WVWe+edd5SUlKSrr746qIM0TXNzsxwOh5qamoI6s/d+9Rea+uu3fbZvX3Qz+7IAALhAwf79fcEV39PS0pSWlubx3PXXX3/RA4JvVY3+92VRwgEAgOhxQcuFiIxndlX6baeEAwAA0YOQZQhKOAAAYBZCliEo4QAAgFkIWYYIVMJhRD8+TxIAgGhCyDLE4JTecibEem1zJsSyVAgAQJQhZBmivP6EGlvavbY1trSrouFkmEcEAAD8IWQZItCerMpjhCwAAKIJIcsQgfZkUb4BAIDoQsgCAAAIAUKWIVguBADALIQsQ7BcCACAWQhZAAAAIUDIMgTLhQAAmIWQZQiWCwEAMAshyxBUfAcAwCyELENQ8R0AALN0u5C1detWjRkzRvHx8UpOTtb06dM92quqqjRlyhQlJiYqOTlZCxYsUFtbm0efAwcOKC8vT/Hx8erXr58effRRWZYVzss4B3uyAAAwS89IDyCYNm/erLlz5+rxxx/XrbfeKsuydODAAXf7mTNnNGnSJKWkpGjnzp06duyYZs+eLcuyVFRUJElqbm7WbbfdpltuuUUlJSUqKytTQUGBEhMTtXDhwkhdGnuyAAAwTLcJWadPn9aPfvQjrVy5UnPmzHE/P2zYMPeft23bpo8++kjV1dVyuVySpNWrV6ugoEDLli1TUlKSnn32WZ06dUobNmyQ3W5XVlaWysrKtGbNGhUWFspms53zvVtbW9Xa2ur+urm5OYRXCgAATNBtlgvfffdd/eMf/1CPHj00atQopaena+LEifrwww/dfXbv3q2srCx3wJKkCRMmqLW1Vfv373f3ycvLk91u9+hTU1OjyspKr997+fLlcjgc7kdGRkbQr4/lQgAAzNJtQlZ5ebkkaenSpfrP//xP/fnPf5bT6VReXp4+//xzSVJtba1SU1M9jnM6nYqLi1Ntba3PPp1fd/Y52+LFi9XU1OR+VFdXB/XaJJYLAQAwTdSHrKVLl8pms/l97Nu3Tx0dHZKkhx56SN/+9reVnZ2t9evXy2az6Y9//KP7fN6W+yzL8nj+7D6dm969HStJdrtdSUlJHo9go4QDAABmifo9WfPnz9eMGTP89snMzNTx48clSSNGjHA/b7fbNXjwYFVVVUmS0tLStHfvXo9jGxsb1d7e7p6tSktLO2fGqq6uTpLOmeEKp66UcCBoAQAQPaI+ZCUnJys5OTlgv+zsbNntdh08eFA33XSTJKm9vV2VlZUaOHCgJCk3N1fLli3T0aNHlZ6eLunLzfB2u13Z2dnuPkuWLFFbW5vi4uLcfVwulzIzM0NwhV3z4VH/m+krjxGyAACIJlG/XNhVSUlJuueee/TII49o27ZtOnjwoH7wgx9Ikm6//XZJUn5+vkaMGKFZs2aptLRUr732mhYtWqS5c+e6l/hmzpwpu92ugoICffDBB9qyZYsef/xxn+8sDJdndlX6bWdPFgAA0SXqZ7LOx8qVK9WzZ0/NmjVL//d//6cxY8bo9ddfl9PplCTFxMRo69atmjdvnsaOHav4+HjNnDlTq1atcp/D4XCouLhY9957r3JycuR0OlVYWKjCwsJIXZbK60+opLLRZ/vogU5msQAAiDI2K9KlzLuh5uZmORwONTU1BWUT/PaDdbp7fYnP9l/fNUqTrnH5bAcAAIEF+/d3t1ku7M4ClW8Y0c8RppEAAICuImQBAACEACHLAFR7BwDAPIQsA1DtHQAA8xCyDEC1dwAAzEPIMkBXqr0DAIDoQsgyQKA9WR/WNIVpJAAAoKsIWQYItCdrY4Bq8AAAIPwIWQYYnNJbOQOdPttLKhtZMgQAIMoQsgxx942Zftsp4wAAQHQhZBliuMt/eX/KOAAAEF0IWYagjAMAAGYhZBmCMg4AAJiFkGUIPloHAACzELIMwUfrAABgFkIWAABACBCyDMFyIQAAZiFkGYLlQgAAzELIAgAACAFCliFYLgQAwCyELEOwXAgAgFkIWYag4jsAAGYhZBmCiu8AAJiFkGUI9mQBAGAWQpYh2JMFAIBZCFkAAAAhQMgyBMuFAACYhZBlCJYLAQAwCyELAAAgBAhZhmC5EAAAsxCyDMFyIQAAZiFkGYKK7wAAmIWQZQgqvgMAYBZCliHYkwUAgFkIWYZgTxYAAGYhZAEAAIQAIcsQLBcCAGAWQpYhAi0X9k2MC9NIAABAVxCyDDE4pbfGD03x2b7q1bIwjgYAAARCyDLIwvyrfLbtOFRPGQcAAKIIIcsgVY3sywIAwBSELIM8s6vSbztlHAAAiB6ELEOU159QSWWjz/bRA518tA4AAFGEkGWIQCUcCm7MDM9AAABAlxCyDBGohMOIfo4wjQQAAHQFIQsAACAECFmGoOI7AABmIWQZgg+IBgDALIQsQwxO6S1nQqzXNmdCLO8sBAAgyhCyDFFef0KNLe1e2xpb2qn2DgBAlCFkGYI9WQAAmIWQZQj2ZAEAYBZCFgAAQAgQsgzBciEAAGYhZBmC5UIAAMxCyDIEJRwAADBLtwpZZWVlmjp1qpKTk5WUlKSxY8dq+/btHn2qqqo0ZcoUJSYmKjk5WQsWLFBbW5tHnwMHDigvL0/x8fHq16+fHn30UVmWFc5LOQclHAAAMEu3ClmTJk3S6dOn9frrr2v//v269tprNXnyZNXW1kqSzpw5o0mTJunkyZPauXOnNm3apM2bN2vhwoXuczQ3N+u2226Ty+VSSUmJioqKtGrVKq1ZsyZSlyVJ+vBos9929mQBABBdbFakp2iCpKGhQSkpKdqxY4fGjRsnSTp+/LiSkpL017/+VV/72tf08ssva/LkyaqurpbL5ZIkbdq0SQUFBaqrq1NSUpLWrVunxYsX67PPPpPdbpckrVixQkVFRTpy5IhsNlvAsTQ3N8vhcKipqUlJSUlBub7bn9ilkspGn+3bF93MkiEAABch2L+/u81MVt++fTV8+HA988wzOnnypE6fPq0nn3xSqampys7OliTt3r1bWVlZ7oAlSRMmTFBra6v279/v7pOXl+cOWJ19ampqVFlZ6fV7t7a2qrm52eMRTOX1J/wGrNEDnQQsAACiTLcJWTabTcXFxSotLVWfPn3Uq1cv/fznP9crr7yiyy67TJJUW1ur1NRUj+OcTqfi4uLcS4re+nR+3dnnbMuXL5fD4XA/MjIygnptgco33JHTP6jfDwAAXLyoD1lLly6VzWbz+9i3b58sy9K8efN0xRVX6K233tI777yjqVOnavLkyTp69Kj7fN6W+yzL8nj+7D6dK6q+lgoXL16spqYm96O6ujoYl+4WqHzDf+8/EtTvBwAALl7PSA8gkPnz52vGjBl++2RmZur111/Xn//8ZzU2NrrXUdeuXavi4mJt3LhRDz74oNLS0rR3716PYxsbG9Xe3u6erUpLSztnxqqurk6Szpnh6mS32z2WF4NtcEpv5Qx0at9h70uGJZWNqmg4yZIhAABRJOpDVnJyspKTkwP2a2n5ckmtRw/PybkePXqoo6NDkpSbm6tly5bp6NGjSk9PlyRt27ZNdrvdvW8rNzdXS5YsUVtbm+Li4tx9XC6XMjMzg3VZ5+3uGzN9hizpy3cXErIAAIgeUb9c2FW5ublyOp2aPXu23n//fZWVlek//uM/VFFRoUmTJkmS8vPzNWLECM2aNUulpaV67bXXtGjRIs2dO9c9+zVz5kzZ7XYVFBTogw8+0JYtW/T444+rsLCwS+8sDJXhLv/vcqDiOwAA0aXbhKzk5GS98sorOnHihG699Vbl5ORo586d+p//+R9dc801kqSYmBht3bpVvXr10tixY3XHHXdo2rRpWrVqlfs8DodDxcXFOnLkiHJycjRv3jwVFhaqsLAwUpcmiYrvAACYptvUyYomoaiTVV5/QreuftNnO3WyAAC4ONTJukQFKuNAxXcAAKILIcsQgco4sCcLAIDoQsgCAAAIAUKWIVguBADALIQsQ7BcCACAWQhZhqCEAwAAZiFkGaK8/oQaW9q9tjW2tKuigeVCAACiCSHLEIH2ZH1Y0xSmkQAAgK4gZBki0J6sjbsqwzMQAADQJYQsQwxO6a2cgU6f7SWVjSwZAgAQRQhZBrn7xky/7SwZAgAQPQhZBhnu8v85SiwZAgAQPQhZBmHJEAAAcxCyDBNoyXBP+bHwDAQAAPhFyDJMoCXDxS8e0HeffkdNPmpqAQCA8CBkGWZwSm+NH5qiGJvNZ5+3P2nQD58rDeOoAADA2QhZBlqYf5WGp/fx2X7GsrTjUD37swAAiKCekR4Auu6LljYteO497ThU36X+lcdO8pmGAABECDNZBlnw3Ht6+5OGLvfP7EvAAgAgUpjJMkR5/Ykuz2DF2GwaOySZWSwAACKImSxDBPqA6K8aOyRZRXeNCuFoAABAIMxkGcIZH+u3/fdzrtfpDkuZfROZwQIAIAoQsgyxpviQz7bxQ1M0bmhKGEcDAAACYbnQAIH2Yy3KvyqMowEAAF1ByDJAoP1Yx1rawjQSAADQVYQsAwy8PMFvO6UaAACIPoQsAACAECBkGSDQcmHlMT4+BwCAaEPIMgDLhQAAmIeQBQAAEAKELAOwXAgAgHkIWQZguRAAAPMQsgwwOKW3nAneP1bHmRDLx+gAABCFCFkGKK8/ocaWdq9tjS3tqmhguRAAgGhDyDIAe7IAADAPIcsA7MkCAMA8hCwAAIAQIGQZgOVCAADMQ8gyAMuFAACYh5BlAEo4AABgHkKWASjhAACAeQhZBgi0J+vDmqYwjQQAAHQVIcsAgfZkbdxVGZ6BAACALiNkGWBwSm/lDHT6bC+pbGTJEACAKEPIMsTdN2b6baeMAwAA0YWQZYjhriS/7ZRxAAAguhCyAAAAQoCQZQiqvgMAYBZCliGo+g4AgFkIWYag6jsAAGYhZBmCqu8AAJiFkGUI9mQBAGAWQpYh2JMFAIBZCFkAAAAhQMgyBMuFAACYhZBlCJYLAQAwCyHLEJRwAADALMaErGXLlunGG29UQkKCLrvsMq99qqqqNGXKFCUmJio5OVkLFixQW1ubR58DBw4oLy9P8fHx6tevnx599FFZluXR580331R2drZ69eqlwYMH64knngjVZXUZJRwAADCLMSGrra1Nt99+u37wgx94bT9z5owmTZqkkydPaufOndq0aZM2b96shQsXuvs0Nzfrtttuk8vlUklJiYqKirRq1SqtWbPG3aeiokLf+MY3NG7cOJWWlmrJkiVasGCBNm/eHPJr9Ic9WQAAmKVnpAfQVT/96U8lSRs2bPDavm3bNn300Ueqrq6Wy+WSJK1evVoFBQVatmyZkpKS9Oyzz+rUqVPasGGD7Ha7srKyVFZWpjVr1qiwsFA2m01PPPGEBgwYoF/84heSpOHDh2vfvn1atWqVvv3tb4fjUr1iTxYAAGYxZiYrkN27dysrK8sdsCRpwoQJam1t1f79+9198vLyZLfbPfrU1NSosrLS3Sc/P9/j3BMmTNC+ffvU3u59ua61tVXNzc0eDwAAcGnrNiGrtrZWqampHs85nU7FxcWptrbWZ5/OrwP1OX36tBoaGrx+7+XLl8vhcLgfGRkZQbmmr2K5EAAAs0Q0ZC1dulQ2m83vY9++fV0+n81mO+c5y7I8nj+7T+em9/Pt81WLFy9WU1OT+1FdXd3lMXcVy4UAAJglonuy5s+frxkzZvjtk5mZ2aVzpaWlae/evR7PNTY2qr293T0zlZaW5p6x6lRXVydJAfv07NlTffv29fq97Xa7xxJkKAxO6a3xQ1P09icNOvOVd0PG2GwaOySZEg4AAESZiIas5ORkJScnB+Vcubm5WrZsmY4ePar09HRJX26Gt9vtys7OdvdZsmSJ2traFBcX5+7jcrncYS43N1f/+7//63Hubdu2KScnR7Gx3utUhUvRXaP0w+dKteNQvfu5sUOSVXTXqAiOCgAAeGOzzi4SFaWqqqr0+eef609/+pNWrlypt956S5I0ZMgQ9e7dW2fOnNG1116r1NRUrVy5Up9//rkKCgo0bdo0FRUVSZKampo0bNgw3XrrrVqyZIkOHTqkgoIC/eQnP3GXeqioqFBWVpa+//3va+7cudq9e7fuuecePffcc11+d2Fzc7McDoeampqUlJQU9NeiouGkKo+dVGbfRGawAAAIkqD//rYMMXv2bEvSOY/t27e7+xw+fNiaNGmSFR8fb11++eXW/PnzrVOnTnmc529/+5s1btw4y263W2lpadbSpUutjo4Ojz5vvPGGNWrUKCsuLs7KzMy01q1bd15jbWpqsiRZTU1NF3y9AAAgvIL9+9uYmSyThHomCwAABF+wf393mxIOAAAA0YSQBQAAEAKELAAAgBAgZAEAAIQAIQsAACAECFkAAAAhQMgCAAAIAUIWAABACBCyAAAAQiCiHxDdXXUW0W9ubo7wSAAAQFd1/t4O1ofhELJC4Pjx45KkjIyMCI8EAACcr+PHj8vhcFz0efjswhDo6OhQTU2N+vTpI5vNdkHnaG5uVkZGhqqrq/n8Q/F6fBWvhSdeD0+8Hv/Ea+GJ1+OffL0WlmXp+PHjcrlc6tHj4ndUMZMVAj169FD//v2Dcq6kpKRL/h/DV/F6/BOvhSdeD0+8Hv/Ea+GJ1+OfvL0WwZjB6sTGdwAAgBAgZAEAAIQAIStK2e12PfLII7Lb7ZEeSlTg9fgnXgtPvB6eeD3+idfCE6/HP4XrtWDjOwAAQAgwkwUAABAChCwAAIAQIGQBAACEACELAAAgBAhZEbR27VoNGjRIvXr1UnZ2tt566y2//d98801lZ2erV69eGjx4sJ544okwjTS0li9frtGjR6tPnz664oorNG3aNB08eNDvMW+88YZsNts5j7///e9hGnVoLF269JxrSktL83tMd70vJCkzM9Pr3/O9997rtX93uy927NihKVOmyOVyyWaz6aWXXvJotyxLS5culcvlUnx8vG6++WZ9+OGHAc+7efNmjRgxQna7XSNGjNCWLVtCdAXB4++1aG9v1wMPPKCRI0cqMTFRLpdL3/3ud1VTU+P3nBs2bPB6v5w6dSrEV3PxAt0bBQUF51zXDTfcEPC8Jt4bUuDXw9vfs81m08qVK32eMxj3ByErQp5//nndd999euihh1RaWqpx48Zp4sSJqqqq8tq/oqJC3/jGNzRu3DiVlpZqyZIlWrBggTZv3hzmkQffm2++qXvvvVd79uxRcXGxTp8+rfz8fJ08eTLgsQcPHtTRo0fdj6FDh4ZhxKH1L//yLx7XdODAAZ99u/N9IUklJSUer0VxcbEk6fbbb/d7XHe5L06ePKlrrrlGv/rVr7y2/9d//ZfWrFmjX/3qVyopKVFaWppuu+029+enerN7927deeedmjVrlt5//33NmjVLd9xxh/bu3RuqywgKf69FS0uL3n33XT388MN699139eKLL6qsrEzf/OY3A543KSnJ4145evSoevXqFYpLCKpA94Ykff3rX/e4rr/85S9+z2nqvSEFfj3O/jv+3e9+J5vNpm9/+9t+z3vR94eFiLj++uute+65x+O5q6++2nrwwQe99v/xj39sXX311R7Pff/737duuOGGkI0xUurq6ixJ1ptvvumzz/bt2y1JVmNjY/gGFgaPPPKIdc0113S5/6V0X1iWZf3oRz+yrrzySqujo8Nre3e9LyzLsiRZW7ZscX/d0dFhpaWlWStWrHA/d+rUKcvhcFhPPPGEz/Pccccd1te//nWP5yZMmGDNmDEj6GMOlbNfC2/eeecdS5J1+PBhn33Wr19vORyO4A4uAry9HrNnz7amTp16XufpDveGZXXt/pg6dap16623+u0TjPuDmawIaGtr0/79+5Wfn+/xfH5+vnbt2uX1mN27d5/Tf8KECdq3b5/a29tDNtZIaGpqkiRdfvnlAfuOGjVK6enp+trXvqbt27eHemhhcejQIblcLg0aNEgzZsxQeXm5z76X0n3R1tamP/zhD/r3f//3gB+83h3vi7NVVFSotrbW4+/fbrcrLy/P588Ryfc94+8YEzU1Nclms+myyy7z2+/EiRMaOHCg+vfvr8mTJ6u0tDQ8AwyDN954Q1dccYWuuuoqzZ07V3V1dX77Xyr3xmeffaatW7dqzpw5Afte7P1ByIqAhoYGnTlzRqmpqR7Pp6amqra21usxtbW1XvufPn1aDQ0NIRtruFmWpcLCQt10003Kysry2S89PV2/+c1vtHnzZr344osaNmyYvva1r2nHjh1hHG3wjRkzRs8884xeffVVPfXUU6qtrdWNN96oY8eOee1/qdwXkvTSSy/piy++UEFBgc8+3fW+8KbzZ8X5/BzpPO58jzHNqVOn9OCDD2rmzJl+Pwj56quv1oYNG/SnP/1Jzz33nHr16qWxY8fq0KFDYRxtaEycOFHPPvusXn/9da1evVolJSW69dZb1dra6vOYS+HekKSNGzeqT58+mj59ut9+wbg/el7sYHHhzv7fuGVZfv+H7q2/t+dNNn/+fP3tb3/Tzp07/fYbNmyYhg0b5v46NzdX1dXVWrVqlcaPHx/qYYbMxIkT3X8eOXKkcnNzdeWVV2rjxo0qLCz0esylcF9I0tNPP62JEyfK5XL57NNd7wt/zvfnyIUeY4r29nbNmDFDHR0dWrt2rd++N9xwg8dm8LFjx+q6665TUVGRfvnLX4Z6qCF15513uv+clZWlnJwcDRw4UFu3bvUbLrrzvdHpd7/7nb7zne8E3FsVjPuDmawISE5OVkxMzDn/O6irqzvnfxGd0tLSvPbv2bOn+vbtG7KxhtMPf/hD/elPf9L27dvVv3//8z7+hhtu6Bb/A/2qxMREjRw50ud1XQr3hSQdPnxYf/3rX/W9733vvI/tjveFJPe7Ts/n50jnced7jCna29t1xx13qKKiQsXFxX5nsbzp0aOHRo8e3S3vl/T0dA0cONDvtXXne6PTW2+9pYMHD17Qz5ILuT8IWREQFxen7Oxs9zulOhUXF+vGG2/0ekxubu45/bdt26acnBzFxsaGbKzhYFmW5s+frxdffFGvv/66Bg0adEHnKS0tVXp6epBHF1mtra36+OOPfV5Xd74vvmr9+vW64oorNGnSpPM+tjveF5I0aNAgpaWlefz9t7W16c033/T5c0Tyfc/4O8YEnQHr0KFD+utf/3pB/8mwLEvvvfdet7xfjh07purqar/X1l3vja96+umnlZ2drWuuuea8j72g++Oits3jgm3atMmKjY21nn76aeujjz6y7rvvPisxMdGqrKy0LMuyHnzwQWvWrFnu/uXl5VZCQoJ1//33Wx999JH19NNPW7GxsdYLL7wQqUsImh/84AeWw+Gw3njjDevo0aPuR0tLi7vP2a/Hz3/+c2vLli1WWVmZ9cEHH1gPPvigJcnavHlzJC4haBYuXGi98cYbVnl5ubVnzx5r8uTJVp8+fS7J+6LTmTNnrAEDBlgPPPDAOW3d/b44fvy4VVpaapWWllqSrDVr1lilpaXud8ytWLHCcjgc1osvvmgdOHDAuuuuu6z09HSrubnZfY5Zs2Z5vGv57bfftmJiYqwVK1ZYH3/8sbVixQqrZ8+e1p49e8J+fefD32vR3t5uffOb37T69+9vvffeex4/R1pbW93nOPu1WLp0qfXKK69Yn376qVVaWmrdfffdVs+ePa29e/dG4hLPi7/X4/jx49bChQutXbt2WRUVFdb27dut3Nxcq1+/ft3y3rCswP9WLMuympqarISEBGvdunVezxGK+4OQFUG//vWvrYEDB1pxcXHWdddd51GyYPbs2VZeXp5H/zfeeMMaNWqUFRcXZ2VmZvq8UUwjyetj/fr17j5nvx4/+9nPrCuvvNLq1auX5XQ6rZtuusnaunVr+AcfZHfeeaeVnp5uxcbGWi6Xy5o+fbr14Ycfutsvpfui06uvvmpJsg4ePHhOW3e/LzpLUpz9mD17tmVZX5ZxeOSRR6y0tDTLbrdb48ePtw4cOOBxjry8PHf/Tn/84x+tYcOGWbGxsdbVV19tRAj191pUVFT4/Dmyfft29znOfi3uu+8+a8CAAVZcXJyVkpJi5efnW7t27Qr/xV0Af69HS0uLlZ+fb6WkpFixsbHWgAEDrNmzZ1tVVVUe5+gu94ZlBf63YlmW9eSTT1rx8fHWF1984fUcobg/bJb1/3fJAgAAIGjYkwUAABAChCwAAIAQIGQBAACEACELAAAgBAhZAAAAIUDIAgAACAFCFgAAQAgQsgAAAEKAkAUAABAChCwAAIAQIGQBwHlob2+P9BAAGIKQBeCS9sorr+imm27SZZddpr59+2ry5Mn69NNPJUmVlZWy2Wz67//+b918883q1auX/vCHP0iS1q9fr+HDh6tXr166+uqrtXbtWo/zPvDAA7rqqquUkJCgwYMH6+GHHyagAZeYnpEeAABE0smTJ1VYWKiRI0fq5MmT+slPfqJvfetbeu+999x9HnjgAa1evVrr16+X3W7XU089pUceeUS/+tWvNGrUKJWWlmru3LlKTEzU7NmzJUl9+vTRhg0b5HK5dODAAc2dO1d9+vTRj3/84whdKYBws1mWZUV6EAAQLerr63XFFVfowIED6t27twYNGqRf/OIX+tGPfuTuM2DAAP3sZz/TXXfd5X7uscce01/+8hft2rXL63lXrlyp559/Xvv27Qv5NQCIDsxkAbikffrpp3r44Ye1Z88eNTQ0qKOjQ5JUVVWlESNGSJJycnLc/evr61VdXa05c+Zo7ty57udPnz4th8Ph/vqFF17QL37xC33yySc6ceKETp8+raSkpDBdFYBoQMgCcEmbMmWKMjIy9NRTT8nlcqmjo0NZWVlqa2tz90lMTHT/uTOEPfXUUxozZozHuWJiYiRJe/bs0YwZM/TTn/5UEyZMkMPh0KZNm7R69eowXBGAaEHIAnDJOnbsmD7++GM9+eSTGjdunCRp586dfo9JTU1Vv379VF5eru985zte+7z99tsaOHCgHnroIfdzhw8fDt7AARiBkAXgkuV0OtW3b1/95je/UXp6uqqqqvTggw8GPG7p0qVasGCBkpKSNHHiRLW2tmrfvn1qbGxUYWGhhgwZoqqqKm3atEmjR4/W1q1btWXLljBcEYBoQgkHAJesHj16aNOmTdq/f7+ysrJ0//33a+XKlQGP+973vqff/va32rBhg0aOHKm8vDxt2LBBgwYNkiRNnTpV999/v+bPn69rr71Wu3bt0sMPPxzqywEQZXh3IQAAQAgwkwUAABAChCwAAIAQIGQBAACEACELAAAgBAhZAAAAIUDIAgAACAFCFgAAQAgQsgAAAEKAkAUAABAChCwAAIAQIGQBAACEwP8DzJ0k2MiySlsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "polygons.plot.scatter(x='area', y='sid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "f2cf3866-659a-4594-a2a7-fb7b78a5e123",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "slope_count=0\n",
    "crest = 3\n",
    "while i < len(polygons):\n",
    "    # row = polygons.iloc[i]\n",
    "    if polygons.iloc[i]['mean_slope'] > 1:\n",
    "        slope_count += 1\n",
    "    else:\n",
    "        slope_count = 0\n",
    "    # print(i,slope_count)\n",
    "    if slope_count >= crest:\n",
    "        non_eddy = polygons.iloc[0:i-(crest-1)]#.dissolve()\n",
    "        \n",
    "        break\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "6658e049-8bf8-4d24-9303-8ad866a4a986",
   "metadata": {},
   "outputs": [],
   "source": [
    "non_eddy.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/non_eddy.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "5ee6d40d-240c-45e5-bf9f-7f2a696c92a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "non_eddy = non_eddy.dissolve().explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "ccc6207b-3bb1-4e18-b9eb-409d9073d897",
   "metadata": {},
   "outputs": [],
   "source": [
    "non_eddy.geometry = non_eddy.geometry.apply(lambda p: close_holes(p))\n",
    "non_eddy = non_eddy.dissolve().explode()\n",
    "non_eddy.to_file(outdir+r\"/non_eddy_nohole1.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19c21886-1ba6-448a-9790-6db775909372",
   "metadata": {},
   "source": [
    "### Land mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "2cb71938-f326-4049-b9df-b9936d6d1baa",
   "metadata": {},
   "outputs": [],
   "source": [
    "land_min=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "71f2e347-451d-44df-bcc5-422370a00278",
   "metadata": {},
   "outputs": [],
   "source": [
    "land_mask(path=rast_path,hmin=land_min,outdir=outdir+r\"/land.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fa60d1-db5e-40f2-b6f3-647508562300",
   "metadata": {},
   "source": [
    "## Find transition zones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "86a3f971-fb54-4cf7-bbc6-bec05c245606",
   "metadata": {},
   "outputs": [],
   "source": [
    "get_survey_shapes(path=slope_path,hmax=50,hmin=0,outdir=outdir+r\"/allslopes.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "459bed94-e0c0-4d5f-bed1-3d515878db54",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/allslopes.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "8a11a412-1a89-44cb-878b-2d5afeead07b",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = slope.loc[slope['sid'] >5].dissolve().explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "2265c864-90c0-48d6-86b1-2dd36a838981",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope.to_file(outdir+r\"/slope_morethan5.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "01e9e346-dfac-4fd6-a197-9cbaecc3d9e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope.geometry = slope.geometry.apply(lambda p: close_holes(p))\n",
    "slope.to_file(outdir+r\"/slope_morethan5_nohole.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d932af2-dc38-4ebe-a1ae-f461a8fcf766",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterstats import zonal_stats\n",
    "\n",
    "# Read the raster\n",
    "with rasterio.open(rast_path) as raster:\n",
    "    affine = raster.transform\n",
    "    array = raster.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceff9c8c-2aa7-46e5-b1ac-c45a607eb5b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = slope.explode()\n",
    "stats = zonal_stats(polygons, array, affine=affine, stats=[\"min\",\"max\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04234f59-bc41-458d-a0e3-f2b69e979036",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate percentage\n",
    "list_min,list_max=[],[]\n",
    "for stat in stats:\n",
    "    list_min.append(stat['min'])\n",
    "    list_max.append(stat['max'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29b07df1-d5a0-4fc0-9263-7307beaf6b0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons['min_elev']=list_min\n",
    "polygons['max_elev']=list_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9918700-985a-4504-a539-0908bcd9ea58",
   "metadata": {},
   "outputs": [],
   "source": [
    "transition = polygons.loc[(polygons['max_elev'] >-100 ) & (polygons['min_elev'] < -500)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99505653-99a7-4df2-b4e0-f01ff46d22a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "transition.to_file(outdir+r\"/transition5.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "8679419b-fb7b-4f0f-a1fd-ae93dc7f2a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "# after this I could fill all holes, and then subtract by the non-eddy regime area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5bebc82-d769-412d-880a-a185ab59e871",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "300aec73-6d83-4bf6-a693-1f63f6027ee3",
   "metadata": {},
   "source": [
    "## Find eddy zones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "a93577af-a40c-4e05-8507-62aaf2a81cce",
   "metadata": {},
   "outputs": [],
   "source": [
    "eddy_max=-450"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "d7351974-1bfe-456d-95e1-03ab066bd7da",
   "metadata": {},
   "outputs": [],
   "source": [
    "eddy(path=rast_path,hmax=eddy_max,outdir=outdir+r\"/eddy450.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "988f25bb-1825-4213-b420-c554794ba5f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "eddy = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/eddy450.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bca34405-4ecc-44a1-a4a5-65857eb25124",
   "metadata": {},
   "source": [
    "## Find transition zones2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "489ebe6a-0d4a-4a46-ba31-6cdfee385199",
   "metadata": {},
   "outputs": [],
   "source": [
    "not_trans= pd.concat([eddy,non_eddy])\n",
    "not_trans = not_trans.dissolve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "4405aca5-1938-478d-9ab5-0dead5ed7c03",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, ymin, xmax, ymax = not_trans['geometry'].bounds.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "8790aed6-2a4e-4bf4-a60f-26643c91e175",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_box = box(xmin, ymin, xmax, ymax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "e5251165-1c73-444e-b6d1-d5065d7c44ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "transition = gpd.GeoDataFrame(geometry=[shape_box.symmetric_difference(not_trans.union_all())],\n",
    "                               #crs=crs\n",
    "                              ).dissolve().explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "f5a1c9b2-ce20-4567-a6c3-e06b97afd4e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\.conda\\envs\\ocsmesh_161\\Lib\\site-packages\\pyogrio\\raw.py:698: UserWarning: 'crs' was not provided.  The output dataset will not have projection information defined and may not be usable in other systems.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "transition.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/transition.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2abeeefa-d826-4220-a5bf-7661ca7dcd8a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06baa6eb-b209-4d0e-b132-a239aef2a1ab",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "29680a5c-8c20-4093-9d8d-37dd37d405d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "not_eddy = pd.concat([transition,non_eddy])\n",
    "not_eddy = not_eddy.dissolve()\n",
    "# xmin, ymin, xmax, ymax = not_eddy['geometry'].bounds.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "47f38fb1-a40a-4f93-87e5-ad21c37c9ebb",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, ymin, xmax, ymax = transition.dissolve()['geometry'].bounds.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "92bc4747-d156-4a8a-9b15-183b10fc04f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_box = box(xmin, ymin, xmax, ymax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "d2d52ea0-67b0-437d-a416-d5b6702f78f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "other_zones = gpd.GeoDataFrame(geometry=[shape_box.symmetric_difference(not_eddy.union_all())],\n",
    "                               #crs=crs\n",
    "                              ).dissolve().explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "5fc53d68-5759-483c-93f6-4cb9bad8fdb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\.conda\\envs\\ocsmesh_161\\Lib\\site-packages\\pyogrio\\raw.py:698: UserWarning: 'crs' was not provided.  The output dataset will not have projection information defined and may not be usable in other systems.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "other_zones.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/other_zones.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "17bc6ce7-92f5-4325-8df1-0bdc1c51069d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterstats import zonal_stats\n",
    "\n",
    "# Read the raster\n",
    "with rasterio.open(rast_path) as raster:\n",
    "    affine = raster.transform\n",
    "    array = raster.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "31eccd9d-a73f-41ed-9ba2-ab17b1362e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = other_zones.explode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "9258add3-0462-4e21-8122-3ae931149fa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = zonal_stats(polygons, array, affine=affine, stats=\"median\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "3d40aada-b1b2-4fcd-8f81-ce049301c049",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate percentage\n",
    "median=[]\n",
    "for stat in stats:\n",
    "    m = stat[\"median\"]\n",
    "    median.append(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "b6fbb0cc-0acc-4ef2-b4ad-35a11d2b84c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons['median']=median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "a1fff095-0bc4-4ba7-851b-20fc300ca2ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "eddy = polygons.loc[(polygons[\"median\"] <-750) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "ad8e8b6c-6cd1-4927-a801-8c4ec36ebf91",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\.conda\\envs\\ocsmesh_161\\Lib\\site-packages\\pyogrio\\raw.py:698: UserWarning: 'crs' was not provided.  The output dataset will not have projection information defined and may not be usable in other systems.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "eddy.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/eddy.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20dfb4fd-d3b3-4564-906e-fd9254bbd053",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c013e43-3a63-44d9-8a56-605d80266fdf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49e2578d-2751-4d79-b84d-ca2681ccc384",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5679ffbc-6bac-48ff-b7a8-396bda62b3b0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "195e0f68-fe71-435b-88ca-88617bf1b229",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75390c25-df39-4d8b-bd7c-f20fe395089e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "be7be38f-1eda-4bc0-9ffb-65db884ddab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.loc[(elev['sid'] > -20) & (elev['sid'] < -10), 'range'] = '-20--10'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "884fb2bb-ec31-49c0-bd17-c152f02c08d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.loc[(elev['sid'] > -30) & (elev['sid'] < -20), 'range'] = '-30--20'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2a82e469-71a6-475a-8246-b69c84007c80",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.loc[(elev['sid'] > -40) & (elev['sid'] < -30), 'range'] = '-40--30'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "2854c4df-5372-4ab1-8bcd-d1083e34a61e",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.loc[(elev['sid'] > -50) & (elev['sid'] < -40), 'range'] = '-50--40'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "156d6dfe-abb8-441a-9e33-e5ef88fc44a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = elev.dissolve(by='range').dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "403fa011-497e-45ed-a3cc-a0fdee4a05d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev.to_file(outdir+r\"/greater_50_by10.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eb194f5-2b89-49b0-a43d-7ad290b8c498",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "b1d6c41f-3d16-4985-8ade-7f865aeff75b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-9.9990000e+03, -9.9990000e+03, -9.9990000e+03, ...,\n",
       "        -9.9990000e+03, -9.9990000e+03, -9.9990000e+03],\n",
       "       [-9.9990000e+03,  4.8916068e-02,  9.7832069e-02, ...,\n",
       "         6.2382107e+00,  6.0996137e+00, -9.9990000e+03],\n",
       "       [-9.9990000e+03,  4.8916068e-02,  9.0196721e-02, ...,\n",
       "         1.6019108e+00,  1.0691147e+00, -9.9990000e+03],\n",
       "       ...,\n",
       "       [-9.9990000e+03,  3.0937243e-02,  2.1875935e-02, ...,\n",
       "         2.4457891e-01,  3.0937243e-02, -9.9990000e+03],\n",
       "       [-9.9990000e+03,  4.8916068e-02,  6.5627784e-02, ...,\n",
       "         5.0739527e-01,  1.4674792e-01, -9.9990000e+03],\n",
       "       [-9.9990000e+03, -9.9990000e+03, -9.9990000e+03, ...,\n",
       "        -9.9990000e+03, -9.9990000e+03, -9.9990000e+03]], dtype=float32)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calculate_slope(rast_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "aab916cc-2f1c-4dd6-9120-c1da809dedcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_path =r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\slope.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "884accd6-f486-4ea8-9a6d-e3ff19478f92",
   "metadata": {},
   "outputs": [],
   "source": [
    "get_survey_shapes(path=slope_path,hmax=.5,hmin=0,outdir=outdir+r\"/slope05.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b5cdaa8f-c4af-4ba3-833d-131a50f6dc1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/slope05.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "c962449d-14fe-411b-ace7-2c926050e78d",
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = gpd.read_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/greater_50.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "fa3f2f0d-ab2c-4503-a2c6-f0fb15f42c1e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\AppData\\Local\\Temp\\1\\ipykernel_10168\\1626583161.py:1: UserWarning: `keep_geom_type=True` in overlay resulted in 594 dropped geometries of different geometry types than df1 has. Set `keep_geom_type=False` to retain all geometries\n",
      "  result = gpd.overlay(elev, slope, how='union')\n"
     ]
    }
   ],
   "source": [
    "result = gpd.overlay(elev, slope, how='union')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "1cf8b235-e53c-4cbb-95f7-28d71247777b",
   "metadata": {},
   "outputs": [],
   "source": [
    "result.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/result.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "8128442f-ab21-48d9-92f7-b9875e7c301f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#combined = pd.concat([elev.explode(),slope.explode()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38a68673-1de9-4775-9ab4-de15491d86ae",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77cb40d3-2e81-407b-a341-6073c8205b7c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "d1694ab6-12f0-45eb-b4d9-d2e871622e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import itertools\n",
    "# geoms = combined['geometry'].tolist()\n",
    "# union_iter = gpd.GeoDataFrame(gpd.GeoSeries([poly[0].union(poly[1]) for poly in  itertools.combinations(geoms, 2) if poly[0].union(poly[1])]), columns=['geometry'])\n",
    "# union_iter.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\OCSMesh_Contributions\\ocean_mesh\\outputs/intersection_iter.shp\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "839218ac-fe89-48ba-a6ce-eefefc5827f5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ocsmesh_161",
   "language": "python",
   "name": "ocsmesh_161"
  },
  "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}