GeoTIFF_tiler_clean.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a3d87967",
   "metadata": {},
   "source": [
    "This script will create a tiled version of a raster dataset, it is expected that the input imagery is either one large raster or multible edge matched rasters of the same resolution. You are able to define the size of the tiles and the amount of overlap between tiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "74725eb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "from pathlib import Path\n",
    "\n",
    "import math\n",
    "from pathlib import Path\n",
    "from tqdm.auto import tqdm\n",
    "from multiprocess import Pool,cpu_count\n",
    "%matplotlib inline\n",
    "\n",
    "from osgeo import gdal,osr\n",
    "import geopandas as gpd\n",
    "from shapely import geometry\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import json\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import figure\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a79e9cac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geotiff folder good? True\n",
      "Output folder good? True\n"
     ]
    }
   ],
   "source": [
    "tile_size_px = [500,500] #x,y\n",
    "tile_oxerlap_px = 150\n",
    "\n",
    "geotiff_folder = '/data/Bike paths/inf image'\n",
    "output_folder = '/data/Bike paths/inf tiles'\n",
    "print('Geotiff folder good?', os.path.isdir(geotiff_folder))\n",
    "print('Output folder good?', os.path.isdir(output_folder))\n",
    "#if you want to use all of your rasters set the below value to None\n",
    "valid_rasters_csv = ''\n",
    "    \n",
    "input_file_ext = '.tif'\n",
    "output_compression = 'JPEG'  #use JPEG for images and LZW for raw data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9e64552c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do you want the output tiles to be nested into row folders?\n",
    "output_to_row_folders = True\n",
    "# do you want to append a string to every tile?\n",
    "raster_name_append = ''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "05380c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# dump tile setting a a json file into output location\n",
    "parent_dir = os.path.dirname(output_folder)\n",
    "settings_file_path = os.path.join(parent_dir,'tiling_settings.json')\n",
    "settings = {'tile_size_px':tile_size_px,\n",
    "           'tile_oxerlap_px':tile_oxerlap_px,\n",
    "           'raster_name_append':raster_name_append,\n",
    "           'tile_format':'.tif'}\n",
    "with open(settings_file_path, 'w') as fp:\n",
    "    json.dump(settings, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "overall-piece",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Path(output_folder).mkdir(parents=True, exist_ok=True)\n",
    "os.path.isdir(geotiff_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ba747de6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We found  1 .tif files\n"
     ]
    }
   ],
   "source": [
    "# search folder and all sub folders for 'input_file_ext' files\n",
    "geo_tiff_list = []\n",
    "for root, dirs, files in os.walk(geotiff_folder):\n",
    "    for file in files:\n",
    "        if file.endswith(input_file_ext):\n",
    "            geo_tiff_list.append(os.path.join(root, file))\n",
    "            \n",
    "print('We found ',len(geo_tiff_list),input_file_ext,'files')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "648aadea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# if you want to filter our input raster use this\n",
    "if valid_rasters_csv:\n",
    "    valid_raster_df = pd.read_csv(valid_rasters_csv)\n",
    "    print(valid_raster_df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6b72565e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert the valid raster names into a set so the lookup if much faster\n",
    "if valid_rasters_csv:\n",
    "    valid_raster_names = set(valid_raster_df.name.tolist())\n",
    "    valid_raster_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4d87f981",
   "metadata": {},
   "outputs": [],
   "source": [
    "if valid_rasters_csv:\n",
    "    short_list = []\n",
    "    for tif in tqdm(geo_tiff_list):\n",
    "    #     my csv file does not have the file extension so we are stripping out the '.tif' to find a match\n",
    "        file_name = os.path.basename(tif)#.replace('.tif','') \n",
    "    #     check if the current raster is in the list, if so add it to the clean list\n",
    "        if file_name in valid_raster_names:\n",
    "            short_list.append(tif)\n",
    "    # reset the full list to the clean list        \n",
    "    geo_tiff_list = short_list\n",
    "    print(len(geo_tiff_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1f9f0746",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_bounds(tif_path):\n",
    "#     open file\n",
    "    data = gdal.Open(tif_path)\n",
    "#     grab bounds\n",
    "    geoTransform = data.GetGeoTransform()\n",
    "    left = geoTransform[0]\n",
    "    top = geoTransform[3]\n",
    "    right = left + geoTransform[1] * data.RasterXSize\n",
    "    bottom = top + geoTransform[5] * data.RasterYSize\n",
    "#     build dict to file bounds\n",
    "    geo_tiff_bounds_dict = {'top':top,'left':left,'bottom':bottom,'right':right,'tif_path':tif_path}\n",
    "    return geo_tiff_bounds_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8448352a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'top': 6458672.2708,\n",
       " 'left': 394793.631,\n",
       " 'bottom': 6457734.1941,\n",
       " 'right': 395827.8306,\n",
       " 'tif_path': '/data/Bike paths/inf image/inf clip.tif'}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test the above function by returning the bounds of the first raster\n",
    "geo_tiff_bounds = get_bounds(geo_tiff_list[0])\n",
    "geo_tiff_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "421f4aa4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b252169a7c54e3d8c7aeb807f735d50",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use multiprocessing to extract raster bounds\n",
    "with Pool() as pool:\n",
    "    geo_tiff_bounds = list(tqdm(pool.imap(get_bounds, geo_tiff_list), total=len(geo_tiff_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8637a06a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make new array with only bounds to extract full raster extents\n",
    "pure_bounds = []\n",
    "for geo_tif_bounds in geo_tiff_bounds:\n",
    "    pure_bounds.append([geo_tif_bounds['top'],geo_tif_bounds['left'],geo_tif_bounds['bottom'],geo_tif_bounds['right']])\n",
    "# convert into numpy array\n",
    "pure_bounds_np = np.array(pure_bounds)\n",
    "# grab max extents\n",
    "bound_y_max = float(pure_bounds_np[:,0].max()) #top\n",
    "bound_x_min = float(pure_bounds_np[:,1].min()) #left\n",
    "bound_y_min = float(pure_bounds_np[:,2].min()) #bottom\n",
    "bound_x_max = float(pure_bounds_np[:,3].max()) #right"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "92b48823",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('top', 6458672.2708, 'left', 394793.631)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'top',bound_y_max,'left',bound_x_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "61e211e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1499926903553289 -0.14999627438440877\n"
     ]
    }
   ],
   "source": [
    "# open one image to get the pixel size, this is necessary to know how to cut the rasters\n",
    "test_raster = gdal.Open(geo_tiff_list[0])\n",
    "test_raster_gt =test_raster.GetGeoTransform()\n",
    "pixel_size_x = test_raster_gt[1]\n",
    "pixel_size_y = test_raster_gt[5]\n",
    "print(pixel_size_x,pixel_size_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "36e2df8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'EPSG:7850'"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proj = osr.SpatialReference(wkt=test_raster.GetProjection())\n",
    "crs = 'EPSG:'+proj.GetAttrValue('AUTHORITY',1)\n",
    "crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1d7c392c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "52.49744162436511 -52.498696034543066\n",
      "74.99634517766444 -74.99813719220438\n"
     ]
    }
   ],
   "source": [
    "# calculate the geographical distance in each direction each tile must be from the last tile\n",
    "x_move = pixel_size_x*(tile_size_px[0]-tile_oxerlap_px)\n",
    "y_move = pixel_size_y*(tile_size_px[1]-tile_oxerlap_px)\n",
    "print(x_move,y_move)\n",
    "\n",
    "# calculate the geographical size of each tile\n",
    "x_tile_size = pixel_size_x*tile_size_px[0]\n",
    "y_tile_size = pixel_size_y*tile_size_px[1]\n",
    "print(x_tile_size,y_tile_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ccd5502b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the number of cols so we can avoid using while loops\n",
    "number_of_cols = math.ceil(abs((bound_x_max-bound_x_min)/x_move))\n",
    "number_of_cols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d8d020a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the number of rows so we can avoid using while loops\n",
    "number_of_rows = math.ceil(abs((bound_y_max-bound_y_min)/y_move))\n",
    "number_of_rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "2f016525",
   "metadata": {},
   "outputs": [],
   "source": [
    "# will return a list of geotiffs which intersect \n",
    "def intersect_tile_with_geotiffs(tile_dict,geo_tiff_bounds):\n",
    "#     setup set to collect results in, a set is used to avoid duplicates\n",
    "    intersecting_geotiffs = set()\n",
    "#     loop over each geotiff\n",
    "    for geo_bounds in geo_tiff_bounds:\n",
    "#         check is tile top or bottom is inside geotiff\n",
    "        if (geo_bounds['top'] > tile_dict['top'] > geo_bounds['bottom'] or \n",
    "            geo_bounds['top'] > tile_dict['bottom'] > geo_bounds['bottom']):\n",
    "#         check if left or right are inside a geotiff\n",
    "            if geo_bounds['right'] > tile_dict['left'] > geo_bounds['left']:\n",
    "                intersecting_geotiffs.add(geo_bounds['tif_path'])\n",
    "            if geo_bounds['right'] > tile_dict['right'] > geo_bounds['left']:\n",
    "                intersecting_geotiffs.add(geo_bounds['tif_path'])\n",
    "    return intersecting_geotiffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "02ce53d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# will take tile bounds and only export them if they fall within a geotiff\n",
    "# this is called row by row by pool below\n",
    "def make_polygons(row):\n",
    "    tile_polygon_list = []\n",
    "    tile_top = bound_y_max + y_move*row\n",
    "    tile_bottom = tile_top + y_tile_size\n",
    "    tile_left = bound_x_min\n",
    "\n",
    "    for col in range(0,number_of_cols):\n",
    "        tile_left = bound_x_min + col*x_move\n",
    "        tile_right = tile_left + x_tile_size\n",
    "        tile_dict = {'top':tile_top,'left':tile_left,'bottom':tile_bottom,'right':tile_right}\n",
    "        tile_list = np.array([tile_top,tile_left,tile_bottom,tile_right])\n",
    "#         check if valid tile\n",
    "        intersect = intersect_tile_with_geotiffs(tile_dict,geo_tiff_bounds)\n",
    "        raster_name = raster_name_append+str(row)+'_'+str(col)+'.tif'\n",
    "        if len(intersect) > 0:\n",
    "            polygon = {'geometry':geometry.Polygon([[tile_left, tile_top], [tile_right, tile_top], [tile_right, tile_bottom], [tile_left, tile_bottom]]),\n",
    "                      'intersect':intersect, 'row':row, 'col':col, 'name':raster_name}\n",
    "            tile_polygon_list.append(polygon)\n",
    "    return tile_polygon_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f27b806b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ebdbdc20556f4dc5a135243c3babe742",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/18 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# multiprocess making polygons\n",
    "with Pool() as pool:\n",
    "    tile_polygon_list = list(tqdm(pool.imap(make_polygons, range(0,number_of_rows)), total=len(range(0,number_of_rows))))\n",
    "\n",
    "# this is returned as a list of list so it must be flattened\n",
    "tile_polygon_list = list(np.concatenate(tile_polygon_list).ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "46639a71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 13 ms, sys: 436 µs, total: 13.4 ms\n",
      "Wall time: 13.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#  convert into geodataframe\n",
    "polygon_tiles_gpd = gpd.GeoDataFrame(tile_polygon_list,geometry='geometry',crs=crs)\n",
    "del polygon_tiles_gpd['intersect']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "1bd40f42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>row</th>\n",
       "      <th>col</th>\n",
       "      <th>name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POLYGON ((394793.631 6458672.271, 394868.627 6...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0_0.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>POLYGON ((394846.128 6458672.271, 394921.125 6...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0_1.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>POLYGON ((394898.626 6458672.271, 394973.622 6...</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0_2.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>POLYGON ((394951.123 6458672.271, 395026.120 6...</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>0_3.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>POLYGON ((395003.621 6458672.271, 395078.617 6...</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0_4.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>355</th>\n",
       "      <td>POLYGON ((395581.093 6457779.793, 395656.089 6...</td>\n",
       "      <td>17</td>\n",
       "      <td>15</td>\n",
       "      <td>17_15.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>356</th>\n",
       "      <td>POLYGON ((395633.590 6457779.793, 395708.586 6...</td>\n",
       "      <td>17</td>\n",
       "      <td>16</td>\n",
       "      <td>17_16.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>357</th>\n",
       "      <td>POLYGON ((395686.088 6457779.793, 395761.084 6...</td>\n",
       "      <td>17</td>\n",
       "      <td>17</td>\n",
       "      <td>17_17.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>358</th>\n",
       "      <td>POLYGON ((395738.585 6457779.793, 395813.581 6...</td>\n",
       "      <td>17</td>\n",
       "      <td>18</td>\n",
       "      <td>17_18.tif</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>359</th>\n",
       "      <td>POLYGON ((395791.082 6457779.793, 395866.079 6...</td>\n",
       "      <td>17</td>\n",
       "      <td>19</td>\n",
       "      <td>17_19.tif</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>360 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              geometry  row  col       name\n",
       "0    POLYGON ((394793.631 6458672.271, 394868.627 6...    0    0    0_0.tif\n",
       "1    POLYGON ((394846.128 6458672.271, 394921.125 6...    0    1    0_1.tif\n",
       "2    POLYGON ((394898.626 6458672.271, 394973.622 6...    0    2    0_2.tif\n",
       "3    POLYGON ((394951.123 6458672.271, 395026.120 6...    0    3    0_3.tif\n",
       "4    POLYGON ((395003.621 6458672.271, 395078.617 6...    0    4    0_4.tif\n",
       "..                                                 ...  ...  ...        ...\n",
       "355  POLYGON ((395581.093 6457779.793, 395656.089 6...   17   15  17_15.tif\n",
       "356  POLYGON ((395633.590 6457779.793, 395708.586 6...   17   16  17_16.tif\n",
       "357  POLYGON ((395686.088 6457779.793, 395761.084 6...   17   17  17_17.tif\n",
       "358  POLYGON ((395738.585 6457779.793, 395813.581 6...   17   18  17_18.tif\n",
       "359  POLYGON ((395791.082 6457779.793, 395866.079 6...   17   19  17_19.tif\n",
       "\n",
       "[360 rows x 4 columns]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polygon_tiles_gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "43e015da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAJlCAYAAABaJA3OAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtTklEQVR4nO3dfbRlZ10n+O+PFAETEoIkIASwgh0yCxkTQhFF2pAQjSRLBTRqMr1aQMfIW5YwbQ/gTHdLa88oags9CLFGQF3NmwRD021IoG1txEXACgRIQjCxOk7KApKAIC9KSNVv/jj73jpcblWdqtzUvU/V57PWWXefZz/75Zy9bvLU9/72s6u7AwAAjOt+630CAADAvWNQDwAAgzOoBwCAwRnUAwDA4AzqAQBgcAb1AAAwuMNmUF9Vb6yqO6rqhgX7/0RV3VRVN1bVW+7r8wMAgPtKHS7z1FfV2Um+nOQPuvsJ++l7apI/TPL07v67qnpYd99xKM4TAADW2mGT1Hf3+5N8fr6tqr6jqq6uquuq6s+r6n+aVv1skt/u7r+btjWgBwBgWIfNoH4vtia5rLuflOQXkrxuan9cksdV1V9U1bVV9Yx1O0MAALiXNq33CdxXqupBSb43yTuqaqn5AdPPTUlOTXJOkkcl+fOqekJ3f+EQnyYAANxrh+2gPrO/Qnyhu89YZd2OJNd299eT/I+q+lRmg/y/PITnBwAAa+KwLb/p7r/PbMD+40lSM6dPq9+V5Nyp/cTMynG2r8d5AgDAvXXYDOqr6q1JPpjktKraUVU/k+SfJfmZqvpYkhuTPHPqfk2Sz1XVTUn+NMm/7O7Prcd5AwDAvXXYTGkJAABHqsMmqQcAgCOVQT0AAAzusJj95sQTT+zNmzev92kAALCGrrvuuru6+6TV1v3gucf25z6/61Cf0rLrPv61a7p7wzzr6LAY1G/evDnbtm1b79MAAGANVdXf7G3d5z6/Kx++5jGH8nS+wVGPuOXEdTv4Kg6LQT0AAEeWTrI7u9f7NDYMNfUAADA4ST0AAAPq7GpJ/RJJPQAADM6gHgAABqf8BgCA4cxulO31Po0NQ1IPAACDk9QDADAkU1ruIakHAIDBSeoBABhOp7Or1dQvkdQDAMDgJPUAAAzJ7Dd7SOoBAGBwBvUAADA45TcAAAynk+xSfrNMUg8AAIOT1AMAMCQ3yu4hqQcAgMFJ6gEAGE4nHj41R1IPAACDk9QDADCk3et9AhuIpB4AAAZnUA8AAINTfgMAwHA67eFTcwzqD9IP3O/HD+nxNj30oYf0eHnIgw/t8ZLs+tZjD+nxvv7gBxzS433twYf21+3u4+uwPt7XH3RID5evH3fo/8dxz4MObbVoHff1Q3q8b3nQ1w7p8U487iuH9HiPOOZLh/R4337M5w7p8ZJk8wPvOqTHO/XozxzS433H/b9wSI+3edNxh/R4B+p+3/ZX630K7INBPQAA4+lkl6B+mZp6AAAYnKQeAIDhdExpOU9SDwAAg5PUAwAwoMquHNpJGjYyST0AAAxOUg8AwHA6yW6z3yyT1AMAwOAM6gEAYHDKbwAAGJIbZfeQ1AMAwOAk9QAADKcjqZ8nqQcAgMFJ6gEAGNLultQvMagHAIA1VlVvTPJDSe7o7idMbW9PctrU5YQkX+juM1bZ9rYkX0qyK8k93b1lf8czqAcAYDgD1NT/XpLXJvmDpYbu/sml5ar6zSRf3Mf253b3XYsezKAeAADWWHe/v6o2r7auqirJTyR5+lodz42yAABwaH1fks929y17Wd9J3ltV11XVpYvsUFIPAMBwOpVd65tPn1hV2+beb+3urQtue0mSt+5j/VO7e2dVPSzJ+6rq5u5+/752aFAPAAAH7q5FbmBdqao2JfnRJE/aW5/u3jn9vKOqrkxyVhKDegAADj+DTmn5/Ulu7u4dq62sqmOT3K+7vzQtn5/k3+5vp2rqAQBgjVXVW5N8MMlpVbWjqn5mWnVxVpTeVNUjq+qq6e3Dk3ygqj6W5MNJ/ri7r97f8ST1AAAMZ6NPadndl+yl/bmrtO1McuG0vD3J6Qd6PEk9AAAMTlIPAMCAKrtaPr3ENwEAAIMzqAcAgMEpvwEAYDidZLd8eplvAgAABiepBwBgSBt5SstDTVIPAACDk9QDADCcblNazvNNAADA4CT1AAAMabea+mULJfVVdUJVXVFVN1fVJ6vqKXvp9+Sq2lVVF8213VZVn6iq66tq21z7GVV17VJ7VZ01t+67quqDVXXjtO0D782HBACAw9miSf1rklzd3RdV1dFJjlnZoaqOSvJrSa5ZZftzu/uuFW2vSvLK7n5PVV04vT+nqjYl+Y9J/nl3f6yqHprk6wueJwAAR4BOsksl+bL9fhNVdXySs5O8IUm6++7u/sIqXS9L8s4kdyx47E5y/LT84CQ7p+Xzk3y8uz82He9z3b1rwX0CAMARZ5F/3jw2yZ1J3lRVH62q362qY+c7VNXJSZ6d5PJVtu8k762q66rq0rn2lyT59aq6PclvJHnF1P64JF1V11TVR6rqfz+wjwQAAEeWRQb1m5KcmeT13f3EJF9J8vIVfV6d5GV7SdSf2t1nJrkgyYuq6uyp/QVJXtrdj07y0kx/CZiO90+T/LPp57Or6ryVO62qS6da/G133nnnAh8DAIDDx2xKy/V6bTSLnNGOJDu6+0PT+ysyG+TP25LkbVV1W5KLkryuqp6VJN29c/p5R5IrkyzdEPucJH80Lb9jrn1Hkv/e3Xd191eTXLXK8dLdW7t7S3dvOemkkxb4GAAAcHja76C+uz+T5PaqOm1qOi/JTSv6nNLdm7t7c2aD/hd297uq6tiqOi5JppKd85PcMG22M8nTpuWnJ7llWr4myXdV1THTTbNPW3k8AACObJ1kd+63bq+NZtHZby5L8uZp5pvtSZ5XVc9Pku5erY5+ycOTXFlVS8d6S3dfPa372SSvmQbu/5jk0ml/f1dV/z7JX2Z2va7q7j8+sI8FAABHjoUG9d19fWYlNvNWHcx393PnlrcnOX0v/T6Q5El7WfcfM5vWEgAAVrWrPXxqycb72wEAAHBAFi2/AQCADaNTHj41xzcBAACDM6gHAIDBKb8BAGBIuzfgQ6DWi28CAAAGJ6kHAGA4nbhRdo5vAgAABiepBwBgOJ3y8Kk5knoAABicpB4AgCHtlk8v800AAMDgDOoBAGBwym8AABhOd7LLw6eW+SYAAGBwknoAAAZU2R1TWi6R1AMAwOAk9QAADKejpn6ebwIAAAYnqQcAYEi75NPLfBMAADA4ST0AAMPpVHa32W+WSOoBAGBwBvUAADA45TcAAAzJjbJ7+CYAAGBwknoAAIbTSXZ7+NQy3wQAAAxOUg8AwIAqu2JKyyWSegAAGJykHgCA4aip/0a+CQAAGJxBPQAADE75DQAAQ3Kj7B6SegAAGJykHgCA4XSXG2Xn+CYAAGBwknoAAIa0S1K/zDcBAACDk9QDADCcTrLb7DfLJPUAADA4g3oAABic8hsAAAZUbpSd45sAAIDBSeoBABhOJ9ndbpRdIqkHAIDBGdQDADCkXbnfur32p6reWFV3VNUNc22/VFV/W1XXT68L97LtM6rqU1V1a1W9fJHvwqAeAADW3u8lecYq7b/V3WdMr6tWrqyqo5L8dpILkjw+ySVV9fj9HUxNPQAAw+nUhq6p7+73V9Xmg9j0rCS3dvf2JKmqtyV5ZpKb9rWRpB4AAA6dF1fVx6fynIessv7kJLfPvd8xte2TQT0AAEPanfut2yvJiVW1be516QKn/Pok35HkjCSfTvKbq/RZ7c8Pvb8dK78BAIADd1d3bzmQDbr7s0vLVfX/Jvkvq3TbkeTRc+8flWTn/vYtqQcAgEOgqh4x9/bZSW5YpdtfJjm1qk6pqqOTXJzk3fvbt6QeAIDhdCe7NvCNslX11iTnZFamsyPJv0lyTlWdkVk5zW1Jfm7q+8gkv9vdF3b3PVX14iTXJDkqyRu7+8b9Hc+gHgAA1lh3X7JK8xv20ndnkgvn3l+V5Jumu9wXg3oAAIa0kae0PNTU1AMAwOAk9QAADGf28Cn59BLfBAAADE5SDwDAkHat+pymI5OkHgAABmdQDwAAg1N+AwDAcDqmtJwnqQcAgMFJ6gEAGJApLef5JgAAYHCSegAAhrTblJbLJPUAADA4ST0AAMPpTnaZ/WaZpB4AAAYnqQcAYEhmv9nDNwEAAIMzqAcAgMEpvwEAYDidym43yi6T1AMAwOAk9QAADMnDp/aQ1AMAwOAk9QAADKcTNfVzJPUAADA4ST0AAEPy8Kk9fBMAADA4g3oAABic8hsAAMbTHj41T1IPAACDk9QDADCcjodPzZPUAwDA4CT1AAAMSU39HpJ6AAAYnKQeAIDhdCT18yT1AAAwOIN6AAAYnPIbAACGpPxmD0k9AAAMTlIPAMBwOiWpnyOpBwCAwUnqAQAY0u5I6pdI6gEAYHCSegAAxtNmv5m3UFJfVSdU1RVVdXNVfbKqnrKXfk+uql1VddFc221V9Ymqur6qts21n1FV1y61V9VZK/b1mKr6clX9wsF+OAAAOBIsmtS/JsnV3X1RVR2d5JiVHarqqCS/luSaVbY/t7vvWtH2qiSv7O73VNWF0/tz5tb/VpL3LHh+AAAcQTqS+nn7HdRX1fFJzk7y3CTp7ruT3L1K18uSvDPJkxc8dic5flp+cJKdc8d8VpLtSb6y4L4AAOCItUj5zWOT3JnkTVX10ar63ao6dr5DVZ2c5NlJLl9l+07y3qq6rqounWt/SZJfr6rbk/xGkldM+zo2ycuSvPJAPwwAAByJFhnUb0pyZpLXd/cTM0vPX76iz6uTvKy7d62y/VO7+8wkFyR5UVWdPbW/IMlLu/vRSV6a5A1T+yuT/FZ3f3lfJ1VVl061+NvuvPPOBT4GAACHk91d6/baaBYZ1O9IsqO7PzS9vyKzQf68LUneVlW3JbkoyeumEpp0987p5x1JrkyydEPsc5L80bT8jrn2707yqmlfL0nyi1X14pUn1d1bu3tLd2856aSTFvgYAABweNpvTX13f6aqbq+q07r7U0nOS3LTij6nLC1X1e8l+S/d/a6plOZ+3f2lafn8JP926rozydOS/FmSpye5ZdrX983t65eSfLm7X3vQnxAAgMNOZ2Mm5utl0dlvLkvy5mnmm+1JnldVz0+S7l6tjn7Jw5NcWVVLx3pLd189rfvZJK+pqk1J/jHJpavvAgAA2JeFBvXdfX1mJTbzVh3Md/dz55a3Jzl9L/0+kORJ+znuLy1yfgAAHHlaUr9soYdPAQAAG9ei5TcAALCh7I6kfomkHgAABmdQDwAAg1N+AwDAcLpjSss5knoAABicpB4AgCGZ0nIPST0AAAxOUg8AwIBKTf0cST0AAAxOUg8AwJDU1O8hqQcAgMEZ1AMAwOCU3wAAMJyOh0/Nk9QDAMAaq6o3VtUdVXXDXNuvV9XNVfXxqrqyqk7Yy7a3VdUnqur6qtq2yPEM6gEAGE8nvY6vBfxekmesaHtfkid093cl+askr9jH9ud29xndvWWRgxnUAwDAGuvu9yf5/Iq293b3PdPba5M8aq2Op6YeAIAh7c7QNfU/neTte1nXSd5bVZ3kd7p76/52ZlAPAAAH7sQV9e5bFxl8J0lV/R9J7kny5r10eWp376yqhyV5X1XdPCX/e2VQDwDAcDrr/vCpuxatd59XVc9J8kNJzutevTq/u3dOP++oqiuTnJVkn4N6NfUAAHAIVNUzkrwsyY9091f30ufYqjpuaTnJ+UluWK3vPEk9AAADqg09T31VvTXJOZmV6exI8m8ym+3mAZmV1CTJtd39/Kp6ZJLf7e4Lkzw8yZXT+k1J3tLdV+/veAb1AACwxrr7klWa37CXvjuTXDgtb09y+oEeT/kNAAAMTlIPAMCQFnwI1BFBUg8AAIOT1AMAMKR1ntJyQ5HUAwDA4CT1AAAMp1tSP09SDwAAg5PUAwAwpI388KlDTVIPAACDM6gHAIDBKb8BAGBIHj61h6QeAAAGJ6kHAGBIprTcQ1IPAACDk9QDADCcTknq50jqAQBgcJJ6AACGZPKbPST1AAAwOIN6AAAYnPIbAADG06a0nCepBwCAwUnqAQAYkztll0nqAQBgcJJ6AACGpKZ+D0k9AAAMTlIPAMCQWk39Mkk9AAAMTlIPAMBwOmrq50nqAQBgcAb1AAAwOOU3AACMp5Mov1kmqQcAgMFJ6gEAGJIpLfeQ1AMAwOAk9QAAjElSv0xSDwAAg5PUAwAwoPLwqTmSegAAGJxBPQAADE75DQAAY3Kj7DJJPQAADE5SDwDAeDpulJ0jqQcAgMFJ6gEAGJOa+mWSegAAGJykHgCAQampXyKpBwCAwRnUAwDA4JTfAAAwJjfKLpPUAwDA4CT1AACMSVK/TFIPAACDk9QDADCeTtKmtFwiqQcAgMFJ6gEAGFKrqV8mqQcAgMFJ6gEAGJOkfpmkHgAABmdQDwAAg1N+AwDAmExpuUxSDwAAg5PUAwAwpHKj7DJJPQAADE5SDwDAeDqmtJwjqQcAgMFJ6gEAGFCZ/WaOpB4AAAZnUA8AAGusqt5YVXdU1Q1zbd9aVe+rqlumnw/Zy7bPqKpPVdWtVfXyRY5nUA8AwJh6HV/793tJnrGi7eVJ/qS7T03yJ9P7b1BVRyX57SQXJHl8kkuq6vH7O5hBPQAArLHufn+Sz69ofmaS35+Wfz/Js1bZ9Kwkt3b39u6+O8nbpu32yY2yAACMaX2ntDyxqrbNvd/a3Vv3s83Du/vTSdLdn66qh63S5+Qkt8+935Hku/d3Mgb1AABw4O7q7i33wX5Xm9Jnv/98MagHAGBM4z186rNV9YgppX9EkjtW6bMjyaPn3j8qyc797VhNPQAAHBrvTvKcafk5Sf7TKn3+MsmpVXVKVR2d5OJpu30yqAcAYDyd2cOn1uu1H1X11iQfTHJaVe2oqp9J8qtJfqCqbknyA9P7VNUjq+qqJOnue5K8OMk1ST6Z5A+7+8b9HW+h8puqOiHJ7yZ5wvQV/nR3f3CVfk9Ocm2Sn+zuK6a225J8KcmuJPcs1R5V1RlJLk/ywCT3JHlhd3+4qpY+4NFJ7k7yL7v7vy1yngAAsBF09yV7WXXeKn13Jrlw7v1VSa46kOMtWlP/miRXd/dF058BjlnZYZpT89cy+1fFSud2910r2l6V5JXd/Z6qunB6f06Su5L8cHfvrKonTPs7ecHzBADgCFHj1dTfZ/Y7qK+q45OcneS5STLNl3n3Kl0vS/LOJE9e8Nid5Php+cGZbgDo7o/O9bkxyQOr6gHd/bUF9wsAAEeURZL6xya5M8mbqur0JNcl+fnu/spSh6o6Ocmzkzw93zyo7yTvrapO8jtz83e+JMk1VfUbmdX2f+8qx/6xJB81oAcAgL1b5EbZTUnOTPL67n5ikq/kmx9p++okL+vuXats/9TuPjOzR92+qKrOntpfkOSl3f3oJC9N8ob5jarqOzMr5/m51U6qqi6tqm1Vte3OO+9c4GMAAHBY6XV8bTCLDOp3JNnR3R+a3l+R2SB/3pYkb5tuir0oyeuq6lnJcuF/uvuOJFdm9ujbZDaNzx9Ny++Ya09VPWrq+1Pd/dernVR3b+3uLd295aSTTlrgYwAAwOFpv4P67v5Mktur6rSp6bwkN63oc0p3b+7uzZkN+l/Y3e+qqmOr6rgkqapjk5yf5IZps51JnjYtPz3JLVO/E5L8cZJXdPdf3IvPBgAAR4RFZ7+5LMmbp5lvtid5XlU9P0m6+/J9bPfwJFdW1dKx3tLdV0/rfjbJa6pqU5J/THLp1P7iJP8kyb+qqn81tZ0/Jf0AAMAKCw3qu/v6zEps5q06mO/u584tb09y+l76fSDJk1Zp/5Ukv7LIeQEAcOQypeUenigLAACDW7T8BgAANpau9T6DDUNSDwAAgzOoBwCAwSm/AQBgPBv0IVDrRVIPAACDk9QDADAmSf0yST0AAAxOUg8AwJA8fGoPST0AAAxOUg8AwJgk9csk9QAAMDiDegAAGJzyGwAAxqT8ZpmkHgAABiepBwBgONWmtJwnqQcAgMFJ6gEAGFPXep/BhiGpBwCAwUnqAQAYk5r6ZZJ6AAAYnKQeAIAhmf1mD0k9AAAMzqAeAAAGp/wGAIAxKb9ZJqkHAIDBSeoBABhPu1F2nqQeAAAGJ6kHAGBMkvplknoAABicpB4AgDFJ6pdJ6gEAYHAG9QAAMDjlNwAADMmUlntI6gEAYHAG9QAAMDiDegAAGJyaegAAxqSmfpmkHgAABiepBwBgPG32m3mSegAAGJxBPQAADE75DQAAY1J+s0xSDwAAg5PUAwAwJkn9Mkk9AAAMTlIPAMBwKqa0nCepBwCAwUnqAQAYk6R+maQeAAAGJ6kHAGA8raZ+nqQeAAAGZ1APAACDU34DAMCYlN8sk9QDAMDgJPUAAIxJUr9MUg8AAIOT1AMAMCRTWu4hqQcAgDVWVadV1fVzr7+vqpes6HNOVX1xrs+/PtjjSeoBABjTBk7qu/tTSc5Ikqo6KsnfJrlyla5/3t0/dG+PJ6kHAID71nlJ/rq7/+a+OoBBPQAAHLgTq2rb3OvSffS9OMlb97LuKVX1sap6T1V958GejPIbAADG01nv8pu7unvL/jpV1dFJfiTJK1ZZ/ZEk397dX66qC5O8K8mpB3MyknoAALjvXJDkI9392ZUruvvvu/vL0/JVSe5fVScezEEk9QAADGmQKS0vyV5Kb6rq25J8tru7qs7KLHD/3MEcxKAeAADuA1V1TJIfSPJzc23PT5LuvjzJRUleUFX3JPmHJBd390H9U8WgHgCAMW3wpL67v5rkoSvaLp9bfm2S167FsdTUAwDA4CT1AAAMaZCa+kNCUg8AAIMzqAcAgMEpvwEAYEzKb5ZJ6gEAYHCSegAAxtOR1M+R1AMAwOAk9QAADKemFzOSegAAGJykHgCAMampXyapBwCAwUnqAQAYUknql0nqAQBgcAb1AAAwOOU3AACMSfnNMkk9AAAMTlIPAMCYJPXLJPUAADA4ST0AAONpU1rOk9QDAMDgJPUAAIxJUr9MUg8AAIMzqAcAgMEpvwEAYEhulN1DUg8AAIOT1AMAMCZJ/TJJPQAADE5SDwDAkNTU7yGpBwCAwUnqAQAYT0dN/RxJPQAADG6hQX1VnVBVV1TVzVX1yap6yl76PbmqdlXVRXNtt1XVJ6rq+qraNtd+RlVdu9ReVWfNrXtFVd1aVZ+qqh+8Nx8QAAAOd4uW37wmydXdfVFVHZ3kmJUdquqoJL+W5JpVtj+3u+9a0faqJK/s7vdU1YXT+3Oq6vFJLk7ynUkemeS/VtXjunvXgucKAMCRQPnNsv0m9VV1fJKzk7whSbr77u7+wipdL0vyziR3LHjsTnL8tPzgJDun5WcmeVt3f627/0eSW5Octcr2AABAFkvqH5vkziRvqqrTk1yX5Oe7+ytLHarq5CTPTvL0JE9esX0neW9VdZLf6e6tU/tLklxTVb+R2T8uvndqPznJtXPb75jaAAAgSVIxpeW8RWrqNyU5M8nru/uJSb6S5OUr+rw6ycv2UiLz1O4+M8kFSV5UVWdP7S9I8tLufnSSl2b6S0Bm12ilb7pkVXXpVIu/7c4771zgYwAAwOFpkUH9jiQ7uvtD0/srMhvkz9uS5G1VdVuSi5K8rqqelSTdvXP6eUeSK7OnlOY5Sf5oWn7HXPuOJI+e2/ejsqc0Z1l3b+3uLd295aSTTlrgYwAAcFjpdXxtMPsd1Hf3Z5LcXlWnTU3nJblpRZ9Tuntzd2/ObND/wu5+V1UdW1XHJUlVHZvk/CQ3TJvtTPK0afnpSW6Zlt+d5OKqekBVnZLk1CQfPtgPCAAAh7tFZ7+5LMmbp5lvtid5XlU9P0m6+/J9bPfwJFdW1dKx3tLdV0/rfjbJa6pqU5J/THLptL8bq+oPM/uHwz1JXmTmGwAAVqregJH5OlloUN/d12dWYjNv1cF8dz93bnl7ktP30u8DSZ60l3X/Lsm/W+TcAADgSLdoUg8AABvHBq1tXy8LPVEWAADYuAzqAQBgcMpvAAAYkodP7SGpBwCAwUnqAQAYk6R+maQeAAAGJ6kHAGBIaur3kNQDAMDgJPUAAIxJUr9MUg8AAIMzqAcAgMEpvwEAYDztRtl5knoAABicpB4AgDFJ6pdJ6gEAYHCSegAAhlNRUz9PUg8AAIOT1AMAMKYW1S+R1AMAwOAk9QAADElN/R6SegAAGJxBPQAADE75DQAA4+l4+NQcST0AAAxOUg8AwJBq93qfwb5V1W1JvpRkV5J7unvLivWV5DVJLkzy1STP7e6PHMyxDOoBAOC+c25337WXdRckOXV6fXeS108/D5hBPQAAYxq/pv6ZSf6guzvJtVV1QlU9ors/faA7UlMPAAD3jU7y3qq6rqouXWX9yUlun3u/Y2o7YJJ6AACGtM4PnzqxqrbNvd/a3VtX9Hlqd++sqocleV9V3dzd759bX6vs96A+lUE9AAAcuLtW3vi6UnfvnH7eUVVXJjkryfygfkeSR8+9f1SSnQdzMspvAABgjVXVsVV13NJykvOT3LCi27uT/FTNfE+SLx5MPX0iqQcAYESdpDf0nbIPT3LlbNbKbErylu6+uqqenyTdfXmSqzKbzvLWzKa0fN7BHsygHgAA1lh3b09y+irtl88td5IXrcXxDOoBABjSOt8ou6GoqQcAgMFJ6gEAGJOkfpmkHgAABiepBwBgOBU19fMk9QAAMDiDegAAGJzyGwAAxtO90R8+dUhJ6gEAYHCSegAAhuRG2T0k9QAAMDhJPQAAY5LUL5PUAwDA4CT1AAAMSU39HpJ6AAAYnKQeAIDxdJLdovolknoAABicQT0AAAxO+Q0AAGNSfbNMUg8AAIOT1AMAMCRTWu4hqQcAgMFJ6gEAGFOL6pdI6gEAYHCSegAAhqSmfg9JPQAADM6gHgAABqf8BgCA8XQ8fGqOpB4AAAYnqQcAYDiVpExpuUxSDwAAg5PUAwAwpt3rfQIbh6QeAAAGJ6kHAGBIaur3kNQDAMDgDOoBAGBwym8AABiPh099A0k9AAAMTlIPAMCAOnGj7DJJPQAADE5SDwDAkEpQv0xSDwAAg5PUAwAwJjX1yyT1AAAwOEk9AADj6aR2r/dJbBySegAAGJxBPQAADE75DQAAY3Kj7DJJPQAADE5SDwDAmAT1yyT1AAAwOEk9AABDKjX1yyT1AAAwOEk9AABjktQvk9QDAMDgDOoBAGBwym8AABhPJ9m93iexcUjqAQBgcJJ6AACGU2lTWs6R1AMAwOAk9QAAjElSv0xSDwAAg5PUAwAwJkn9Mkk9AAAMbqFBfVWdUFVXVNXNVfXJqnrKXvo9uap2VdVFc223VdUnqur6qto21/72qe36qc/1U/v9q+r3p20+WVWvuJefEQAADmuLlt+8JsnV3X1RVR2d5JiVHarqqCS/luSaVbY/t7vvmm/o7p+c2/Y3k3xxevvjSR7Q3f9zVR2T5Kaqemt337bguQIAcLjz8KlvsN+kvqqOT3J2kjckSXff3d1fWKXrZUnemeSOAzmBqqokP5HkrVNTJzm2qjYl+ZYkdyf5+wPZJwAArKeqenRV/elUeXJjVf38Kn3OqaovzlWv/OuDPd4iSf1jk9yZ5E1VdXqS65L8fHd/Ze6ETk7y7CRPT/LkFdt3kvdWVSf5ne7eumL99yX5bHffMr2/Iskzk3w6s78IvLS7P39gHwsAgMPdBn/41D1J/kV3f6SqjktyXVW9r7tvWtHvz7v7h+7twRapqd+U5Mwkr+/uJyb5SpKXr+jz6iQv6+5dq2z/1O4+M8kFSV5UVWevWH9J9qT0SXJWkl1JHpnklCT/oqoeu3KnVXVpVW2rqm133nnnAh8DAAAOje7+dHd/ZFr+UpJPJjn5vjreIoP6HUl2dPeHpvdXZDbIn7clyduq6rYkFyV5XVU9K0m6e+f0844kV2Y2aE+STCU2P5rk7XP7+l8yq9//+rTNX0z7/wbdvbW7t3T3lpNOOmmBjwEAwGGle/1eB6CqNid5YpIPrbL6KVX1sap6T1V958F+Ffsd1Hf3Z5LcXlWnTU3nJblpRZ9Tuntzd2/ObND/wu5+V1UdO/25IVV1bJLzk9wwt+n3J7m5u3fMtf1/SZ5eM8cm+Z4kNx/cxwMAgPvEiUtVI9Pr0tU6VdWDMrvv9CXdvfI+0Y8k+fbuPj3J/5PkXQd7MovOfnNZkjdPM99sT/K8qnp+knT35fvY7uFJrpzdC5tNSd7S3VfPrb8431h6kyS/neRNmQ3+K8mbuvvjC54nAABHhANPzNfYXd39TdUk86rq/pkN6N/c3X+0cv38IL+7r6qq11XViStnjVzEQoP67r4+31wCs+pgvrufO7e8Pcnp+9jvc1dp+3Jm01oCAMCQphke35Dkk9397/fS59symzCmq+qszKpoPncwx1s0qQcAgI2js95J/f48Nck/T/KJpYesJvnFJI9JlqtdLkrygqq6J8k/JLm4++A+lEE9AACsse7+QGal5Pvq89okr12L4y0y+w0AALCBSeoBABjT7vU+gY1DUg8AAIOT1AMAMKTa2DfKHlKSegAAGJykHgCAMUnql0nqAQBgcJJ6AADG00l2S+qXSOoBAGBwBvUAADA45TcAAAyo3Sg7R1IPAACDk9QDADAmSf0yST0AAAxOUg8AwJgk9csk9QAAMDhJPQAA4/HwqW8gqQcAgMEZ1AMAwOCU3wAAMKBOevd6n8SGIakHAIDBSeoBABiTKS2XSeoBAGBwknoAAMZjSstvIKkHAIDBSeoBABiTmvplknoAABicpB4AgDFJ6pdJ6gEAYHAG9QAAMDjlNwAADKiV38yR1AMAwOAk9QAAjKeT7N693mexYUjqAQBgcJJ6AADGpKZ+maQeAAAGJ6kHAGBMkvplknoAABicQT0AAAxO+Q0AAAPqZLfymyWSegAAGJykHgCA8XTS7eFTSyT1AAAwOEk9AABjUlO/zKD+IL1v9zvW+xQAACCJQT0AAKPy8KllauoBAGBwBvUAADA45TcAAIynO9ltSsslknoAABicpB4AgDG5UXaZpB4AAAYnqQcAYEitpn6ZpB4AAAYnqQcAYECtpn6OpB4AAAYnqQcAYDydZLekfomkHgAABmdQDwAAg1N+AwDAmNqUlksk9QAAMDhJPQAAw+kk7UbZZZJ6AAAYnKQeAIDxdKupnyOpBwCAwUnqAQAYkpr6PST1AACwxqrqGVX1qaq6tapevsr6qqr/MK3/eFWdeW+OZ1APAABrqKqOSvLbSS5I8vgkl1TV41d0uyDJqdPr0iSvvzfHVH4DAMCYNu6NsmclubW7tydJVb0tyTOT3DTX55lJ/qC7O8m1VXVCVT2iuz99MAeU1AMAwNo6Ocntc+93TG0H2mdhh0VSf911191VVX+zzqdxYpK71vkcWIxrNQ7Xahyu1Vhcr3Ec6dfq2/e24kv5u2v+a19x4qE8mRUeWFXb5t5v7e6t03Kt0n/lXb2L9FnYYTGo7+6T1vscqmpbd29Z7/Ng/1yrcbhW43CtxuJ6jcO12rvufsZ6n8M+7Ejy6Ln3j0qy8yD6LEz5DQAArK2/THJqVZ1SVUcnuTjJu1f0eXeSn5pmwfmeJF882Hr65DBJ6gEAYKPo7nuq6sVJrklyVJI3dveNVfX8af3lSa5KcmGSW5N8Ncnz7s0xDerXztb9d2GDcK3G4VqNw7Uai+s1DtdqUN19VWYD9/m2y+eWO8mL1up4NdsfAAAwKjX1AAAwuCN6UF9VD6yqD1fVx6rqxqp65dR+elV9sKo+UVX/uaqOX7HdY6rqy1X1C3Ntl0z9P15VV1fViVP7A6rq7dMjgD9UVZvntnlOVd0yvZ5ziD72kA70WlXV5qr6h6q6fnpdPrevJ039b50ez1xTu2u1RtbqelXVMVX1x1V187SfX507huu1Btbyd2tun++uqhvm3rtWa2CN/zt4dFVtraq/mn6/fmxqd63WwBpfK+MLFtPdR+wrs/lBHzQt3z/Jh5J8T2Z3LD9tav/pJL+8Yrt3JnlHkl+Y3m9KckeSE6f3r0ryS9PyC5NcPi1fnOTt0/K3Jtk+/XzItPyQ9f5ONurrQK9Vks1JbtjLvj6c5CnTPt+T5ALXamNeryTHJDl3Wj46yZ+7XhvzWs3t70eTvGW+j2u18a5Vklcm+ZVp+X7Z8/8v12oDXasYX3gdwOuITup75svT2/tPr05yWpL3T+3vS/JjS9tU1bMy+wW5cW5XNb2OrapKcnz2zDP6zCS/Py1fkeS8qc8PJnlfd3++u/9uOs5Gnm91XR3MtVpNVT0iyfHd/cHu7iR/kORZ02rXao2s1fXq7q92959Oy3cn+Uhm8/gmrteaWKtrlSRV9aAk/1uSX1mxyrVaA2t5rTIbUP7f0353d/fSw41cqzWwhtfK+IKFHdGD+iSpqqOq6vrM/iX8vu7+UJIbkvzI1OXHMz0YoKqOTfKyzBKOZd399SQvSPKJzH7ZHp/kDdPq5UcAd/c9Sb6Y5KFZ40cDHwkO5FpNTqmqj1bVf6+q75vaTs7su14y/727Vmtoja7X/P5OSPLDSf5kanK91sgaXqtfTvKbmU3NNs+1WiNrca2m36Uk+eWq+khVvaOqHj61uVZrZC2ulfEFB+KIH9R3967uPiOz9O+sqnpCZgnGi6rquiTHJbl76v7KJL8196/vJElV3T+zX7onJnlkko8necXS6tUOu4929uIAr9Wnkzymu5+YWXL4lql2cV/fu2u1htboeiVJqmpTkrcm+Q/dvX2pebXD7qOdvViLa1VVZyT5J9195SqHcK3WyBr9Xm2atv+L7j4zyQeT/Ma0jWu1Rtbo98r4goUd8YP6Jd39hSR/luQZ3X1zd5/f3U/KbCDx11O3707yqqq6LclLkvxizR4scMa0j7+eSjr+MMn3TtssPwJ4Gpg8OMnns8aPBj6SLHKtuvtr3f25afm6qf1xmX3vj5rb3fz37lrdB+7l9VqyNckt3f3quTbXa43dy2v1lCRPmv77+IEkj6uqP5t27VqtsXt5rT6X2V9Tlv4B9o4kZ07LrtUau5fX6oypzfiC/TqiB/VVddLSnyGr6luSfH+Sm6vqYVPb/ZL8n0kuT5Lu/r7u3tzdm5O8Osn/1d2vTfK3SR5fVSdNu/6BJJ+clt+dZOnO84uS/LfpF/OaJOdX1UOq6iFJzp/aWMWBXqup/1HT8mOTnJpke88ev/ylqvqeqfbwp5L8p+kwrtUaWavrNb3/lcz+Z/WSFYdxvdbAGv5uvb67Hzn99/GfJvmr7j5nOoxrtQbW8Fp1kv+c5Jxp1+cluWladq3WwBr+N9D4gsX1Brhbd71eSb4ryUcz+3PWDUn+9dT+80n+anr9aqaHdK3Y9pcyzX4zvX9+Zr9oH8/sP5YPndofmFkKcmtms648dm6bn57ab03yvPX+Pjby60CvVWY3H92Y5GOZ3Vz5w3P72jLt46+TvHZuG9dqg12vzBKmnn63rp9e/6vrtfGu1Yp9bs43zn7jWm2wa5Xk2zO7YfPjmd2n8hjXasNeK+MLr4VenigLAACDO6LLbwAA4HBgUA8AAIMzqAcAgMEZ1AMAwOAM6gEAYHAG9QAAMDiDegAAGJxBPQAADO7/B4elUqU0gW+eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1332x756 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot up a bunch of polygons\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(18.5, 10.5)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.1)\n",
    "polygon_tiles_gpd.head(10000).plot(column='col', ax=ax, legend=True,cax=cax);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "b450090c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAABYCAYAAADskazQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAANO0lEQVR4nO3df5BV9XnH8fdHFjCAKAg6KJKFVNvSNiFxBWzGaYSJok2CndSKTic2mGptzKTTmlabmQbHtjGYjBnbTNQ2GbWZJKj5pWkqUpuaToYgixIERFmJUXQnoqhoaSDo0z/Os3Jc9xd379675/J5zZzZc59zvt/vee6B55793r33KCIwM7PqOaLZB2BmZrVxATczqygXcDOzinIBNzOrKBdwM7OKcgE3M6uoUVfAJX1C0mOStkhaOcB+YyQ9LOn7pdgKSc9I2pjLuRkfK+k2SY9IelTS1aU24yTdIulxSdskfTjj4yWtktQlaZ2k9lKbiyVtz+XiUnx27rs9247LuCTdmH1tkvSeUpslmW+XpKtK8amS1mRfayRNGcJzN0vSfZnj1vIxm1nraVoBl/Q+Sbf2ip0JLAXeGRG/BXx+gC4+CTzaR/yGiJiXyw8ydj4wPiJ+BzgVuKxU3D4NPBcRpwBzgQcyfgnwYkT8GnAD8Lk8xqnAZ4AFwHzgM6Xi+rkc/2TgxewD4Bzg5FwuBb6cfY0BvpTb5wIXSpqbba4C7s++7s/Hg7kduD4ifjOP7bkhtDGzihptV+CXA9dFxD6AiOizAEmaCfw+8K9D7DeAiZLagLcB+4E9uW058Nkc7/WIeD7jS4Hbcv0uYLEkAWcDayJid0S8CKwBluS2Rbkv2fa8Ul+3R+EnwDGSZlAU2a6I2BER+4Fv5r69x3+jr/zN43pJ6/Nq/rKMzwXaImJN5vJqROwd4vNjZhU02gr4KcAZOQ3xgKTT+tnvi8BfA6/3se2KLGxfLV0Z3wX8L9ANPAV8PiJ2Szomt18r6SFJd0o6PmMnAk8DRMQB4GXg2HI87czYscBLuW85zgBt+osDHB8R3Tl+N3Bcxi8BXo6I04DTgD+VNJviuXtJ0rdzaun6vMI3sxbV8AKexXkjxdXzh0rz1WcDbcAUYCHwKeCOvLItt/8AxZTHhj66/zLwDmAeRbH+QsbnA68BJwCzgb+SNCfHmwn8OCLeA6zl4LTNm8ZNUUO8lr4GchbwkXwO11G8cJxMkcsZwJUUhX0O8CeD9GVmFdbwAh4RCyJiHvAx4O7SfPVqiivQb+dUw4MUV9jTenXxXorC/yTFlMMiSV/Lvn8REa9FxOvAv1AUboCLgHsj4lc5LfNjoAN4AdgLfCf3uxPoeYNxJ3ASQE69HA3sLsfTTOBZ4HmKqZG2XnEGaNNfHOAXOc1C/uyZThLwidLzNjsi7su+Hs7pmAPAd0u5mFkLGm1TKN+lmEdG0inAOIrC+IaIuDoiZkZEO7AM+K+I+ONsM6O06x8Am3P9KYpCL0kTKa7wt0XxTV73AO/L/RYDW3P9bqDnL0z+MMcJYDVwlqQpOUVzFrA6t/0w9yXbfq/U10dy/IUUUyDdwHrg5PzrlXGZz919jF/uazVwuaSxPc9T5rQemCJpeu63qJSLmbWiiGjKQlE0b+0VGwd8jaLwPgQsyvgJwA/66eP7pcf/BjwCbKIogDMyPoni6noLRVH7VKnN24EfZZv7gVkZPzLbdAEPAnNKbZZnvAv4aCk+J/ftyrbjMy6KvzZ5Io+vo9TmXODx3PbpUvzYPJ7t+XNqxo8A/jH72UzxonF0bnt/5vEIcCswrlnn14sXLyO/KMJfJ2tmVkWjbQrFzMyGyAXczKyi2gbfpX6mTZsW7e3tjRzSzKzyNmzY8HxETO8db2gBb29vp7Ozs6a2F9y8ts5HM7it3XsG36mXaZPGc9xR4xs6rsccuTHrMe7hMmYt41Z1zFosnXciFy2YVVNbST/vK16JKZQLbl7Lup/tbuiYW7v38MovDwy+Y8nefQd4/tV9DR3XY47cmPUY93AZs5ZxqzpmLbZ27+F7G5+pe78NvQIfjgWzp7LqstMbNl7PFf+hjFlLm+H24TFHbsx6jHu4jFlLH1UdsxYjNYNQiStwMzN7KxdwM7OKcgE3M6soF3Azs4pyATczqygXcDOzinIBNzOrqCH9HXjePOEVirvaHIiIjry57yqgHXgS+KMo7hFpZmYNcChX4GdGcQeYjnxcy13TzcysToYzhdLnXdPNzKwxhlrAA7hP0gZJl2asv7umv4mkSyV1SurctWvX8I/YzMyAoX8Xynsj4llJxwFrJG0b6gARcQtwC0BHR4dv/2NmVidDugKPiGfz53MUd3CfT/93TTczswYYtIBLmijpqJ51iruwb6b/u6abmVkDDGUK5XjgO5J69v96RNwraT1wh6RLgKeA80fuMM3MrLdBC3hE7ADe1Uf8BWDxSByUmZkNzp/ENDOrKBdwM7OKcgE3M6soF3Azs4pyATczqygXcDOzinIBNzOrKBdwM7OKcgE3M6soF3Azs4pyATczqygXcDOzinIBNzOrKBdwM7OKcgE3M6soF3Azs4pyATczqygXcDOzihpWAZe0RNJjkrokXVWvgzIzs8HVXMAljQG+BJwDzAUulDS3XgdmZmYDG84V+HygKyJ2RMR+4JvA0voclpmZDWY4BfxE4OnS450ZMzOzBhhOAVcfsXjLTtKlkjolde7atWsYw5mZWZki3lJzh9ZQOh1YERFn5+OrASLis/216ejoiM7OzkMe64Kb17K1ew9zZ0yu6VhrsbV7D3v3v8aEcWOG3GbvvgNMGN82rOM81HE95siNWY9xD5cxaxm3qmPWoqd+rbrs9JraS9oQER2948O5Al8PnCxptqRxwDLg7mH016+l805s+BM+beK4Q/oPADBhfBvTJo1v6Lgec+TGrMe4h8uYtYxb1TFrMXfGZJbOq/8Mc1utDSPigKQrgNXAGOCrEbGlbkdWctGCWVy0YNZIdG1mVlk1T6HUNJi0C/h5nbqbBjxfp75Gm1bODVo7v1bODZxfs7w9Iqb3Dja0gNeTpM6+5oRaQSvnBq2dXyvnBs5vtPFH6c3MKsoF3MysoqpcwG9p9gGMoFbODVo7v1bODZzfqFLZOXAzs8Ndla/AzcwOaw0t4JKOlPSgpJ9K2iLpmoy/S9JaSY9IukfS5F7tZkl6VdKVpdiFuf8mSfdKmpbx8ZJW5VfcrpPUXmpzsaTtuVw8yvO7IHPbImllKV6Z/CS1S/o/SRtzuanU16m5f5ekGyWpmfnVObd/kPS0pFd7jVH5cydpgqR/l7Qt+7mu2fnV+dzdW+rnJhXfutrUczegiGjYQvH9KZNyfSywDlhI8anO38v4cuDaXu2+BdwJXJmP24DngGn5eCXFx/oB/hy4KdeXAatyfSqwI39OyfUpozS/Y4GngOn5+DZgcdXyA9qBzf309SBwevb5H8A5zcyvzrktBGYAr/aKV/7cAROAM3N9HPA/LXbuJpf6/BawrNnnbqCloVfgUei5KhmbSwC/Dvwo42uAD/e0kXQexZNS/pSncpmYV26TgWdz21KKggdwF7A49zkbWBMRuyPixRxnySjNbw7weET0fPvXf5baVCq/vkiaQfEfZW0U/wtuB87LzU3Jr165ZV8/iYjuPjZV/txFxN6I+GGu7wceAmbm5lY4d3tytY3iBarnTcKmnbuBNHwOXNIYSRsprqDXRMQ6YDPwodzlfOCk3Hci8DfANeU+IuJXwOXAIxSFey7wldz8xtfcRsQB4GWKK9qGfP1tPfIDuoDfyF/12iiK20m5rTL5pdmSHpb0gKQzSjns7OdYm5ZfnXIbSCucu3J/xwAfBO7PUEucO0mrs59XKIo1NPnc9afhBTwiXouIeRSv2vMl/TbFrzcfl7QBOArYn7tfA9xQenUFQNJYigL+buAEYBNwdc/mvoYdIF5X9cgvX8kvB1ZR/Ir6JHAgN1cpv25gVkS8G/hL4Os5DznQsTYtvzrlNpBWOHcA5IXFN4AbI2JHT7ifPCp17qL4htUZwHhgUYabeu7607S/QomIl4D/BpZExLaIOCsiTqX4R/FE7rYAWCnpSeAvgL9V8QVa87KPJ/JX8DuA3802Ozl4hdsGHA3sLsfTTA5Ou9TdMPMjIu6JiAURcTrwGLC9avlFxL6IeCHXN2T8lDzWmaXuysfa9PyGmdtAmp4b1C2/W4DtEfHFUqzp+dXr3EXELym+XbXnLmNNz61P0aDJ9qLOMh04JtffRnF1+QHguIwdQTEfuryPtis4+CbfCRSvoj1v8l0LfCHXP86b32y4Iw6+2fAzijcapuT61NGYXz7uaTMF2AicUrX8cv8xuT4HeKbnmCjeYFrIwTcxz21mfvXMrdRn7zcxW+Xc/T3FG3xHjIb86pUbMAmYkfE2it+Ar2j2uRsw90YNlMm+E3iYYspjM/B3Gf8k8Hgu15EfMOrVdgVvLnB/Bjyafd0DHJvxIyn+oqOL4i8d5pTaLM94F/DRUZ7fN4CtuSwrxSuTH8WbRluAn1K82fXBUl8d2ccTwD+X2jQlvzrntpLiyuz1/LmiVc4dxRVmUPzf25jLx1rh3AHHU1xYbMrt/wS0NfvcDbT4k5hmZhXlT2KamVWUC7iZWUW5gJuZVZQLuJlZRbmAm5lVlAu4mVlFuYCbmVWUC7iZWUX9P878DdU17LE8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot a small number of polygon bounds\n",
    "polygon_tiles_gpd.head(10).boundary.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "1c5a3b95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/data/Bike paths/inf tiles_polygons.gpkg'"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polygon_export_loc = output_folder+'_polygons.gpkg'\n",
    "polygon_export_loc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "8fd9b7a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon_tiles_gpd.to_file(polygon_export_loc, driver=\"GPKG\")    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "d366200c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a list of which tiles are within which geotiffs\n",
    "def intersector(geo_tiff):\n",
    "    tiles_inside_geo_tiff = []\n",
    "#     loop over each tile and check if the geotiff is the the intersect list\n",
    "    for tile in tile_polygon_list:\n",
    "        if geo_tiff in tile['intersect']:\n",
    "#             count this so we know if the tile will be incomplete or not\n",
    "            incomplete = len(tile['intersect'])>1\n",
    "#             build dict with geom the current row and col for naming\n",
    "            tiles_inside_geo_tiff.append({'geometry':tile['geometry'],'row':tile['row'],'col':tile['col'],'name':tile['name'],'incomplete':incomplete})\n",
    "    return([geo_tiff,tiles_inside_geo_tiff])  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "803304c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f9791650132449d299e99cc530ac0149",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with Pool() as pool:\n",
    "    geo_tiff_with_tiles = list(tqdm(pool.imap(intersector, geo_tiff_list), total=len(geo_tiff_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "de06b782",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# cut tiles from rasters\n",
    "def cut_tiles(geotiff):\n",
    "#     grab path to to file and open it\n",
    "    geotiff_open = gdal.Open(geotiff[0])\n",
    "#     grab the filename and strip the extension\n",
    "    geo_tiff_filename = os.path.basename(geotiff[0]).replace(input_file_ext,'')\n",
    "    incomplete_tile_list = []\n",
    "    for tile in geotiff[1]:\n",
    "        tile_geometry = tile['geometry']\n",
    "#         shapely bounds returns \"minx, miny, maxx, maxy\" but we need minx, maxy, maxx, miny\n",
    "        top = list(tile_geometry.bounds)[3]\n",
    "        bottom = list(tile_geometry.bounds)[1]\n",
    "        left = list(tile_geometry.bounds)[0]\n",
    "        right =list(tile_geometry.bounds)[2]\n",
    "        \n",
    "#         make row folder path\n",
    "        if output_to_row_folders:\n",
    "            output_row_folder = os.path.join(output_folder,str(tile['row']))\n",
    "            Path(output_row_folder).mkdir(parents=True, exist_ok=True)\n",
    "        else:\n",
    "            output_row_folder = output_folder\n",
    "#       make row folder if necessary\n",
    "        \n",
    "        export_file_name = str(tile['name'])#str(tile['row'])+'_'+str(tile['col'])+'.tif'\n",
    "        \n",
    "#         check if tile is incomplete if so append the getiff name so that it is unique\n",
    "        if tile['incomplete']:\n",
    "            append_name = '-'+geo_tiff_filename+'_incomplete.tif'\n",
    "            export_file_name = export_file_name.replace('.tif',append_name)\n",
    "#             add tile to list so we dont need to refind them to compile incomplete tiles\n",
    "            export_file_path = os.path.join(output_row_folder,export_file_name)\n",
    "            incomplete_tile_list.append(export_file_path)\n",
    "        else:\n",
    "            export_file_path = os.path.join(output_row_folder,export_file_name)\n",
    "        \n",
    "#         check if already done\n",
    "        if not os.path.isfile(export_file_path):\n",
    "\n",
    "    #     clip the data\n",
    "    #         make a string of tile dims to pass as a command line arg, this is kinda of hacky, would like a better option\n",
    "            tile_clip_string = str(left) +' '+str(top) +' '+str(right) +' '+str(bottom)\n",
    "\n",
    "            translate_options = gdal.TranslateOptions(gdal.ParseCommandLine(\"-projwin \"+tile_clip_string)\n",
    "                                                     ,creationOptions=['COMPRESS='+output_compression])\n",
    "\n",
    "            tile_clip = gdal.Translate(export_file_path, geotiff_open, options = translate_options)\n",
    "    #     close the tile\n",
    "            tile_clip = None\n",
    "    return incomplete_tile_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "2f07cdb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# used to slip a list into n chunks of roughly the same size\n",
    "def split(a, n):\n",
    "    k, m = divmod(len(a), n)\n",
    "    return (a[i*k+min(i, m):(i+1)*k+min(i+1, m)] for i in range(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "64bb740c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set this to cpu_count() unless you are running out of RAM in which case decrease it\n",
    "concurrent_threads = cpu_count()\n",
    "concurrent_threads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "b34ebdf6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "360"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tile_count = 0\n",
    "for raster,tiles in geo_tiff_with_tiles:\n",
    "    tile_count+=len(tiles)\n",
    "average_tile_count = round(tile_count/len(geo_tiff_with_tiles))\n",
    "average_tile_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "fbc498ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You have more threads than rasters so we each thread will get the same raster\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "18d3ff05e18c41748306721fb7981180",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# if you have more cores than rasters and more cores than average tile count then give every core a copy of the same raster to work on\n",
    "if concurrent_threads > len(geo_tiff_with_tiles) and concurrent_threads < average_tile_count:\n",
    "    print('You have more threads than rasters so we each thread will get the same raster')\n",
    "    incomplete_tile_list = []\n",
    "\n",
    "    for raster,tiles in tqdm(geo_tiff_with_tiles):\n",
    "        geo_tiff_with_tile_list = []\n",
    "\n",
    "        for tile_split in split(tiles,concurrent_threads):\n",
    "            geo_tiff_with_tile_list.append([raster,tile_split])\n",
    "\n",
    "        pool = Pool(concurrent_threads)\n",
    "        with pool:\n",
    "            one_raster_incomplete_tile_list = list(pool.imap(cut_tiles,geo_tiff_with_tile_list))\n",
    "\n",
    "        incomplete_tile_list.append(np.concatenate(one_raster_incomplete_tile_list))\n",
    "    \n",
    "else:\n",
    "    print('you have more rasters than threads so each thread will get its own raster')\n",
    "    pool = Pool()\n",
    "    with pool:\n",
    "        incomplete_tile_list = list(tqdm(pool.imap(cut_tiles,geo_tiff_with_tiles), total=len(geo_tiff_with_tiles)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "ca407735",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 incomplete tiles\n"
     ]
    }
   ],
   "source": [
    "#  the pool above returns all the incomplete tiles as a list of list, we need a flat list\n",
    "flat_incomplete_tile_list = np.concatenate(incomplete_tile_list).ravel()\n",
    "print(len(flat_incomplete_tile_list),'incomplete tiles')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "c387973c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from the flat list make a new list with just rows and cols to filter by\n",
    "incomplete_tile_file_names = []\n",
    "for incomplete_tile in flat_incomplete_tile_list:\n",
    "    incomplete_tile_file_names.append(os.path.basename(incomplete_tile).split('-')[0].replace(raster_name_append,''))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "19367e0f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>incomplete_tiles</th>\n",
       "      <th>row_col</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [incomplete_tiles, row_col]\n",
       "Index: []"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make a new dataframe with paths and row/cols\n",
    "incomplete_tile_df = pd.DataFrame(\n",
    "    {'incomplete_tiles': flat_incomplete_tile_list,\n",
    "     'row_col': incomplete_tile_file_names\n",
    "    })\n",
    "\n",
    "unique_tiles_list = incomplete_tile_df.row_col.unique()\n",
    "\n",
    "incomplete_tile_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "13279c29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "incomplete_tile_df['incomplete_tiles'].to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "8da67c28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# handy gdal tool to merge files\n",
    "def merge_tiles(merge_imgs, output_path):\n",
    "    merge_command = ['gdal_merge.py', '-o', output_path, '-quiet', '-co','COMPRESS='+output_compression]\n",
    "    for name in merge_imgs:\n",
    "        merge_command.append(name)\n",
    "    subprocess.run(merge_command)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "90092fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def join_incomplete_tile(unique_tiles):\n",
    "#     filter the dataframe down to only one row and col\n",
    "    df_filt = incomplete_tile_df[incomplete_tile_df['row_col']==unique_tiles]\n",
    "#     get paths as list\n",
    "    combine_these = df_filt['incomplete_tiles'].tolist()\n",
    "#     get export name by removing the geotiff names from end\n",
    "    export_file_name = os.path.basename(combine_these[0]).split('-')[0]+'.tif'\n",
    "#     grab the folder from the first tile\n",
    "    export_dir_path = os.path.dirname(combine_these[0])\n",
    "#     build full export path\n",
    "    export_full_path = os.path.join(export_dir_path,export_file_name)\n",
    "#     use gdal gdal_merge.py to merge the tiles\n",
    "    merge_tiles(combine_these,export_full_path)\n",
    "#     remove the incomplete tiles and msk files\n",
    "    for incomplete_tile in combine_these:\n",
    "        try:\n",
    "            os.remove(incomplete_tile)\n",
    "        except:\n",
    "            print('could not remove')\n",
    "        if os.path.isfile(incomplete_tile+'.msk'):\n",
    "            try:\n",
    "                os.remove(incomplete_tile+'.msk')\n",
    "            except:\n",
    "                print('could not remove')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "e29d8828",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b03d05804c9424193aa320fde457fe9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "0it [00:00, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pool = Pool()\n",
    "with pool:\n",
    "    list(tqdm(pool.imap(join_incomplete_tile,unique_tiles_list), total=len(unique_tiles_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "22a039ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# if you want you can validate the output of the above code with the below cells, this is probibly not needed in most cases\n",
    "# from https://github.com/fastai/fastai/blob/master/fastai/vision/utils.py\n",
    "from PIL import Image\n",
    " \n",
    "def verify_image(fn):\n",
    "    try:\n",
    "        im = Image.open(fn)\n",
    "        im.draft(im.mode, (32,32))\n",
    "        im.load()\n",
    "        return True\n",
    "    except: return fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "b6e5ca2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "360"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_tiles = []\n",
    "for root, dirs, files in os.walk(output_folder):\n",
    "    for file in files:\n",
    "        if file.endswith('.tif'):\n",
    "            all_tiles.append(os.path.join(root, file))\n",
    "len(all_tiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "7adc632a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2057804f1ffc4a909e4ed781f475a8ec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/360 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with Pool():\n",
    "    failed_tiles = list(tqdm(Pool().imap(verify_image,all_tiles), total=len(all_tiles)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "c00a5a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in set(failed_tiles):\n",
    "    if not i:\n",
    "        print(i)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}