From 7c9ded2b17ea1453dae5beb64ec4a80f3d183936 Mon Sep 17 00:00:00 2001 From: uceshaf Date: Wed, 8 Jan 2025 04:27:37 +0500 Subject: [PATCH] Expand grid functionality and add demo notebook --- setup.py | 2 + src/extras/grid_functionality_demo.ipynb | 1016 ++++++++++++++++++++++ src/grid.py | 250 +++++- 3 files changed, 1264 insertions(+), 4 deletions(-) create mode 100644 src/extras/grid_functionality_demo.ipynb diff --git a/setup.py b/setup.py index 55b9e2b..72296db 100644 --- a/setup.py +++ b/setup.py @@ -11,6 +11,8 @@ install_requires=[ "torch>=1.10.0", "torchvision", + "torchgeo", # Added New + "open_clip_torch", # Added New "pandas", "geopandas", "numpy", diff --git a/src/extras/grid_functionality_demo.ipynb b/src/extras/grid_functionality_demo.ipynb new file mode 100644 index 0000000..57968cd --- /dev/null +++ b/src/extras/grid_functionality_demo.ipynb @@ -0,0 +1,1016 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# # Grid System Added Functionality Demonstration\n", + "# This notebook demonstrates the usage and functionality of the updated `grid.py` script." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# ## Imports and Setup\n", + "# Import required modules and initialize the grid.\n", + "\n", + "import os\n", + "import geopandas as gpd\n", + "import matplotlib.pyplot as plt\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), '..')))\n", + "\n", + "# Import the grid module\n", + "from grid import * #Grid, merge_utm_files_to_wgs84, ensure_file_extension" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# provide a folder path to the output folder\n", + "# ### Create Output Directory\n", + "# Set up directories to save output files.\n", + "output_folder = \"test_output\"\n", + "os.makedirs(output_folder, exist_ok=True)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid initialized with resolution: 10 km\n" + ] + } + ], + "source": [ + "\n", + "# ### Initialize Grid\n", + "# Create a `Grid` object with a 10 km grid resolution.\n", + "dist = 10 # Grid resolution in kilometers\n", + "grid = Grid(dist)\n", + "print(\"Grid initialized with resolution:\", dist, \"km\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid points (sample):\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
grid_cellrowcolrow_idxcol_idxutm_zoneutm_crsgeometry
0946D_176L946D176L0032701EPSG:32701POINT (-180.00000 -84.97006)
1946D_175L946D175L0132701EPSG:32701POINT (-178.97727 -84.97006)
2946D_174L946D174L0232701EPSG:32701POINT (-177.95455 -84.97006)
3946D_173L946D173L0332701EPSG:32701POINT (-176.93182 -84.97006)
4946D_172L946D172L0432701EPSG:32701POINT (-175.90909 -84.97006)
\n", + "
" + ], + "text/plain": [ + " grid_cell row col row_idx col_idx utm_zone utm_crs \\\n", + "0 946D_176L 946D 176L 0 0 32701 EPSG:32701 \n", + "1 946D_175L 946D 175L 0 1 32701 EPSG:32701 \n", + "2 946D_174L 946D 174L 0 2 32701 EPSG:32701 \n", + "3 946D_173L 946D 173L 0 3 32701 EPSG:32701 \n", + "4 946D_172L 946D 172L 0 4 32701 EPSG:32701 \n", + "\n", + " geometry \n", + "0 POINT (-180.00000 -84.97006) \n", + "1 POINT (-178.97727 -84.97006) \n", + "2 POINT (-177.95455 -84.97006) \n", + "3 POINT (-176.93182 -84.97006) \n", + "4 POINT (-175.90909 -84.97006) " + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# ## Generate Grid and Points\n", + "# Display the first few rows of the grid points.\n", + "grid_points = grid.points\n", + "print(\"Grid points (sample):\")\n", + "grid_points.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Footprint for grid cell 946D_176L:\n", + "POLYGON ((480986.97927458223 562504.1200605556, 480986.97927458223 573184.1200605556, 470306.97927458223 573184.1200605556, 470306.97927458223 562504.1200605556, 480986.97927458223 562504.1200605556))\n", + "Subset of grid points with footprints:\n", + " grid_cell utm_footprint\n", + "2389 322D_635R POLYGON ((721428.378041 6798369.184628, 721428...\n", + "2902 0U_898R POLYGON ((2507401.749273 9999660, 2507401.7492...\n", + "85 632U_1013L POLYGON ((3321495.912167 22592161.866856, 3321...\n", + "665 121U_1303L POLYGON ((748303.70546 28793323.915129, 748303...\n", + "782 15U_1221L POLYGON ((-304183.851526 29845446.444146, -304...\n", + "3436 169U_1502R POLYGON ((11557199.604514 15606109.726045, 115...\n", + "339 449D_1189L POLYGON ((2477258.195161 -5267955.752287, 2477...\n", + "273 188U_1644L POLYGON ((4738947.216293 27681556.244205, 4738...\n", + "2582 118U_612R POLYGON ((-266123.978933 11180138.349109, -266...\n", + "2802 59D_807R POLYGON ((1600783.111277 9405228.375646, 16007...\n" + ] + } + ], + "source": [ + "# ## Test Footprint for a Single Grid Cell\n", + "# Demonstrate how to retrieve the footprint for a single grid cell.\n", + "sample_grid_cell = grid.points['grid_cell'].iloc[0] # Use the first grid cell\n", + "footprint = grid.get_product_outline_for_cell(grid_cell=sample_grid_cell)\n", + "\n", + "print(f\"Footprint for grid cell {sample_grid_cell}:\")\n", + "print(footprint)\n", + "\n", + "# ## Generate Footprints for Multiple Grid Cells\n", + "# Add a `utm_footprint` column to a subset of the grid points.\n", + "subset_gdf = grid.points.sample(10) # Take a random sample of grid points\n", + "subset_with_footprints = generate_product_outlines(\n", + " subset_gdf, get_footprints=True\n", + ")\n", + "\n", + "print(\"Subset of grid points with footprints:\")\n", + "print(subset_with_footprints[[\"grid_cell\", \"utm_footprint\"]])" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# For each UTM Zone, we can run the product outline creation and footprint generation process." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['EPSG:32701' 'EPSG:32702' 'EPSG:32703' 'EPSG:32704' 'EPSG:32705'\n", + " 'EPSG:32706' 'EPSG:32707' 'EPSG:32708' 'EPSG:32709' 'EPSG:32710'\n", + " 'EPSG:32711' 'EPSG:32712' 'EPSG:32713' 'EPSG:32714' 'EPSG:32715'\n", + " 'EPSG:32716' 'EPSG:32717' 'EPSG:32718' 'EPSG:32719' 'EPSG:32720'\n", + " 'EPSG:32721' 'EPSG:32722' 'EPSG:32723' 'EPSG:32724' 'EPSG:32725'\n", + " 'EPSG:32726' 'EPSG:32727' 'EPSG:32728' 'EPSG:32729' 'EPSG:32730'\n", + " 'EPSG:32731' 'EPSG:32732' 'EPSG:32733' 'EPSG:32734' 'EPSG:32735'\n", + " 'EPSG:32736' 'EPSG:32737' 'EPSG:32738' 'EPSG:32739' 'EPSG:32740'\n", + " 'EPSG:32741' 'EPSG:32742' 'EPSG:32743' 'EPSG:32744' 'EPSG:32745'\n", + " 'EPSG:32746' 'EPSG:32747' 'EPSG:32748' 'EPSG:32749' 'EPSG:32750'\n", + " 'EPSG:32751' 'EPSG:32752' 'EPSG:32753' 'EPSG:32754' 'EPSG:32755'\n", + " 'EPSG:32756' 'EPSG:32757' 'EPSG:32758' 'EPSG:32759' 'EPSG:32760'\n", + " 'EPSG:32601' 'EPSG:32602' 'EPSG:32603' 'EPSG:32604' 'EPSG:32605'\n", + " 'EPSG:32606' 'EPSG:32607' 'EPSG:32608' 'EPSG:32609' 'EPSG:32610'\n", + " 'EPSG:32611' 'EPSG:32612' 'EPSG:32613' 'EPSG:32614' 'EPSG:32615'\n", + " 'EPSG:32616' 'EPSG:32617' 'EPSG:32618' 'EPSG:32619' 'EPSG:32620'\n", + " 'EPSG:32621' 'EPSG:32622' 'EPSG:32623' 'EPSG:32624' 'EPSG:32625'\n", + " 'EPSG:32626' 'EPSG:32627' 'EPSG:32628' 'EPSG:32629' 'EPSG:32630'\n", + " 'EPSG:32631' 'EPSG:32632' 'EPSG:32633' 'EPSG:32634' 'EPSG:32635'\n", + " 'EPSG:32636' 'EPSG:32637' 'EPSG:32638' 'EPSG:32639' 'EPSG:32640'\n", + " 'EPSG:32641' 'EPSG:32642' 'EPSG:32643' 'EPSG:32644' 'EPSG:32645'\n", + " 'EPSG:32646' 'EPSG:32647' 'EPSG:32648' 'EPSG:32649' 'EPSG:32650'\n", + " 'EPSG:32651' 'EPSG:32652' 'EPSG:32653' 'EPSG:32654' 'EPSG:32655'\n", + " 'EPSG:32656' 'EPSG:32657' 'EPSG:32658' 'EPSG:32659' 'EPSG:32660']\n" + ] + } + ], + "source": [ + "print(grid.points['utm_crs'].unique())" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using UTM Zone: EPSG:32633\n" + ] + } + ], + "source": [ + "utm_zone = grid.points['utm_crs'].unique()[92] # Use any available UTM zone\n", + "print(\"Using UTM Zone:\", utm_zone)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "shift=340 # Sentinel-2 grid extended on each side to create 1068*1068 raster\n", + "pixel_size=10 # Sentinel-2 pixel size in meters\n", + "raster_width=1068 # Sentinel-2 tile size in Major TOM\n", + "raster_height=1068 # Sentinel-2 tile size in Major TOM\n", + "driver=\"ESRI Shapefile\" # Use \"GeoJSON\" or \"GPKG\"\n", + "get_footprints=False # Set to True to get utm_footprints in UTM coordinates" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using UTM Zone: EPSG:32633\n", + "Product outlines generated for EPSG:32633 and saved to test_output\\product_outlines_sample_utmzone.shp\n" + ] + } + ], + "source": [ + "# ## Generate Product Outlines for a UTM Zone\n", + "# Choose a UTM zone and generate product outlines. Default values for Sentinel-2 L2A products are used.\n", + "# The product outlines are saved to a shapefile.\n", + "utm_zone = grid.points['utm_crs'].unique()[92] # Use the first available UTM zone\n", + "print(\"Using UTM Zone:\", utm_zone)\n", + "\n", + "product_outlines_file = os.path.join(output_folder, \"product_outlines_sample_utmzone.shp\")\n", + "\n", + "# generate_product_outlines_for_utm_zone(self, utm_zone, shift=340, pixel_size=10, raster_width=1068, raster_height=1068, output_file=None, driver=\"ESRI Shapefile\", get_footprints=False)\n", + "\n", + "\n", + "product_outlines = grid.generate_product_outlines_for_utm_zone(\n", + " utm_zone=utm_zone,\n", + " output_file=product_outlines_file,\n", + " driver=\"ESRI Shapefile\",\n", + " get_footprints=False\n", + ")\n", + "print(f\"Product outlines generated for {utm_zone} and saved to {product_outlines_file}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAANRCAYAAABOfBGKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIjUlEQVR4nO3deXxTVcI+8Cdt2nRvgUI3SktLoewwgLIvKm1ZVUZHGVFAmBlZ3ECcATfwBSqMMvg6guMCDCKgzIAiQmlZBNkUEGX5AUa0FoFSqNAC1dLl/P7om0BImnRJcm5ynu/n04+S3p77FMpDzs3JuTohhAARkYb4yA5ARHQrFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItKcWhXTsmXLoNPpzB96vR5NmzbF2LFjcebMGVdltKDT6TBz5kyXjb9o0SIsW7asVl9z7do1vPLKK+jcuTNCQkIQHByMTp06Ye7cubh27Vqds+zZswczZ87E5cuXrT7Xv39/9O/f3+IxV//e1Ncbb7yBFi1awN/fHzqdzub35Sy3/qze+vH555+bj01MTLT4XEhICG6//XYsX77cYsyysjL861//Qrdu3dCwYUMEBQUhISEBd999N9atW2eV4fz585gxYwY6deqEsLAw+Pv7o2nTphgxYgTWr1+PiooKh9/Hli1bMHDgQMTGxsJgMKBJkya44447sHHjRovjiouLMWfOHPTv3x/R0dEICQlB+/btMW/ePPz22282xz569Cjuv/9+NG7cGAaDAYmJiZg4caLFMatWrULfvn0RFRUFg8GA2NhYDBs2DHv27LEab/z48WjXrh0iIiIQGBiIli1bYtq0abh48aLD7/NW+lp/BYClS5ciNTUVv/76K3bu3InMzEzs2LEDR44cQXBwcF2G1IxFixYhMjISY8aMqdHx58+fx1133YVTp07hiSeewPz58wEA27Ztw+zZs7Fq1Sps2bIFUVFRtc6yZ88ezJo1C2PGjEFERITD4/fu3YumTZvW+jzu8M033+CJJ57A+PHjMXr0aOj1eoSGhrr8vKaf1Vu1adPG4te9evXCq6++CgD4+eef8eqrr2L06NG4du0aJkyYAAB4+OGHsXbtWjz11FOYNWsWDAYDfvjhB2RlZWHz5s249957zePt27cPw4cPhxACEyZMQPfu3RESEoK8vDx8+umnGDFiBP71r39h3LhxdvMXFhaibdu2GD9+PKKjo/HLL7/grbfewpAhQ/D+++9j1KhRAIC8vDwsXLgQDz/8MKZMmYKQkBB88cUXmDlzJnJycpCTkwOdTmced/v27RgyZAj69OmDt956C5GRkcjLy8OhQ4eszt+rVy88+eSTiIyMxLlz57BgwQL07dsXW7duRb9+/czHXrt2DX/+85/RokULBAQE4MCBA5gzZw42btyIQ4cOwd/fvyZ/ZFVELSxdulQAEPv377d4/IUXXhAAxIoVK6r92mvXrtXmVNUCIF566SWnjGVL27ZtRb9+/Wp8fFpamtDr9eKLL76w+twXX3wh9Hq9SE9Pr1OWv//97wKA+PHHH60+169fv1rllG3FihUCgPjyyy+dNqa9n6nqflZtSUhIEEOGDLF47NKlSyIsLEy0aNFCCCHEDz/8IACIF1980eYYFRUVFl8bFRUlmjdvLs6ePWvz+G+//VZs27bNYTZbrl+/LuLi4kSfPn3Mj129elVcvXrV6ljTz9DNP5/Xrl0TMTExYsiQIaKysrLW5798+bLw8/MTDz/8sMNjFy1aJACIrVu31uocTrnG1L17dwDATz/9BAAYM2YMQkJCcOTIEaSlpSE0NBR33nknAOCXX37BxIkTERcXB39/fyQlJeG5555DaWmpxZjFxcX405/+hEaNGiEkJAQZGRn47rvvrM49ZswYJCYmWj0+c+ZMi38hAKCyshJvvPEGOnXqhMDAQERERKB79+5Yv349gKqn9MeOHcOOHTvMT+ttjW1y4MABZGdnY9y4cejdu7fV53v37o1HH30UmzdvxsGDBwEAubm50Ol0NqeLN0/FZs6ciWnTpgEAmjdvbnMKYu/rgRvTme3bt2PChAmIjIxEo0aNMGLECJw9e9bq6z/88EP06NEDwcHBCAkJQXp6utW/oD/88AMefPBB89QiKioKd955J7755ptqc/Xv39/8L/vtt98OnU5n8Yx0yZIl6NixIwICAtCwYUPce++9OH78uMUY9n6mXCEiIgKtWrUy/0wXFhYCAGJiYmwe7+Nz46/SO++8g/Pnz2P+/PnVHt+hQwcMGDCgTtn8/PwQEREBvf7GhCc4ONjmbOW2224DAJw+fdr82Jo1a3Du3DlMmzbN6u9ITYSGhiIgIMDi/NVp3LgxANTo2Js5pZi+//57ixAAcP36dQwfPhx33HEHPvnkE8yaNQu//fYbBgwYgOXLl2PKlCn47LPPMGrUKMyfPx8jRowwf60QAvfccw/ef/99TJ06FevWrUP37t0xaNCgeuUcM2YMnnzySXTr1g0ffvghVq9ejeHDhyM3NxcAsG7dOiQlJaFz587Yu3cv9u7da/PagUlOTg4A4J577qn2GNPnTMfW1Pjx4/H4448DANauXWvO87vf/a5W45jG8vPzw8qVKzF//nx8/vnn5qIwmTt3LkaOHIk2bdrgo48+wvvvv48rV66gT58++H//7/+Zjxs8eDAOHjyI+fPnIycnB4sXL0bnzp3tXi9atGgRnn/+eQBVU6u9e/fihRdeAABkZmZi3LhxaNu2LdauXYvXX38dhw8fRo8ePWA0Gi3GsfUz5UhFRQXKy8stPmpybaesrAw//fST+We6devWiIiIwKxZs/D222+bf2ZsycnJga+vLwYPHuzwPDdLTEys9h/CyspKlJeX4+zZs3jppZfw3XffYerUqQ7H3LZtGwCgbdu25sd27twJoOr3pnfv3vD390eDBg0wcuRIm/9gmY4tKytDbm4uJkyYACEEJk2aZPPY8vJyXLt2Dbt378YLL7yA3r17o1evXg6zWqjN0yvT0+N9+/aJsrIyceXKFbFhwwbRuHFjERoaKvLz84UQQowePVoAEEuWLLH4+rfeeksAEB999JHF4/PmzRMARHZ2thBCiE2bNgkA4vXXX7c4bs6cOVZTudGjR4uEhASrrC+99JK4+dvbuXOnACCee+45u99jbaZyjz32mAAgTpw4Ue0xx48fFwDEhAkThBBC/PjjjwKAWLp0qdWxt35vtZ3K3fr1pj+viRMnWhw3f/58AUCcO3dOCCFEXl6e0Ov14vHHH7c47sqVKyI6Olr84Q9/EEIIcfHiRQFALFy4sNrvtzq2plaXLl0SgYGBYvDgwRbH5uXlCYPBIP74xz+aH6vuZ8rR+Wx9+Pr6WhybkJAgBg8eLMrKykRZWZn48ccfzeebNm2a+bjPPvtMREZGmsdp1KiRuP/++8X69estxktNTRXR0dFWmSoqKsznKCsrs5j+CSFEcnKySE5Otvn9pKenm88bFhYm1q5d6/D34NtvvxWBgYHi3nvvtTlWRESEePbZZ8W2bdvEW2+9JRo1aiRatGhhc4rcqlUr8/ljYmLErl27bJ5z7969Fr/XgwcPFsXFxQ6z3qpOz5i6d+8OPz8/hIaGYujQoYiOjsamTZusLvD+/ve/t/j1tm3bEBwcjPvuu8/icdPT+q1btwKoujAHAA899JDFcX/84x/rEhcAsGnTJgCotuVdRfzfdld1ecrsLMOHD7f4dYcOHQDcmHpv3rwZ5eXleOSRRyyeWQQEBKBfv37m6WPDhg2RnJyMv//971iwYAEOHTqEysrKOufau3cvfv31V6sXGuLj43HHHXeYfx5uduvPlCPLly/H/v37LT6+/PJLq+M2btwIPz8/+Pn5oXnz5vjoo4/w+OOPY/bs2eZjBg8ejLy8PKxbtw7PPPMM2rZti48//hjDhw/H5MmTHWaZMmWK+Rx+fn5Wfy7ff/+9efZxqzfeeANfffUVPvnkE6Snp+OBBx7AqlWrqj1Xbm4uhg4divj4eLz77rsWnzP9mT3wwAOYN28eBgwYgL/85S9477338P3332PlypVW4/33v//Fl19+iTVr1qBNmzYYNGiQzcsK7du3x/79+7Fjxw68/vrrOHToEAYOHIiSkhJHvz0W6vSq3PLly9G6dWvo9XpERUXZnEcHBQUhLCzM4rHCwkJER0db/SVt0qQJ9Hq9eR5fWFgIvV6PRo0aWRwXHR1dl7gAgAsXLsDX17deY9yqWbNmAIAff/wRrVq1snmM6Sl/fHy8085bW7f+PhoMBgDAr7/+CqDqlUUA6Natm82vN10/0el02Lp1K15++WXMnz8fU6dORcOGDfHQQw9hzpw5tX6Vzd51m9jYWKvpr62fKUdat26Nrl27Ojyud+/e+Mc//gGdToegoCAkJyfbfBUpMDAQ99xzj3mKnpeXh0GDBuHNN9/EhAkT0LZtWzRr1gxGoxElJSUICgoyf+3UqVPNU+hbS8mRlJQU8/8PHz4cgwYNwqRJk/DAAw9YXN8Cqv7BGTBgAPR6PbZu3YqGDRtafN7085Cenm7xeHp6OnQ6Hb7++mur85umgrfddhvuuecedO7cGU8++SS+/fZbi+OCg4PNv999+/bF7bffju7du+Nf//oXnn766Rp/v3V6xmT6w+7UqVO1F/dsPUNo1KgRzp8/b34WYVJQUIDy8nJERkaajysvLzf/4Jrk5+dbjRkQEGB14RyA1dqJxo0bo6KiwuYYdTVw4EAAwMcff1ztMabPmY4NCAgAAKvMt36v7mT6ff/Pf/5j9ezi1mcYCQkJeO+995Cfn4+TJ0/i6aefxqJFi8wX6mvD9Bfk3LlzVp87e/asOZeJK591hoeHo2vXrujSpQtat25d45e2mzVrhj//+c8AgGPHjgGo+rOuqKiwWmsUHx+Prl27omvXrrV76dyG2267DZcuXcKFCxcsHv/pp5/Qv39/CCGwfft2m8tHTM+Yq3Nr0d1Kr9fjd7/7nc0Xo27VtWtX+Pj41OhYiwy1Orqe7rzzTly9etXqL7JpIZvpVRbTqxUffPCBxXG2nmImJiaioKDA/K8+UHWRdPPmzRbHmS6cL1682G5Gg8FgfibhSNeuXZGWlob33nsPu3fvtvr8rl27sGTJEmRkZKBLly4AgKioKAQEBODw4cMWx37yySc2swCocZ66Sk9Ph16vx6lTp8x/cW79sKVly5Z4/vnn0b59e5v/yjrSo0cPBAYGYsWKFRaP//zzz9i2bZtLX3WrrStXruDq1as2P2d6BTE2NhZA1YsNUVFRePbZZ22Wbn0JIbBjxw5ERERYPBvOy8tD//79UVFRgW3btiEhIcHm1997773Q6XTmyxsmmzZtghDC/Cp7dX777Tfs27cPLVq0cJh1x44dqKysrNGxN6vTVK6uHnnkEbz55psYPXo0cnNz0b59e+zatQtz587F4MGDcddddwEA0tLS0LdvXzz77LO4du0aunbtit27d+P999+3GvOBBx7Aiy++iAcffBDTpk3Db7/9hv/93/+1euWlT58+ePjhhzF79mycP38eQ4cOhcFgwKFDhxAUFGR+Bax9+/ZYvXo1PvzwQyQlJSEgIADt27ev9ntavnw57rrrLqSlpeGJJ54w/2Xatm0bXn/9daSmplosDdDpdBg1ahSWLFmC5ORkdOzYEV999ZXN0jWd9/XXX8fo0aPh5+eHVq1aOX1hYmJiIl5++WU899xz+OGHH5CRkYEGDRrg/Pnz+OqrrxAcHIxZs2bh8OHDmDx5Mu6//36kpKTA398f27Ztw+HDh/G3v/2t1ueNiIjACy+8gBkzZuCRRx7ByJEjUVhYiFmzZiEgIAAvvfRSvb+3o0ePory83Orx5ORki1eRHTl58iTS09Px4IMPol+/foiJicGlS5fw2Wef4e2330b//v3Rs2dP8/f18ccfY9iwYejYsaPFAsvCwkLs3LkT+fn55uNNTH95b77OdPfdd6Njx47o1KkTGjVqhLNnz2LZsmXYsWMH3nzzTfPL8AUFBRgwYADOnTuH9957DwUFBSgoKDCP07RpU/Ozp9TUVEyaNAmLFi1CaGgoBg0ahO+++w7PP/88OnfujD/84Q/mr+vZsyeGDx+O1q1bIzw8HLm5uVi8eDFOnTpl8Yr1hg0b8M4772D48OFISEhAWVkZDhw4gIULF6JFixYYP358jX+vAThngeWtRo8eLYKDg21+rrCwUDz22GMiJiZG6PV6kZCQIKZPny5+++03i+MuX74sHn30URERESGCgoLEwIEDxYkTJ2wusNy4caPo1KmTCAwMFElJSeKf//yn1atyQlS9KvKPf/xDtGvXTvj7+4vw8HDRo0cP8emnn5qPyc3NFWlpaSI0NFQAsPmK362uXr0q5s6dKzp16iSCgoJEUFCQ6NChg5g9e7bNRW9FRUVi/PjxIioqSgQHB4thw4aJ3Nxcm9/b9OnTRWxsrPDx8REAxPbt24UQtXtV7tY/r+3bt1uMZfLxxx+LAQMGiLCwMGEwGERCQoK47777xJYtW4QQQpw/f16MGTNGpKamiuDgYBESEiI6dOgg/vGPf4jy8nK7v0f2fnbeffdd0aFDB/Ofyd133y2OHTtmcYy9nyl756vu45133jEfa2uB5a0uXbokZs+eLe644w4RFxcn/P39RXBwsOjUqZOYPXu2KCkpsfqa/Px8MX36dNGhQwcRHBws/Pz8RGxsrBg2bJhYvny5KCsrszg+ISHB6udt3rx5olu3bqJBgwbC19dXNGrUSKSnp4sNGzZYHGf6M63u49afq/LycvHKK6+IFi1aCD8/PxETEyMmTJggLl26ZHHc1KlTRceOHUV4eLjQ6/UiOjpa3HvvvWL37t0Wxx0/flzcd999IiEhQQQEBIiAgACRmpoqpk2bJgoLC+3+3tqiE4J3SSEibeHuAkSkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTArYuXMnhg0bhtjYWOh0OrtvOq6OEAKvvvoqWrZsCYPBgPj4eMydO9f5YYng5vfKkRzXrl1Dx44dMXbs2FrvZ2Ty5JNPIjs7G6+++irat2+PoqKiOt39gqgm+JYUxeh0Oqxbt85iO+Dr16/j+eefxwcffIDLly+jXbt2mDdvnvn2UMePH0eHDh1w9OjRavedInImTuUIY8eOxe7du7F69WocPnwY999/PzIyMsx7bn/66adISkrChg0b0Lx5cyQmJmL8+PH45ZdfJCcnb8ViUtypU6ewatUqrFmzBn369EFycjKeeeYZ9O7dG0uXLgVQdWeUn376CWvWrMHy5cuxbNkyHDx40GqLZCJn4TUmxX399dcQQqBly5YWj5eWlpo3IausrERpaSmWL19uPu69995Dly5dcPLkSU7vyOlYTIqrrKyEr68vDh48CF9fX4vPhYSEAKjak1uv11uUV+vWrQFU7ZrIYiJnYzEprnPnzqioqEBBQQH69Olj85hevXqhvLwcp06dQnJyMgCY93CubvtWovrgq3IKuHr1qnm71s6dO2PBggUYMGAAGjZsiGbNmmHUqFHYvXs3XnvtNXTu3BkXL17Etm3b0L59ewwePBiVlZXo1q0bQkJCsHDhQlRWVmLSpEkICwtDdna25O+OvFKt97wkj1PdtqujR48WQghx/fp18eKLL4rExETh5+dn3j718OHD5jHOnDkjRowYIUJCQkRUVJQYM2ZMnbZMJaoJPmMiIs3hcgEi0hwWExFpDl+V81KVlZU4e/YsQkNDXXoHW1mEELhy5QpiY2Md3jmWPA+LyUudPXsW8fHxsmO43OnTp23eBps8G4vJS5nu1nv69GmEhYVJTuN8xcXFiI+Pd/pdiUkbWExeyjR9CwsL88piMvHGaSrx4jcRaRCLiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOXrZAcjzZWUBvr5V/5+XB/TtC6SkyM1Eno3FRPXyyivA9OnWj3/3HcuJ6o5TOaozo9F2KQFAQYF7s5B3YTGRSxw4IDsBeTIWE9VZSgowYoTtz5WUuDcLeRcWE9WZ0QisXWv7czNmVH2eqC5YTFRnubnVf87XFwgJcVsU8jIsJqqzxMTqP1dRAVy96rYo5GVYTFQnOTnAf/9r/5j33uN0jupGJ4QQskOQ8xUXFyM8PBxFRUUICwtz6tg5OUBaWs2Pd8WaJld+fyQfnzFRrb39du2OP3jQNTnIe7GYqFZycoD//Kd2XzNyJKd0VDssJnILrgSn2mAxUa0MHAh07So7BXk7FhPVSlZW3d5ucvq087OQ92IxUa2YtjeprcBA5+Yg78Ziolo5c6ZuXxcZ6dwc5N1YTFQrZ8/KTkAqYDFRrXTrJjsBqYDFRLVS18WS3J+JaoPFRDVmb8dKR7g/E9UGi4lqzN42J45wfyaqDRYT1Zi9bU4c4f5MVBssJnIL7s9EtcFiohqrz1QO4PvlqOZYTFRjdV1cSVRbLCaqMS6uJHdhMVGN1Xdx5f79zslB3o/FRDVW350on36aSwaoZlhMVGN79tTv6/V6LhmgmmExUY0YjcCnn9ZvjFGjgJgY5+Qh78ZiIrdZtoxTOaoZFhO5jY8Pp3JUMywmqpH6Lq4EgMpKrv6mmmExUY3k5TlnHK7+pppgMVGN1HWvb6K6YDERkeawmKhGKiqcMw5v40Q1wWKiGrlwwTnjXLrknHHIu7GYqEa+/NI540ycyLVM5BiLiRw6cQL4+GPnjOXnx7VM5BiLiRxKTQXuu885Y5WVcS0TOcZiIoeMRuA//3HeeFzLRI6wmMihnTtlJyDVsJjIoWbNZCcg1bCYyCFnvR3FhGuZyBEWE7ldYKDsBKR1LCZyyNlTOWfsVEDejcVEDq1e7dzxnD01JO/DYiK7jEZgyRLnjrlgAVd/k30sJnI73pSAHGExkduVl3P1N9nHYiK7XHWhmqu/yR4WE7mdjw+QlCQ7BWkZi4nscsVtvXlTAnKExUR2xca6ZlxO5cgeFhPZFRfnmnH5thSyh8VEdmVluWbckSO5lomqx2Iiu1w15eJaJrKHxUTVMhqBFStcMzbXMpE9LCaShhfAqTosJpKCa5nIHhYTVWv3bteNXVkJHD3quvHJs7GYSJrGjWUnIK1iMZE0wcGyE5BWsZioWvn5rh2fF7+pOiwmqlZ0tOwEpCoWE1XL1c+YTp507fjkuXRCCCE7BDlfcXExwsPDUVRUhLCwsFp/vdEItGzpgmC3qOtPX32/P9I2PmMiaSZNkp2AtIrFRDalpABTp7r2HG++yTfykm0sJrLJaARee8315+EbeckWFhPZ5K6bUvKNvGQLi4mk4lomsoXFRDbxbrkkE4uJbGrWTHYCUhmLiaTi3t9kC4uJbFq82D3nuXTJPechz8JiIitGI7BunXvONXEi1zKRNRYTSeXry7VMZI3FRFZSUoC0NPecq6KCa5nIGouJrBiNQHa2+87HtUx0KxYTWdm5U3YCUh2LiaxwDRPJxmIiIs1hMZEVvh2FZGMxkRVO5Ug2FhNZeecd957vwAH3no+0j8VEFrKygDVr3HvOkhL3no+0j8VEFpKT3X/OGTP4thSyxGIiCykpQEaGe8/Jt6XQrVhMZCErq+rDnfi2FLoVi4ks+PrKOS/flkI3YzERkeawmMgCF1eSFrCYyIKsxZXcYpduxmIiC5s3yzkvt9ilm7GYyMxdd9+15cIFOeclbWIxkZm77r5rS1GRvHOT9rCYyCwxUd65Fyzg6m+6gcVEmsDV33QzFhOZyZzKcfU33YzFRGay1zBx9TeZsJjIjBvEkVawmDSqvLwczz//PJo3b47AwEAkJSXh5ZdfRmVlpexoRC6nlx2AbJs3bx7eeust/Pvf/0bbtm1x4MABjB07FuHh4XjyySddck7ZU7kDB4BeveRmIG1gMWnU3r17cffdd2PIkCEAgMTERKxatQoHXLgPraydBUxCQ+Wen7SDUzmN6t27N7Zu3YrvvvsOAPDtt99i165dGDx4sMvOuXGjy4YmqhU+Y9Kov/71rygqKkJqaip8fX1RUVGBOXPmYOTIkTaPLy0tRWlpqfnXxcXFtTrfiRPu3+v7VuPGAY8+KjcDaQOLSaM+/PBDrFixAitXrkTbtm3xzTff4KmnnkJsbCxGjx5tdXxmZiZmzZpV5/OlpgJbtgC7dwMHDwLr19cnfe0MHgz07QsMHOi+c5K26YQQQnYIshYfH4+//e1vmDRpkvmx2bNnY8WKFThx4oTV8baeMcXHx6OoqAhhYWG1Pr/RaLmu6NNPgXnzaj2MlUWLgA4dbvw6KQmIian9OMXFxQgPD6/z90faxmdMGlVSUgIfH8tLgL6+vtUuFzAYDDAYDE47f0pK1ceNPPUvJh8f4J576lZEpBYWk0YNGzYMc+bMQbNmzdC2bVscOnQICxYswKMefBGmspJvO6GaYTFp1BtvvIEXXngBEydOREFBAWJjY/GXv/wFL774opQ8zlrjVFBg+UyMyBYWk0aFhoZi4cKFWLhwoewoAPh2FXIvrmOiGjl4UHYCUgmLiWqkpMQ547hw4Tp5ERYT1UhSknPG4dtOqCZYTOSQ0QiMHeucsYKCnDMOeTcWE7nV/731j8guFhM5lJICZGY6Zyw/P+eMQ96NxUQOGY3A9OnOGWvGDN4NhRxjMZFDzrxJgV7Pu6GQYywmcsiZ95sbNYrvlSPHWEzkVsuWcSpHjrGYyCFO5cjdWEzkVpzKUU2wmMitOJWjmmAxkVv5+nIqR46xmMitKiq4WRw5xmIih5x9I8yb9xInsoXFRA45c5M4Hx/n7VRA3ovFRG7Ffb+pJlhM5ND+/c4dj1M5coTFRA7FxspOQKphMZFDcXGyE5BqWEzkkLNfleO+3+QIi4kcunDBueNx329yhMVEDoWFOXe8Vq2cOx55HxYTOWQwOHc8TuXIERYTObRvn3PHc9Y96sh7sZjIrhMngHffde6Y3PebHGExkV2pqcC0ac4dk5vFkSMsJnLImTtYAsAzz3CzOLKPxUR2ZWUBa9Y4d8xXXuFUjuxjMZFdycnOH9PPj1M5so/FRHalpAAZGc4dc9gwTuXIPhYT2ZWVVfXhTGvXcipH9rGYyC5fX+ePyVflyBEWE7nd8OGcypF9LCayy9k7CwCcypFjLCayy5n7fZvwFk7kCIuJ3I63cCJHWExklyumcgD3/Sb7WExklyumcgBw+rRrxiXvwGIiKQIDZScgLWMxkV3OXlxpEhnpmnHJO7CYyC5XvFcO4FSO7GMxkV0pKa4Zl1M5sofFRHatXOmacZ29xxN5FxYTVctoBJYtc83YnMqRPSwmkkII2QlIy1hMVK2UFODRR10z9oIFfL8cVY/FRHb98otrxuXWJ2QPi4mqdeIE8PHHrhm7vJzvl6PqsZioWqmpwPjxrhuf75ej6rCYyK6ff3bNuL6+QFKSa8Ymz8diomq5Yr9vE259QvawmKhartjv+2acylF1WExUrcRE143t48OpHFWPxUTVOnXKdWNXVnIqR9VjMVG1OJUjWVhMJAWncmQPi4mq5ar9vgFO5cg+FhNVy1X7fZtwKkfVYTERkeawmKharpzKEdnDYqJquXoqR1QdFhNJc+CA7ASkVSwmqpar9vs2CQ117fjkuVhMZJMr9/s2adXKteOT52IxkTSbNslOQFrFYiKbXLnft0lFhWvHJ8/FYqJquWq/b5NXXuENCcg2FhPZ5Mr9vk38/HhDArKNxUQ2uXq/bwB4+WUgJsa15yDPxGKiarlqv2+T6dM5lSPbWExkkyv3+zbhVI6qw2Iim5KTXX+OsjJufUK2sZjIppQUICPD9efh1idkC4uJbHLHVA7gLpZkG4uJbHL1ft8mnMqRLSwmkopTObKFxURScSpHtrCYSCpO5cgWFhNJxakc2cJiIpu43zfJxGIim7jfN8nEYiIizWExkU3797vnPKdPu+c85FlYTGSTqzeJMwkMdM95yLOwmMgmd73r/6uv3HMe8iwsJrIpMdE95yksdM95yLOwmMim3bvdc56EBPechzwLi4msnDgBvPuue841YwZ3sSRrLCaykpoKTJvmnnPp9dzFkqyxmMim3Fz3nOeZZ3hDArLGYiIrWVnAmjXuORfvLUe2sJjIirs2iQN4QwKyjcVEVty1VAAAhg3jVI6ssZhIqrVrOZUjaywmsuKuC98m3JOJbsViIql8fLi9LlljMZFUlZXcXpessZhIOk7l6FYsJpKKUzmyhcVEVty53zencmQLi4msuHu/b07l6FYsJiLSHBYTWeGtm0g2FhNZced75YhsYTGRlfx8955v0yb3no+0j8VEVqKj3Xu+hg3dez7SPhYTSXf77bITkNawmMiC0QiMHevec7rr5prkOVhMJF14uOwEpDUsJrKQkgIsXerec7Zs6d7zkfaxmMjK+vXuPd+BA+49H2kfi4ksGI3AunXuPWdoqHvPR9rHYiIL7t69ksgWFhNZcOeNCEzGjXP/OUnbWEwk3UMPyU5AWsNiIgsypnIffMA7pZAlFhNZqKiQc17e9JJuxmIiCxkZQJ8+7j1nWhpvekmWWExkwWgEvvjCvefMzuZUjiyxmMjCzp1yzsvtdelmLCay4O79vgHeKYWssZjIgox1TLxTCt2KxUQWTp2Sc15O5ehmLCaywP2+SQtYTGRh82Y55z19Ws55SZtYTGRB1kLHwEA55yVtYjGRhV695Jw3MlLOeUmbWExkgTe7JC1gMZGF48dlJyBiMdFNjEbgtddkpyBiMRGRBrGYyCwlBcjMlHNuLhegm7GYyMLXX8s578iR3GGAbmAxkdmJE8CaNXLO7efHzeLoBhYTmaWmAvPmyTn3yy9zszi6gcWkYWfOnMGoUaPQqFEjBAUFoVOnTjh48KBLzynrzbTTp3MqRzfoZQcg2y5duoRevXphwIAB2LRpE5o0aYJTp04hIiLCZeeUvVyAUzkyYTFp1Lx58xAfH4+lS5eaH0t08WZJMm92OWYMp3J0A6dyGrV+/Xp07doV999/P5o0aYLOnTvjnXfecek5ZWwSZ7JsGadydAOLSaN++OEHLF68GCkpKdi8eTMee+wxPPHEE1i+fLnN40tLS1FcXGzx4Wk4lSMTTuU0qrKyEl27dsXcuXMBAJ07d8axY8ewePFiPPLII1bHZ2ZmYtasWfU6p8ypXO/enMrRDXzGpFExMTFo06aNxWOtW7dGXjVv/58+fTqKiorMH6frsJT6zJk6RXWKXbuAnBx55ydt4TMmjerVqxdOnjxp8dh3332HhIQEm8cbDAYYDIZ6nTMurl5fXm8lJXLPT9rBZ0wa9fTTT2Pfvn2YO3cuvv/+e6xcuRJvv/02Jk2a5LJzyrz4rdMBt90m7/ykLSwmjerWrRvWrVuHVatWoV27dvif//kfLFy4EA899JDLzrl7t8uGdkgI3sKJbuBUTsOGDh2KoUOHyo7hNgUFVTscEPEZE5nJvsZEZMJiIjPu901awWIis2bNZCcgqsJiIiLNYTGRWVaW7AREVVhMZJacLDsBURUWE5nxpXrSChYTme3fLzsBURUWE5mdOiU7AVEVFhMBqNqkbckSuRluec8yKYzFRJoxbpzsBKQVLCYCIPcuvCYu3DiBPAyLicx27JB7/jff5L7fVIXFRACqFlfKXmCp13Pfb6rCYiIA2lhced993PebqrCYCEDVNabf/15uhtWrOZWjKiwmAlBVCP/9r+wU8m5RTtrCYiIAcm/dZOLrCyQlyU5BWsBiIgBARYXsBFUZuO83ASwm+j8ZGUCPHrJTcCpHVVhMBKDqGtPevbJTEFVhMREAYOdO2QmIbmAxEQDu903awmIiAHLvwkt0KxYTAdDOXkxcLkAAi4n+j6+v7ARVuFyAABYTaQyXCxDAYqL/I3tnAaKbsZgIgDZ2FyAyYTERAN66ibSFxUQAgLw82QmIbmAxEQDgwgXZCaocOCA7AWkBi4kAAMeOyU5QpaREdgLSAhYTAdDOxe82bWQnIC1gMREAoFcv2QmqaGHDOpKPxUQAqvbb1gJO5QhgMRG0cXtwk4gI2QlIC1hMpCkTJ/JOKcRiIlQtrhw7VnaKG3jTS2IxEQCgrEx2gipPP82bXhKLiVA1dVqxQnaKKo0by05AWsBiIk3t9z1jBq8xEYuJAPTtKzvBDXo9rzERi4lQdfG7Z0/ZKaqMGsVrTMRiIlRtErdnj+wUVZYt41SOWEwE7ez3bcLtdYnFRJri48M7pRCLiTSmspJ3SiEWE2kQp3LEYiKcOSM7gaXTp2UnINlYTKQ5gYGyE5BsLCbSnMhI2QlINhYTIT9fdgJLnMoRi4kQHS07gSVO5YjFRIiLk53AEqdyxGIifPCB7AREllhMhNtuk53AEm96SSwmQkqK7ASWQkNlJyDZWEyEvDzZCSy1aiU7AcnGYiIcPy47AZElFpPijEbgtddkp7DEa0zEYlLcqVOyE1jjNSZiMSkuIwOYNk12CktBQbITkGwsJkJlpewElkaO5Pa6qmMxKU6L15j8/HinFNWxmBSXmys7gbWxY3mnFNWxmBSXmCg7gbW33wZycmSnIJlYTIpLSQEGDpSdwlpJiewEJBOLSXFGo/aenfj4aO/9e+ReLCbF7d4tO4E13imFWEykSbxTitpYTIrr1Ut2AiJrLCYi0hwWk+K0uI6JiMVERJrDYlKc1jaJIwJYTMpr1kx2AiJrLCYi0hwWk+I4lSMtYjEpjlM50iIWExFpDotJcZzKkRaxmBTHqRxpEYuJNGnTJtkJSCYWk+L275edwLaGDWUnIJlYTIr75hvZCWy7/XbZCUgmFpPCTpwA1qyRncK2jRtlJyCZWEwKS03V3s0uTQoLZScgmVhMitPijQgA7vmtOhaT4laulJ3ANt4mXG0sJoUZjcCyZbJT2MbbhKuNxaQwLe9eyduEq43FpDAt3oXXZOpU3iZcZSwmhaWkVBWAFr3yCqdyKmMxKe7bb2UnqB7vLacuFpPCsrKALVtkp7BNrweSkmSnIFlYTArz9ZWdoHrPPMNrTCpjMSlMyxe/eY1JbSwm0ixeY1IXi0lhWl7H5OPDa0wqYzEp7MwZ2QmqV1kJXL0qOwXJwmIizeJUTl0sJoX16iU7AZFtLCaFafkaE6mNxaQwLV9jIrWxmBQWFyc7AZFtLCaF8WaXpFUsJoVduCA7AZFtLCaFhYXJTkBkG4tJYQaD7AREtrGYFNasmewERLaxmIhIc1hMCsvKkp2AyDYWk8KSk2UnILKNxaSwlBTZCYhsYzEpTOsLLA8ckJ2AZGExKez4cdkJ7CspkZ2AZGExKcpoBF57TXYK+yIiZCcgWVhMiuKWJ6RlLCZFafkOKSYTJ/JOKapiMSkqJQV49FHZKezT64GQENkpSAYWk8IuXpSdwL7nnuNNL1XFYlKU0QisXy87hX18k7G6WEyK2rlTdgLHZszgNSZVsZgU1bev7ASO8RqTulhMikpJAaZOlZ3Cvvvu4zUmVbGYFOUJCyxXr+ZUTlUsJkV5ygJL3o1XTSwmItIcFpOitL6zAAD4+gJJSbJTkAwsJkV5wqtyFRXA1auyU5AMLCbSNF5jUhOLSVGesMCS1MViUhRv3URaxmJS1JkzshPUzOnTshOQDCwm0rTAQNkJSAYWk6Li4mQnqJnISNkJSAYWk6I8YR0TqYvF5CEyMzOh0+nw1FNPOWU8XvwmLWMxeYD9+/fj7bffRocOHWRHIXILFpPGXb16FQ899BDeeecdNGjQwGnj7t/vtKGInI7FpHGTJk3CkCFDcNddd9k9rrS0FMXFxRYf9vzyizNTEjkXi0nDVq9eja+//hqZmZkOj83MzER4eLj5Iz4+3u7xnvLmWN4mXE0sJo06ffo0nnzySaxYsQIBAQEOj58+fTqKiorMH6cdrEz0lI3+Q0NlJyAZ9LIDkG0HDx5EQUEBunTpYn6soqICO3fuxD//+U+UlpbC19fX/DmDwQBDLdrGU16Va9VKdgKSgcWkUXfeeSeOHDli8djYsWORmpqKv/71rxalRORtWEwaFRoainbt2lk8FhwcjEaNGlk9XhcrV9Z7CLc4cADo1Ut2CnI3XmNSkNEILFsmO0XNlJTITkAy8BmTB/n8889lR3C7Nm1kJyAZ+IxJQSkpwJw5slPUzK+/yk5AMrCYFBUbKztBzYwcyXvLqYjFpKj33pOdoGZ8fXmbcBWxmBRkNAK7dslOUTO//z1vE64iFpOCPOUuvADQv7/sBCQDi0lBFRWyE9TcxIm8xqQiFpOCMjKA8eNlp6g53ltOPSwmRVVWyk5QM3q95+yEQM7DYlKQ0QgsWSI7Rc088wwvfquIxaQgT7r4/corvMakIhYTaR6vMamHxaQgT3pVjtTEYlJQRgbQo4fsFDXj68uL3ypiMSnIaAT27pWdombuvJMXv1XEYlLQzp2yE9RcdjYvfquIxaQgT9nv24QXv9XDYlLQmTOyE9SOgxu+kBdiMZHmpabKTkDuxmJSUFyc7AS1ExwsOwG5G4tJQZ42leM1JvWwmIhIc1hMCvK0qRyph8WkoLw82QmI7GMxKYh3FyetYzEp6OhR2QmI7GMxKSgoSHYCIvtYTArq00d2AiL7WEykeQcOyE5A7sZiUpCnvSoXGio7Abkbi0lBFy7ITlA7rVrJTkDuxmJS0Jdfyk5AZB+LSTEnTgAffyw7Re3wGpN6WEyKSU0F5s2TnaJ2uB+TelhMCmrSRHaC2mnUSHYCcjcWk4I8beV3mzayE5C7sZgUdP687AS1U1goOwG5G4tJMUYjsGKF7BS1M26c7ATkbiwmxeTmyk5Qe6NGyU5A7sZiUkxiouwEtcdrTOphMSkmJQWYM0d2itqJiJCdgNyNxaQgHw/7U584kXfjVY2H/YiSM2zdKjtB7fFOKWphMSkmKwvYskV2itpLSpKdgNyJxaQYT9zve8QIICZGdgpyJxYTaV7XrrITkLuxmBRTUSE7Qe3NmMGL36phMSkmI8Mz9/zmxW+1sJgUYzQCX3whOwWRfSwmxXjiW1J8ffmqnGpYTKR5d9/NV+VUw2JSjKfdIQUA1q7lxW/VsJgU06yZ7AR1w4vfamExKcYTdxcg9bCYyCPw4rdaWEyK8cRX5dLSePFbNSwm0rzsbF78Vg2LSTGe+KocwIvfqmExKcZTX5UjtbCYiEhzWEyK8dSpHKmFxaQYTuXIE7CYiEhzWEyK4VSOPAGLSTGcypEnYDEp5swZ2QnqZuNG2QnInVhMisnPl52gbiorZScgd2IxKebiRdkJ6qZ7d9kJyJ1YTIoJCpKdoG4iI2UnIHdiMSnGE++QQuphMRGR5rCYFOOp65gOHJCdgNyJxaQYX1/ZCeqmpER2AnInFpNijh6VnaBu2rSRnYDcicWkmJMnZSeoG74qpxYWk0KMRmDDBtkp6ub0adkJyJ1YTOQRRo7kvt8qYTEpJCUFmDNHdoq68fMDQkJkpyB3YTEppls32Qnq5tFHeQsnlbCYFLN6tewEdRMcLDsBuROLSTGeeMNLAGjXTnYCcicWk0KMRmDbNtkp6sZTt2uhumExKWTnTtkJ6m7GDL4qpxIWk0I8eVtdvZ6vyqmExaSQxETZCepu3Di+KqcSFpNCUlKAp5+WnaJuEhJkJyB3YjEpxsdD/8QjImQnIHfy0B9Tqqtvv5WdoG4mTuTFb5WwmBSSlQVs2SI7Rd0VFMhOQO7CYlKIp24SR+phMZFH0OuBpCTZKchdWEwK8dT9vgEuF1ANi0khnrzA8vp12QnInVhMCvHkBZZLlwI5ObJTkLuwmBSSkuK5+zEBvFOKSlhMCjEagf37ZaeoO96QQB0sJoV48u4CpBYWk0I8+eK3jw+XC6iExUQe4a67uFxAJSwmhXjyOqbsbL5XTiUsJoV48lQO4HvlVMJiIiLNYTEpxJOncqQWFpNCPH0qR+pgMRGR5rCYFMKpHHkKFpNCOJUjT8FiUsiZM7IT1M/p07ITkLuwmMhjBAbKTkDuwmIij8HdBdTBYlJIfr7sBPWzcaPsBOQuLCaFREfLTlA/lZWyE5C7sJgUEhcnO0H9dO8uOwG5C4tJIZ6+jonXmNTBYlII1zGRp2AxEZHmsJgU4ulTOVIHi0khFy7ITlA/Bw7ITkDuwmJSyLFjshPUD9+Sog4WkyKMRmDFCtkp6qdNG9kJyF1YTIrIzZWdoP5atZKdgNyFxaSIxETZCYhqjsWkiJQUYM4c2SmIaobFpJBu3WQnqB/e4lwdLCaFePpGcTNm8KaXqmAxKeToUdkJ6sfXFwgJkZ2C3IHFpJDz52UnqJ/x44GYGNkpyB1YTIrwhnVMt90mOwG5C4tJEVzHRJ6ExUREmsNiUoQ37CzA98qpg8WkiL59ZSeov5EjuVxAFSwmRaSkAEuXyk5RfwUFshOQO7CYFHL2rOwE9ZeUJDsBuQOLSSG//CI7Qf2MGMF1TKpgMSnE0xdYNm8uOwG5C4tJEd6wwJIbxamDxaQIb1hgWVwsOwG5C4tJozIzM9GtWzeEhoaiSZMmuOeee3Dy5Mk6j+fpOwsAwNNPAzk5slOQO7CYNGrHjh2YNGkS9u3bh5ycHJSXlyMtLQ3Xrl2THU2qkhLZCcgd9LIDkG1ZWVkWv166dCmaNGmCgwcPom8dVkvGxTkrmVy8TbgaWEweoqioCADQsGFDm58vLS1FaWmp+dfFt1yQqahwXTYiZ+NUzgMIITBlyhT07t0b7dq1s3lMZmYmwsPDzR/x8fEWn8/IqNrPyJP5+nKBpSpYTB5g8uTJOHz4MFatWlXtMdOnT0dRUZH547SNd7yGh7sypetxozh1cCqncY8//jjWr1+PnTt3omnTptUeZzAYYDAY7I5Vjxf1NIHb6qqDxaRRQgg8/vjjWLduHT7//HM0r+eyZ6MR2LDBSeEkadRIdgJyF07lNGrSpElYsWIFVq5cidDQUOTn5yM/Px+//vprncbzhgWWvEuKOlhMGrV48WIUFRWhf//+iImJMX98+OGHdRrPGzaKA7jtiSo4ldMoIYRTx2vWzKnDEbkUnzERkeawmBThLVM5UgOLSRGcypEnYTEpwht2F9DpuPJbFSwmRXjDm3iFAK5elZ2C3IHFpIiBA4HOnWWnqD8uF1ADi0kRWVnAoUOyUxDVDItJEb6+shMQ1RyLiYg0h8WkCK5jIk/CYlIE1zGRJ2ExEZHmsJgU4S1TOS6wVAOLSRHeMpXjAks1sJjIo3CBpRpYTIrwlqkcqYHFpAgusCRPwmJSxBdfyE7gHDbuSkVeiMWkiFvuf+mxAgNlJyB3YDEpolcv2QmcIzJSdgJyBxaTInjxmzwJi0kR3rKOidTAYiKPcuCA7ATkDiwmRXjLVC40VHYCcgcWkyK8ZSrXqpXsBOQOLCYi0hwWkyK8ZSpHamAxKYJvSSFPwmJSRH6+7ARENcdiUsR338lOQFRzLCZFBAfLTkBUcywmRVy+LDuBc2zcKDsBuQOLSRHJybITOEdhoewE5A4sJkV4y+4CISGyE5A7sJgU4S3rmBo1kp2A3IHFpAhveUtKmzayE5A7sJjIo3CjODWwmBThLVM5UgOLSRHeMpUjNbCYiEhzWEyK4MJE8iQsJkVwdwHyJCwmRZSWyk7gHCdPyk5A7sBiUsSlS7ITOMe4cbITkDuwmBRx9arsBM4xapTsBOQOLCYF5OQAn3wiO4VzXLwoOwG5A4tJAWfOyE7gPEFBshOQO7CYFBAXJzuB8zRvLjsBuQOLSQEVFbITOA93F1ADi0kBGRnAtGmyUzgHdxdQA4tJEQMHyk7gHNxdQA0sJkVwdwHyJCwmRXB3AfIkLCbyKKdPy05A7sBiUoS37C4wciRgNMpOQa7GYlKEN+0uUFAgOwG5GotJEd6yuwCpgcWkCG/ZXYDUwGJShE4nOwFRzbGYFGA0AitWyE7hHD4+QFKS7BTkaiwmBeTmyk7gPOnpQEyM7BTkaiwm8ih8r5waWEwK4NtRyNOwmBTgTW9H4bYnamAxkUeJiJCdgNyBxaQAb5rKxcbKTkDuwGJSgDdN5bgfkxpYTORRuLuAGlhMCvCmPb+5u4AaWEwK8KY9vwHuLqACFpMivGXPb1IDi0kR+/fLTkBUcywmRfjwT5o8CH9cFfH557ITOA9fmfN+LCZFeNMzpsBA2QnI1bzox5XsSU6WncB5uMjS+7GYFOFNW+tyozjvx2JSRIMGshM4x4gR3ChOBSwmRZw4ITuBczRvLjsBuQOLSQFZWcCWLbJTENUci0kB3nSzS1IDi4mINIfFpABv2iiO1MBiUoA3bRTHPb/VwGIij8I9v9XAYlIAp3LkaVhMCvCmqdzEidzBUgUsJgWcOSM7gXMdPCg7Abkai4k8DncX8H4sJgXExclO4FzcXcD7sZgUwIvf5GlYTArwpovfpAYWExFpDotJAZzKkadhMSmAUznyNCwmItIcFpMCOJUjT8NiUkDfvrITENUOi0kBKSnAqFGyUxDVHItJEULITkBUcywmRXjL7ZtIDSwmRRw/LjuB85w+LTsBuRqLSRH+/rITOE9qquwE5GosJkUkJ8tO4ByjRgGdOslOQa7GYlJEbq7sBM4RFSU7AbkDi0kBRiOwYYPsFEQ1x2JSgLc8WyJ1sJgUwLekkKdhMSmAuwuQp2ExEZHmsJgUwKkceRoWkwK8aSr300+yE5A7sJgU4E03vPzPf4CcHNkpyNVYTORxSkpkJyBXYzEpgDe8JE/DYlIAL36Tp2ExKcCbLn6TGlhMRKQ5LCYFcCpHnobFpABO5cjTsJiISHNYTArgVI48DYtJAZzKkadhMRGR5rCYFMCpHHkaFpMCOJUjT8NiUoA37S4A8IaXKmAxkccJDJSdgFyNxaQA7i5AnobFpABe/CZPw2JSAC9+k6dhMRGR5rCYFMCpHHkaFpMCOJUjT8NiIiLNYTEpgFM58jQsJgVwKkeehsWkcYsWLULz5s0REBCALl264IsvvpAdicjlWEwa9uGHH+Kpp57Cc889h0OHDqFPnz4YNGgQ8mo5N+NUjjwNi0nDFixYgHHjxmH8+PFo3bo1Fi5ciPj4eCxevLhW43AqR56GxaRR169fx8GDB5GWlmbxeFpaGvbs2WN1fGlpKYqLiy0+TLi7AHkaFpNGXbx4ERUVFYiKirJ4PCoqCvn5+VbHZ2ZmIjw83PwRHx9v/py3vYk3NVV2AnI1vewAZJ9Op7P4tRDC6jEAmD59OqZMmWL+dXFxsbmcBg4EsrOB77+v3bmPHwdCQoCbOq5enDFejx5Ap07ATU8IyQuxmDQqMjISvr6+Vs+OCgoKrJ5FAYDBYIDBYKh2vIEDqz6IPAGnchrl7++PLl26ICcnx+LxnJwc9OzZU1IqIvfgMyYNmzJlCh5++GF07doVPXr0wNtvv428vDw89thjsqMRuRSLScMeeOABFBYW4uWXX8a5c+fQrl07bNy4EQkJCbKjEbmUTgghZIcg5ysuLkZ4eDiKiooQFhYmO47Tefv3pzpeYyIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHNYTESkOSwmItIcFhMRaQ6LiYg0h8VERJrDYiIizWExEZHmsJiISHP0sgOQawghAADFxcWSk7iG6fsyfZ/kXVhMXurKlSsAgPj4eMlJXOvKlSsIDw+XHYOcTCf4T45XqqysxNmzZxEaGgqdTmfzmOLiYsTHx+P06dMICwur1/mcNVZNxxFC4MqVK4iNjYWPD69IeBs+Y/JSPj4+aNq0aY2ODQsLq3cxOXusmozDZ0rei//UEJHmsJiISHNYTAozGAx46aWXYDAYNDOWMzOR5+LFbyLSHD5jIiLNYTERkeawmIhIc1hMRKQ5LCaNyszMRLdu3RAaGoomTZrgnnvuwcmTJy2OGTNmDHQ6ncVH9+7dLY4pLS3FbbfdZnVcdHS0+ZhffvkFHTp0gI+PD3Q6HaKiorB3716LcYxGI5o3b241julj0qRJ2L17NwIDA+1mysvLw+DBg+Hn5wcfHx/4+flh6NCh+Pnnn83HHDlyBD179oSvry98fHwQEBCAhx9+GJcvX7bIlJeXh2HDhiE4OBiRkZF44okncP36dYtjjhw5gn79+iEwMBBxcXF4+eWX+f46TyBIk9LT08XSpUvF0aNHxTfffCOGDBkimjVrJq5evWo+ZvTo0SIjI0OcO3fO/FFYWGgxzmOPPSZCQ0NFYmKiyM7OFj179hRt2rQR586dMx/TsmVL4ePjI+bOnSvef/99ER4eLgwGgyguLhZCCFFeXi4aNmwoDAaDmD9/vnjllVeEn5+faNiwocjKyhIAxIYNG4S/v7/w8/MTHTp0ENOnTxc6nU6kpKSYM5WXl4t27dqJ2NhY0aRJEzFv3jzRuHFjERcXJzp27CjKy8tFUVGRiIqKEjExMaJFixZi7ty5IigoSERHR4uhQ4eaM5vGGjBggPj6669FTk6OiI2NFZMnTzYfYxrrwQcfFEeOHBH//e9/RWhoqHj11Vdd8mdGzsNi8hAFBQUCgNixY4f5sdGjR4u777672q+5fPmy8PPzE7///e9Fx44dhRBCnDlzRvj4+IisrCwhhBDHjh0TAMTEiRPNX7djxw4BQMycOVMIIcSaNWsEALFo0SLzMW+++aYAIIYMGSKSk5PF888/LwCIjIwMc6YJEyYIAOL48eNCCCE2btwodDqd8PPzE6tXrxZCCLFq1Srh7+9vzrRo0SIREhIiAIh9+/YJIYTIzMwUkZGRAoA4ceKEeSwfHx9x5swZc6ZVq1YJg8EgioqKhBBCLFq0SISHh4vffvvNfExmZqaIjY0VlZWVNf/NJ7fjVM5DFBUVAQAaNmxo8fjnn3+OJk2aoGXLlvjTn/6EgoIC8+cOHjyIsrIyJCcnw2g0IjY2Fr169UJoaCg+++wzAMCnn34KABg/frz56/r27Qu9Xo/s7GwAwCeffAIAePDBB83HjBw5EgCwZcsWPProo9i6dSv0ej2ioqLMmUznWLlyJQBg7969aN68OcrKypCWlgYASE9Px/Xr15GYmIg9e/Zg7969aNGiBcLDw3H77bebj7l48SJCQ0OxZ88e81jt2rVDbGysOVN6ejpKS0tx8OBB8zH9+vWzWKyZnp6Os2fPIjc3t5Z/AuROLCYPIITAlClT0Lt3b7Rr1878+KBBg/DBBx9g27ZteO2117B//37ccccdKC0tBQDk5+fD398f/fv3x/Lly7F582a888478PHxwbvvvovCwkL88MMPAICoqCiLcwYFBeHChQsAgDNnzkCn06FBgwbmz5v+//r16xgzZgwuXLiAoKAgi0yzZs0CAPzzn/9EaWkp8vPzERQUBH9/f/PXN2jQAP7+/ggODkZ+fj7y8/Oh1+vRpEkT87lM2SIiIpCfn2/+3m7NbBrL3jGmX5uOIW3i7gIeYPLkyTh8+DB27dpl8fgDDzxg/v927dqha9euSEhIwGeffYYRI0aYPzdo0CDz/7dv3x6dOnXCvn378O9//9v8eHVbozjSqFEj87MWnU5nkclUPkVFReZnT7bOI6ouKZg/Z7pofvPnTf+9+fHqxrJ3jGmsun6/5B58xqRxjz/+ONavX4/t27c73MYkJiYGCQkJMBqNAIDo6Ghcv34dly5dsjiusLAQ0dHRMBqNSEpKAmD9DKKkpASNGzcGAMTFxUEIYTHO4cOHAcD8DK5x48a4du2axRimaWWDBg1gNBoRHR2Na9euWWS6dOkSysrKUFJSgqioKERHR6OsrAznz5+3GqeoqMj8jCc6Otoqs2kse8eYxrr1mRRpC4tJo4QQmDx5MtauXYtt27ahefPmDr+msLAQp0+fRkxMDACgS5cu8PPzQ05OjvmYc+fO4ciRI7hy5QpiYmIwbNgwAMCSJUvMx+zatQvl5eXm60B33303AOCjjz4yHzNjxgwAwMSJEwEAd911F8rLy7F7927zMe+++y4AmM/Vo0cP/PjjjxaZsrOz4e/vj9zcXPTs2RM9evTAqVOnUFRUhK+++sp8TGRkJK5cuYKePXsCAHr06IGjR4/i3Llz5vNlZ2fDYDCgS5cu5mN27txpsYQgOzsbsbGxSExMdPj7SRJJu+xOdk2YMEGEh4eLzz//3GI5QElJiRBCiCtXroipU6eKPXv2iB9//FFs375d9OjRQ8TFxZlf5heiarlASEiIeO2118Snn34qunTpIsLCwkRoaKjIzc0VQtxYLpCZmSlWrFhhd7nA4sWLza/IBQYGivLyciHEjVcADQaDeOmll8Szzz4rdDqdCAwMNGcyvcQfFxcnoqKixPz580WTJk1E06ZNzcsFLl++bF4u0LJlS5GZmSmCgoJETEyMzeUCd955p/j666/Fli1bRNOmTS2WC5jGGjlypDhy5IhYu3atCAsL43IBD8Bi0igANj+WLl0qhBCipKREpKWlicaNGws/Pz/RrFkzMXr0aJGXl2cxzq+//ipSUlKETqcTAITBYBCDBg0Sx44dMx9z8eJF0b59e/MxTZo0Ebt377YY5+TJkyIhIcEiS3Z2tsUxW7ZsEX5+fubP+/n5iZEjR1pk+umnn0RGRobw9fUVOp1O6PV6MXjwYItjDh8+LLp37y58fHyETqcT/v7+4qGHHhKXLl2yON9PP/0khgwZIgIDA0XDhg3F5MmTLZYGmMbq06ePMBgMIjo6WsycOZNLBTwAtz0hIs3hNSYi0hwWExFpDouJiDSHxUREmsNiIiLNYTERkeawmIhIc1hMRKQ5LCYi0hwWExFpDouJiDSHxUREmvP/AX/H5lDXvM9/AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# ## Visualize Product Outlines\n", + "# Plot the generated product outlines to visualize.\n", + "product_outlines_gdf = gpd.read_file(product_outlines_file)\n", + "fig, ax = plt.subplots(figsize=(4, 10))\n", + "product_outlines_gdf.plot(ax=ax, edgecolor=\"blue\", facecolor=\"none\")\n", + "plt.title(f\"Product Outlines for {utm_zone}\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Global product outlines saved to test_output\\global_product_outlines_utm_zones\n" + ] + } + ], + "source": [ + "# ## Generate Global Product Outlines Grouped by UTM Zones\n", + "# Create global product outlines split by UTM zones.\n", + "global_output_folder = os.path.join(output_folder, \"global_product_outlines_utm_zones\")\n", + "grid.generate_global_product_outlines_by_utm(\n", + " output_folder=global_output_folder,\n", + " driver=\"ESRI Shapefile\",\n", + " get_footprints=False,\n", + " naming_convention=\"zone\"\n", + ")\n", + "\n", + "print(f\"Global product outlines saved to {global_output_folder}\")\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "processing file: raster_outlines_UTM_zone_10N.shp\n", + "processing file: raster_outlines_UTM_zone_10S.shp\n", + "processing file: raster_outlines_UTM_zone_11N.shp\n", + "processing file: raster_outlines_UTM_zone_11S.shp\n", + "processing file: raster_outlines_UTM_zone_12N.shp\n", + "processing file: raster_outlines_UTM_zone_12S.shp\n", + "processing file: raster_outlines_UTM_zone_13N.shp\n", + "processing file: raster_outlines_UTM_zone_13S.shp\n", + "processing file: raster_outlines_UTM_zone_14N.shp\n", + "processing file: raster_outlines_UTM_zone_14S.shp\n", + "processing file: raster_outlines_UTM_zone_15N.shp\n", + "processing file: raster_outlines_UTM_zone_15S.shp\n", + "processing file: raster_outlines_UTM_zone_16N.shp\n", + "processing file: raster_outlines_UTM_zone_16S.shp\n", + "processing file: raster_outlines_UTM_zone_17N.shp\n", + "processing file: raster_outlines_UTM_zone_17S.shp\n", + "processing file: raster_outlines_UTM_zone_18N.shp\n", + "processing file: raster_outlines_UTM_zone_18S.shp\n", + "processing file: raster_outlines_UTM_zone_19N.shp\n", + "processing file: raster_outlines_UTM_zone_19S.shp\n", + "processing file: raster_outlines_UTM_zone_1N.shp\n", + "processing file: raster_outlines_UTM_zone_1S.shp\n", + "processing file: raster_outlines_UTM_zone_20N.shp\n", + "processing file: raster_outlines_UTM_zone_20S.shp\n", + "processing file: raster_outlines_UTM_zone_21N.shp\n", + "processing file: raster_outlines_UTM_zone_21S.shp\n", + "processing file: raster_outlines_UTM_zone_22N.shp\n", + "processing file: raster_outlines_UTM_zone_22S.shp\n", + "processing file: raster_outlines_UTM_zone_23N.shp\n", + "processing file: raster_outlines_UTM_zone_23S.shp\n", + "processing file: raster_outlines_UTM_zone_24N.shp\n", + "processing file: raster_outlines_UTM_zone_24S.shp\n", + "processing file: raster_outlines_UTM_zone_25N.shp\n", + "processing file: raster_outlines_UTM_zone_25S.shp\n", + "processing file: raster_outlines_UTM_zone_26N.shp\n", + "processing file: raster_outlines_UTM_zone_26S.shp\n", + "processing file: raster_outlines_UTM_zone_27N.shp\n", + "processing file: raster_outlines_UTM_zone_27S.shp\n", + "processing file: raster_outlines_UTM_zone_28N.shp\n", + "processing file: raster_outlines_UTM_zone_28S.shp\n", + "processing file: raster_outlines_UTM_zone_29N.shp\n", + "processing file: raster_outlines_UTM_zone_29S.shp\n", + "processing file: raster_outlines_UTM_zone_2N.shp\n", + "processing file: raster_outlines_UTM_zone_2S.shp\n", + "processing file: raster_outlines_UTM_zone_30N.shp\n", + "processing file: raster_outlines_UTM_zone_30S.shp\n", + "processing file: raster_outlines_UTM_zone_31N.shp\n", + "processing file: raster_outlines_UTM_zone_31S.shp\n", + "processing file: raster_outlines_UTM_zone_32N.shp\n", + "processing file: raster_outlines_UTM_zone_32S.shp\n", + "processing file: raster_outlines_UTM_zone_33N.shp\n", + "processing file: raster_outlines_UTM_zone_33S.shp\n", + "processing file: raster_outlines_UTM_zone_34N.shp\n", + "processing file: raster_outlines_UTM_zone_34S.shp\n", + "processing file: raster_outlines_UTM_zone_35N.shp\n", + "processing file: raster_outlines_UTM_zone_35S.shp\n", + "processing file: raster_outlines_UTM_zone_36N.shp\n", + "processing file: raster_outlines_UTM_zone_36S.shp\n", + "processing file: raster_outlines_UTM_zone_37N.shp\n", + "processing file: raster_outlines_UTM_zone_37S.shp\n", + "processing file: raster_outlines_UTM_zone_38N.shp\n", + "processing file: raster_outlines_UTM_zone_38S.shp\n", + "processing file: raster_outlines_UTM_zone_39N.shp\n", + "processing file: raster_outlines_UTM_zone_39S.shp\n", + "processing file: raster_outlines_UTM_zone_3N.shp\n", + "processing file: raster_outlines_UTM_zone_3S.shp\n", + "processing file: raster_outlines_UTM_zone_40N.shp\n", + "processing file: raster_outlines_UTM_zone_40S.shp\n", + "processing file: raster_outlines_UTM_zone_41N.shp\n", + "processing file: raster_outlines_UTM_zone_41S.shp\n", + "processing file: raster_outlines_UTM_zone_42N.shp\n", + "processing file: raster_outlines_UTM_zone_42S.shp\n", + "processing file: raster_outlines_UTM_zone_43N.shp\n", + "processing file: raster_outlines_UTM_zone_43S.shp\n", + "processing file: raster_outlines_UTM_zone_44N.shp\n", + "processing file: raster_outlines_UTM_zone_44S.shp\n", + "processing file: raster_outlines_UTM_zone_45N.shp\n", + "processing file: raster_outlines_UTM_zone_45S.shp\n", + "processing file: raster_outlines_UTM_zone_46N.shp\n", + "processing file: raster_outlines_UTM_zone_46S.shp\n", + "processing file: raster_outlines_UTM_zone_47N.shp\n", + "processing file: raster_outlines_UTM_zone_47S.shp\n", + "processing file: raster_outlines_UTM_zone_48N.shp\n", + "processing file: raster_outlines_UTM_zone_48S.shp\n", + "processing file: raster_outlines_UTM_zone_49N.shp\n", + "processing file: raster_outlines_UTM_zone_49S.shp\n", + "processing file: raster_outlines_UTM_zone_4N.shp\n", + "processing file: raster_outlines_UTM_zone_4S.shp\n", + "processing file: raster_outlines_UTM_zone_50N.shp\n", + "processing file: raster_outlines_UTM_zone_50S.shp\n", + "processing file: raster_outlines_UTM_zone_51N.shp\n", + "processing file: raster_outlines_UTM_zone_51S.shp\n", + "processing file: raster_outlines_UTM_zone_52N.shp\n", + "processing file: raster_outlines_UTM_zone_52S.shp\n", + "processing file: raster_outlines_UTM_zone_53N.shp\n", + "processing file: raster_outlines_UTM_zone_53S.shp\n", + "processing file: raster_outlines_UTM_zone_54N.shp\n", + "processing file: raster_outlines_UTM_zone_54S.shp\n", + "processing file: raster_outlines_UTM_zone_55N.shp\n", + "processing file: raster_outlines_UTM_zone_55S.shp\n", + "processing file: raster_outlines_UTM_zone_56N.shp\n", + "processing file: raster_outlines_UTM_zone_56S.shp\n", + "processing file: raster_outlines_UTM_zone_57N.shp\n", + "processing file: raster_outlines_UTM_zone_57S.shp\n", + "processing file: raster_outlines_UTM_zone_58N.shp\n", + "processing file: raster_outlines_UTM_zone_58S.shp\n", + "processing file: raster_outlines_UTM_zone_59N.shp\n", + "processing file: raster_outlines_UTM_zone_59S.shp\n", + "processing file: raster_outlines_UTM_zone_5N.shp\n", + "processing file: raster_outlines_UTM_zone_5S.shp\n", + "processing file: raster_outlines_UTM_zone_60N.shp\n", + "processing file: raster_outlines_UTM_zone_60S.shp\n", + "processing file: raster_outlines_UTM_zone_6N.shp\n", + "processing file: raster_outlines_UTM_zone_6S.shp\n", + "processing file: raster_outlines_UTM_zone_7N.shp\n", + "processing file: raster_outlines_UTM_zone_7S.shp\n", + "processing file: raster_outlines_UTM_zone_8N.shp\n", + "processing file: raster_outlines_UTM_zone_8S.shp\n", + "processing file: raster_outlines_UTM_zone_9N.shp\n", + "processing file: raster_outlines_UTM_zone_9S.shp\n", + "test_output\\merged_global_outlines.shp\n", + "test_output\\merged_global_outlines.shp\n" + ] + }, + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[14], line 4\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m# ## Merge UTM Zone Files into WGS84\u001b[39;00m\n\u001b[0;32m 2\u001b[0m \u001b[38;5;66;03m# Merge all the generated UTM zone files into a single WGS84 file.\u001b[39;00m\n\u001b[0;32m 3\u001b[0m merged_file \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(output_folder, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmerged_global_outlines.shp\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m----> 4\u001b[0m \u001b[43mmerge_utm_files_to_wgs84\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m 5\u001b[0m \u001b[43m \u001b[49m\u001b[43mutm_folder\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mglobal_output_folder\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 6\u001b[0m \u001b[43m \u001b[49m\u001b[43moutput_file\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmerged_file\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 7\u001b[0m \u001b[43m \u001b[49m\u001b[43mdriver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mESRI Shapefile\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\n\u001b[0;32m 8\u001b[0m \u001b[43m)\u001b[49m\n\u001b[0;32m 10\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mMerged global outlines saved to \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmerged_file\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", + "File \u001b[1;32md:\\Major-TOM-main\\MajorTOM\\grid.py:385\u001b[0m, in \u001b[0;36mmerge_utm_files_to_wgs84\u001b[1;34m(utm_folder, output_file, ext, driver)\u001b[0m\n\u001b[0;32m 383\u001b[0m export_columns \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrid_cell\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mutm_crs\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[0;32m 384\u001b[0m \u001b[38;5;66;03m# Drop temporary columns and return the GeoDataFrame\u001b[39;00m\n\u001b[1;32m--> 385\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m reprojected_gdf[export_columns]\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\geopandas\\geodataframe.py:1264\u001b[0m, in \u001b[0;36mGeoDataFrame.to_file\u001b[1;34m(self, filename, driver, schema, index, **kwargs)\u001b[0m\n\u001b[0;32m 1173\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Write the ``GeoDataFrame`` to a file.\u001b[39;00m\n\u001b[0;32m 1174\u001b[0m \n\u001b[0;32m 1175\u001b[0m \u001b[38;5;124;03mBy default, an ESRI shapefile is written, but any OGR data source\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 1260\u001b[0m \n\u001b[0;32m 1261\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 1262\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mgeopandas\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mfile\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m _to_file\n\u001b[1;32m-> 1264\u001b[0m \u001b[43m_to_file\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdriver\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mschema\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mindex\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\geopandas\\io\\file.py:612\u001b[0m, in \u001b[0;36m_to_file\u001b[1;34m(df, filename, driver, schema, index, mode, crs, engine, **kwargs)\u001b[0m\n\u001b[0;32m 609\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmode\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m should be one of \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mw\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m or \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124ma\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m, got \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmode\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m instead\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 611\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfiona\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m--> 612\u001b[0m \u001b[43m_to_file_fiona\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdriver\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mschema\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcrs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 613\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpyogrio\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m 614\u001b[0m _to_file_pyogrio(df, filename, driver, schema, crs, mode, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\geopandas\\io\\file.py:641\u001b[0m, in \u001b[0;36m_to_file_fiona\u001b[1;34m(df, filename, driver, schema, crs, mode, **kwargs)\u001b[0m\n\u001b[0;32m 637\u001b[0m crs_wkt \u001b[38;5;241m=\u001b[39m crs\u001b[38;5;241m.\u001b[39mto_wkt(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mWKT1_GDAL\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 638\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m fiona\u001b[38;5;241m.\u001b[39mopen(\n\u001b[0;32m 639\u001b[0m filename, mode\u001b[38;5;241m=\u001b[39mmode, driver\u001b[38;5;241m=\u001b[39mdriver, crs_wkt\u001b[38;5;241m=\u001b[39mcrs_wkt, schema\u001b[38;5;241m=\u001b[39mschema, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs\n\u001b[0;32m 640\u001b[0m ) \u001b[38;5;28;01mas\u001b[39;00m colxn:\n\u001b[1;32m--> 641\u001b[0m \u001b[43mcolxn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwriterecords\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdf\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43miterfeatures\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\fiona\\collection.py:558\u001b[0m, in \u001b[0;36mCollection.writerecords\u001b[1;34m(self, records)\u001b[0m\n\u001b[0;32m 556\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmode \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m (\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124ma\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mw\u001b[39m\u001b[38;5;124m\"\u001b[39m):\n\u001b[0;32m 557\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mOSError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcollection not open for writing\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m--> 558\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msession\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwriterecs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrecords\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m 559\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_len \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msession\u001b[38;5;241m.\u001b[39mget_length()\n\u001b[0;32m 560\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_bounds \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[1;32mfiona\\\\ogrext.pyx:1392\u001b[0m, in \u001b[0;36mfiona.ogrext.WritingSession.writerecs\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\geopandas\\geodataframe.py:951\u001b[0m, in \u001b[0;36mGeoDataFrame.iterfeatures\u001b[1;34m(self, na, show_bbox, drop_id)\u001b[0m\n\u001b[0;32m 949\u001b[0m feature[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mtype\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFeature\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 950\u001b[0m feature[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mproperties\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m properties_items\n\u001b[1;32m--> 951\u001b[0m feature[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgeometry\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m mapping(geom) \u001b[38;5;28;01mif\u001b[39;00m geom \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m 953\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m show_bbox:\n\u001b[0;32m 954\u001b[0m feature[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbbox\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m geom\u001b[38;5;241m.\u001b[39mbounds \u001b[38;5;28;01mif\u001b[39;00m geom \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\shapely\\geometry\\base.py:113\u001b[0m, in \u001b[0;36mBaseGeometry.__bool__\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 112\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__bool__\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m--> 113\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mis_empty\u001b[49m \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mFalse\u001b[39;00m\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\shapely\\geometry\\base.py:629\u001b[0m, in \u001b[0;36mBaseGeometry.is_empty\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 626\u001b[0m \u001b[38;5;129m@property\u001b[39m\n\u001b[0;32m 627\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mis_empty\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m 628\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"True if the set of points in this geometry is empty, else False\"\"\"\u001b[39;00m\n\u001b[1;32m--> 629\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mbool\u001b[39m(\u001b[43mshapely\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mis_empty\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m)\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\shapely\\decorators.py:77\u001b[0m, in \u001b[0;36mmultithreading_enabled..wrapped\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m arr \u001b[38;5;129;01min\u001b[39;00m array_args:\n\u001b[0;32m 76\u001b[0m arr\u001b[38;5;241m.\u001b[39mflags\u001b[38;5;241m.\u001b[39mwriteable \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[1;32m---> 77\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 78\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[0;32m 79\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m arr, old_flag \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(array_args, old_flags):\n", + "File \u001b[1;32mc:\\Users\\stays\\anaconda3\\envs\\PytKerasNewHQ\\Lib\\site-packages\\shapely\\predicates.py:162\u001b[0m, in \u001b[0;36mis_empty\u001b[1;34m(geometry, **kwargs)\u001b[0m\n\u001b[0;32m 137\u001b[0m \u001b[38;5;129m@multithreading_enabled\u001b[39m\n\u001b[0;32m 138\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mis_empty\u001b[39m(geometry, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m 139\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns True if a geometry is an empty point, polygon, etc.\u001b[39;00m\n\u001b[0;32m 140\u001b[0m \n\u001b[0;32m 141\u001b[0m \u001b[38;5;124;03m Parameters\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 160\u001b[0m \u001b[38;5;124;03m False\u001b[39;00m\n\u001b[0;32m 161\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[1;32m--> 162\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mlib\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mis_empty\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgeometry\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[1;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "\n", + "# ## Merge UTM Zone Files into WGS84\n", + "# Merge all the generated UTM zone files into a single WGS84 file.\n", + "merged_file = os.path.join(output_folder, \"merged_global_outlines.shp\")\n", + "merge_utm_files_to_wgs84(\n", + " utm_folder=global_output_folder,\n", + " output_file=merged_file,\n", + " driver=\"ESRI Shapefile\"\n", + ")\n", + "\n", + "print(f\"Merged global outlines saved to {merged_file}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ## Visualize Merged Global Outlines\n", + "# Load and plot the merged global outlines.\n", + "merged_outlines_gdf = gpd.read_file(merged_file)\n", + "\n", + "merged_outlines_gdf.head()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save grid points to a GeoJSON file.\n", + "grid_output_file = os.path.join(output_folder, f\"grid_{dist}km.shp\")\n", + "grid_points_export = grid_points[[\"grid_cell\", \"geometry\"]].copy()\n", + "grid_points_export.to_file(grid_output_file, driver=\"ESRI Shapefile\")\n", + "print(f\"Grid points saved to {grid_output_file}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ## Clean-Up Outputs\n", + "# Optionally, clean up generated files.\n", + "clean_up = False\n", + "if clean_up:\n", + " import shutil\n", + " shutil.rmtree(output_folder)\n", + " print(\"Cleaned up all generated files.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Now, we can similarly work from the metadata parquet files" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# provide a folder path to the output folder\n", + "# ### Create Output Directory\n", + "# Set up directories to save output files.\n", + "output_folder = \"test_filtered_output\"\n", + "os.makedirs(output_folder, exist_ok=True)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Grid initialized with resolution: 10 km\n" + ] + } + ], + "source": [ + "# ### Initialize Grid\n", + "# Create a `Grid` object with a 10 km grid resolution.\n", + "dist = 10 # Grid resolution in kilometers\n", + "grid = Grid(dist)\n", + "print(\"Grid initialized with resolution:\", dist, \"km\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Metadata loaded successfully!\n" + ] + } + ], + "source": [ + "# Get the S2L2A metadata file\n", + "\n", + "metadata_file = \"metadata_s2l2.parquet\"\n", + "metadata = pd.read_parquet(metadata_file)\n", + "print(\"Metadata loaded successfully!\")\n", + "metadata['timestamp'] = pd.to_datetime(metadata.timestamp)\n", + "gdf = gpd.GeoDataFrame(\n", + " metadata, geometry=gpd.points_from_xy(metadata.centre_lon, metadata.centre_lat), crs=metadata.crs.iloc[0]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "def filter_metadata(df,\n", + " region=None,\n", + " daterange=None,\n", + " cloud_cover=(0,100),\n", + " nodata=(0, 1.0)\n", + " ):\n", + " \"\"\"Filters the Major-TOM dataframe based on several parameters\n", + "\n", + " Args:\n", + " df (geopandas dataframe): Parent dataframe\n", + " region (shapely geometry object) : Region of interest\n", + " daterange (tuple) : Inclusive range of dates (example format: '2020-01-01')\n", + " cloud_cover (tuple) : Inclusive percentage range (0-100) of cloud cover\n", + " nodata (tuple) : Inclusive fraction (0.0-1.0) of no data allowed in a sample\n", + "\n", + " Returns:\n", + " df: a filtered dataframe\n", + " \"\"\"\n", + " # temporal filtering\n", + " if daterange is not None:\n", + " assert (isinstance(daterange, list) or isinstance(daterange, tuple)) and len(daterange)==2\n", + " df = df[df.timestamp >= daterange[0]]\n", + " df = df[df.timestamp <= daterange[1]]\n", + " \n", + " # spatial filtering\n", + " if region is not None:\n", + " idxs = df.sindex.query(region)\n", + " df = df.take(idxs)\n", + " # cloud filtering\n", + " if cloud_cover is not None:\n", + " df = df[df.cloud_cover >= cloud_cover[0]]\n", + " df = df[df.cloud_cover <= cloud_cover[1]]\n", + "\n", + " # spatial filtering\n", + " if nodata is not None:\n", + " df = df[df.nodata >= nodata[0]]\n", + " df = df[df.nodata <= nodata[1]]\n", + "\n", + " return df\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
grid_cellgrid_row_ugrid_col_rproduct_idtimestampcloud_covernodatacentre_latcentre_loncrsparquet_urlparquet_rowgeometry
1591452453U_119R453119S2A_MSIL2A_20210510T100031_N0500_R122_T33TVF_2...2021-05-10 10:00:310.0000000.040.73410914.155341EPSG:32633https://huggingface.co/datasets/Major-TOM/Core...33POINT (14.155 40.734)
1593374454U_120R454120S2B_MSIL2A_20220719T095559_N0400_R122_T33TVF_2...2022-07-19 09:55:590.0000000.040.82386114.292709EPSG:32633https://huggingface.co/datasets/Major-TOM/Core...455POINT (14.293 40.824)
1595266455U_120R455120S2B_MSIL2A_20220719T095559_N0400_R122_T33TVF_2...2022-07-19 09:55:590.0000000.040.91367114.311585EPSG:32633https://huggingface.co/datasets/Major-TOM/Core...347POINT (14.312 40.914)
1595265455U_119R455119S2A_MSIL2A_20200113T095351_N0500_R079_T33TVF_2...2020-01-13 09:53:510.0000000.040.91373114.192730EPSG:32633https://huggingface.co/datasets/Major-TOM/Core...346POINT (14.193 40.914)
1595264455U_118R455118S2A_MSIL2A_20230318T095031_N0509_R079_T33TVF_2...2023-03-18 09:50:312.4085060.040.91379114.073875EPSG:32633https://huggingface.co/datasets/Major-TOM/Core...345POINT (14.074 40.914)
\n", + "
" + ], + "text/plain": [ + " grid_cell grid_row_u grid_col_r \\\n", + "1591452 453U_119R 453 119 \n", + "1593374 454U_120R 454 120 \n", + "1595266 455U_120R 455 120 \n", + "1595265 455U_119R 455 119 \n", + "1595264 455U_118R 455 118 \n", + "\n", + " product_id \\\n", + "1591452 S2A_MSIL2A_20210510T100031_N0500_R122_T33TVF_2... \n", + "1593374 S2B_MSIL2A_20220719T095559_N0400_R122_T33TVF_2... \n", + "1595266 S2B_MSIL2A_20220719T095559_N0400_R122_T33TVF_2... \n", + "1595265 S2A_MSIL2A_20200113T095351_N0500_R079_T33TVF_2... \n", + "1595264 S2A_MSIL2A_20230318T095031_N0509_R079_T33TVF_2... \n", + "\n", + " timestamp cloud_cover nodata centre_lat centre_lon \\\n", + "1591452 2021-05-10 10:00:31 0.000000 0.0 40.734109 14.155341 \n", + "1593374 2022-07-19 09:55:59 0.000000 0.0 40.823861 14.292709 \n", + "1595266 2022-07-19 09:55:59 0.000000 0.0 40.913671 14.311585 \n", + "1595265 2020-01-13 09:53:51 0.000000 0.0 40.913731 14.192730 \n", + "1595264 2023-03-18 09:50:31 2.408506 0.0 40.913791 14.073875 \n", + "\n", + " crs parquet_url \\\n", + "1591452 EPSG:32633 https://huggingface.co/datasets/Major-TOM/Core... \n", + "1593374 EPSG:32633 https://huggingface.co/datasets/Major-TOM/Core... \n", + "1595266 EPSG:32633 https://huggingface.co/datasets/Major-TOM/Core... \n", + "1595265 EPSG:32633 https://huggingface.co/datasets/Major-TOM/Core... \n", + "1595264 EPSG:32633 https://huggingface.co/datasets/Major-TOM/Core... \n", + "\n", + " parquet_row geometry \n", + "1591452 33 POINT (14.155 40.734) \n", + "1593374 455 POINT (14.293 40.824) \n", + "1595266 347 POINT (14.312 40.914) \n", + "1595265 346 POINT (14.193 40.914) \n", + "1595264 345 POINT (14.074 40.914) " + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "filtered_df = filter_metadata(gdf,\n", + " cloud_cover = (0,10), # cloud cover between 0% and 10%\n", + " region=box(14.011710578,40.7015558593,14.423765416,41.1019258062), # you can try with different bounding boxes, like in the cell above\n", + " daterange=('2020-01-01', '2025-01-01'), # temporal range\n", + " nodata=(0.0,0.0) # only 0% of no data allowed\n", + " )\n", + "\n", + "filtered_df.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 9 unique grid cells in metadata.\n", + "Filtered grid points to 9 matching records.\n" + ] + } + ], + "source": [ + "filtered_points = grid.filter_gridpoints_from_metadata(filtered_df)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Subset of grid points with footprints:\n", + " grid_cell utm_crs utm_footprint\n", + "1638 453U_119R EPSG:32633 POLYGON ((434017.794305 4504244.40132, 434017....\n", + "1637 454U_120R EPSG:32633 POLYGON ((445697.420125 4514104.878801, 445697...\n", + "1633 455U_118R EPSG:32633 POLYGON ((427349.551347 4524259.936869, 427349...\n", + "1634 455U_119R EPSG:32633 POLYGON ((437358.568485 4524154.14666, 437358....\n", + "1635 455U_120R EPSG:32633 POLYGON ((447367.539888 4524061.936163, 447367...\n", + "1631 456U_118R EPSG:32633 POLYGON ((429020.402583 4534213.513238, 429020...\n", + "1632 456U_119R EPSG:32633 POLYGON ((439029.079893 4534109.665314, 439029...\n", + "1629 457U_118R EPSG:32633 POLYGON ((430691.345607 4544167.511129, 430691...\n", + "1631 457U_120R EPSG:32633 POLYGON ((450707.928244 4543977.390889, 450707...\n" + ] + } + ], + "source": [ + "subset_with_footprints = generate_product_outlines(\n", + " filtered_points, get_footprints=True\n", + ")\n", + "\n", + "print(\"Subset of grid points with footprints:\")\n", + "print(subset_with_footprints[[\"grid_cell\", \"utm_crs\", \"utm_footprint\"]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Conclusion\n", + "### This notebook shows how to get the actual product polygons of the dataset from the grid class without having to download the data, including generating outlines, working with UTM zones, saving files, and retrieving utm footprints.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PytKerasNewHQ", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/grid.py b/src/grid.py index f27b644..9a62a2d 100644 --- a/src/grid.py +++ b/src/grid.py @@ -1,8 +1,9 @@ import numpy as np +import os import math import pandas as pd import geopandas as gpd -from shapely.geometry import LineString, Polygon +from shapely.geometry import LineString, Polygon, box from tqdm import tqdm import re @@ -105,14 +106,14 @@ def get_points(self): c_idx += 1 points_by_row[r_idx] = gpd.GeoDataFrame({ - 'name':point_names, + 'grid_cell':point_names, 'row':grid_row_names, 'col':grid_col_names, 'row_idx':grid_row_idx, 'col_idx':grid_col_idx, 'utm_zone':utm_zones, - 'epsg':epsgs - },geometry=gpd.points_from_xy(grid_lons,grid_lats)) + 'utm_crs':epsgs + },geometry=gpd.points_from_xy(grid_lons,grid_lats), crs='EPSG:4326') r_idx += 1 points = gpd.GeoDataFrame(pd.concat(points_by_row)) # points.reset_index(inplace=True,drop=True) @@ -202,7 +203,238 @@ def get_bounded_footprint(self,point,buffer_ratio=0): bbox = Polygon([(new_left,new_bottom),(new_left,new_top),(new_right,new_top),(new_right,new_bottom)]) return bbox + + + + def generate_product_outlines_for_utm_zone(self, utm_zone, shift=340, pixel_size=10, raster_width=1068, raster_height=1068, + output_file=None, driver="ESRI Shapefile", get_footprints=False): + """ + Generate the Major-TOM grid and product outlines for a given UTM zone. + + Args: + utm_zone (str): UTM zone (e.g., "EPSG:32633"). + shift (int): Shift applied to the bottom-left corner (default is 340 meters). + pixel_size (int): Pixel size in meters (default is 10). + raster_width (int): Raster width in pixels (default is 1068). + raster_height (int): Raster height in pixels (default is 1068). + output_file (str): Path to save the resulting shapefile or GeoJSON. + driver (str): File driver (e.g., "ESRI Shapefile", "GeoJSON"). + get_footprints (bool): Whether to include the `utm_footprint` column. + + Returns: + GeoDataFrame: Product outlines for the specified UTM zone. + """ + grid_points_utm = self.points[self.points['utm_crs'] == utm_zone] + product_outlines = generate_product_outlines( + grid_points_utm, + shift=shift, + pixel_size=pixel_size, + raster_width=raster_width, + raster_height=raster_height, + get_footprints=get_footprints, + ) + + if output_file: + output_file = ensure_file_extension(output_file, driver) + product_outlines.to_file(output_file, driver=driver) + return product_outlines + + + + def get_product_outline_for_cell(self, grid_cell=None, lat_lon=None, shift=340, pixel_size=10, raster_width=1068, raster_height=1068): + """ + Retrieve the product outline for a single grid cell. + + Args: + grid_cell (str): The grid cell name (optional). + lat_lon (tuple): A tuple of (latitude, longitude) (optional). + shift (int): Shift applied to the bottom-left corner (default is 340 meters). + pixel_size (int): Pixel size in meters (default is 10). + raster_width (int): Raster width in pixels (default is 1068). + raster_height (int): Raster height in pixels (default is 1068). + + Returns: + shapely.geometry.Polygon: Product outline for the specified grid cell. + """ + if grid_cell: + grid_row = self.points[self.points['grid_cell'] == grid_cell] + elif lat_lon: + lat, lon = lat_lon + grid_row_idx = self.latlon2rowcol([lat], [lon], return_idx=True)[-1][0] + grid_row = self.points.iloc[[grid_row_idx]] + else: + raise ValueError("Either 'grid_cell' or 'lat_lon' must be provided.") + + return generate_product_outlines(grid_row, shift=shift, pixel_size=pixel_size, raster_width=raster_width, raster_height=raster_height).iloc[0].geometry + + def generate_global_product_outlines_by_utm(self, output_folder, shift=340, pixel_size=10, raster_width=1068, raster_height=1068, + naming_convention="zone", driver="ESRI Shapefile", get_footprints=False): + """ + Generate global product outlines grouped by UTM zones. + + Args: + output_folder (str): Directory to save the shapefiles. + shift (int): Shift applied to the bottom-left corner (default is 340 meters). + pixel_size (int): Pixel size in meters (default is 10). + raster_width (int): Raster width in pixels (default is 1068). + raster_height (int): Raster height in pixels (default is 1068). + naming_convention (str): Naming convention for filenames ("epsg" or "zone"). + - "epsg": Use EPSG code in the filename. + - "zone": Use UTM zone with hemisphere in the filename. + + Returns: + None + """ + os.makedirs(output_folder, exist_ok=True) + + for utm_zone in self.points['utm_crs'].unique(): + print("Processing utm_zone... ", utm_zone) + # Determine naming convention + if naming_convention == "zone": + if "EPSG:326" in utm_zone: + zone_number = int(utm_zone.split(":326")[1]) # UTM Zone = EPSG - 32600 + hemisphere = "N" # Northern Hemisphere + elif "EPSG:327" in utm_zone: + zone_number = int(utm_zone.split(":327")[1]) # UTM Zone = EPSG - 32700 + hemisphere = "S" # Southern Hemisphere + else: + print(f"Skipping EPSG code: {utm_zone} (not a valid UTM zone)") + continue + output_filename = f"raster_outlines_UTM_zone_{zone_number}{hemisphere}.shp" + else: + # Default to using EPSG in the filename + output_filename = f"raster_outlines_{utm_zone.replace(':', '_')}.shp" + + # Construct full path to output file + output_file = os.path.join(output_folder, output_filename) + + # Generate the product outlines for the current UTM zone + self.generate_product_outlines_for_utm_zone( + utm_zone, + shift=shift, + pixel_size=pixel_size, + raster_width=raster_width, + raster_height=raster_height, + output_file=output_file, + driver=driver, + get_footprints=get_footprints + ) + + def filter_gridpoints_from_metadata(self, metadata): + """ + Filter grid_points using the metadata.parquet file. + + Args: + metadata: Dataframe with the contents of the metadata file (Parquet format) containing a 'grid_cell' column. + Returns: + Filtered grid points that match the metadata rows + """ + + if 'grid_cell' not in metadata.columns: + raise ValueError("The metadata file must contain a 'grid_cell' column.") + + unique_grid_cells = metadata['grid_cell'].unique() + print(f"Found {len(unique_grid_cells)} unique grid cells in metadata.") + + # Filter the grid points to include only the matching grid cells + filtered_grid_points = self.points[self.points['grid_cell'].isin(unique_grid_cells)] + print(f"Filtered grid points to {len(filtered_grid_points)} matching records.") + + return filtered_grid_points + + + +def generate_product_outlines(grid_points_gdf, shift=340, pixel_size=10, raster_width=1068, raster_height=1068, get_footprints=False): + """ + Optimized generation of product outlines from grid points by processing all points in a UTM zone together. + + Args: + grid_points_gdf (GeoDataFrame): Grid points GeoDataFrame for a single UTM zone. + shift (int): Distance in meters to shift the bottom-left point (default is 340m). + pixel_size (int): Pixel size in meters (e.g., 10m for Sentinel-2). + raster_width (int): Raster width in pixels (e.g., 1068). + raster_height (int): Raster height in pixels (e.g., 1068). + + Returns: + GeoDataFrame: GeoDataFrame with raster outlines as polygons. + """ + # Reproject all points in the GeoDataFrame to the UTM CRS + utm_crs = grid_points_gdf["utm_crs"].iloc[0] # Get the CRS (all points in this group share the same CRS) + reprojected_gdf = grid_points_gdf.to_crs(utm_crs) + # Shift the bottom-left points + reprojected_gdf["shifted_x"] = reprojected_gdf.geometry.x - shift + reprojected_gdf["shifted_y"] = reprojected_gdf.geometry.y - shift + + # Vectorized bounding box creation + reprojected_gdf["geometry"] = reprojected_gdf.apply( + lambda row: box( + row["shifted_x"], + row["shifted_y"], + row["shifted_x"] + raster_width * pixel_size, + row["shifted_y"] + raster_height * pixel_size, + ), + axis=1, + ) + if get_footprints: + reprojected_gdf["utm_footprint"] = reprojected_gdf["geometry"].astype('string') + export_columns = ["grid_cell", "utm_crs", "geometry", "utm_footprint"] + else: + export_columns = ["grid_cell", "utm_crs", "geometry"] + # Drop temporary columns and return the GeoDataFrame + return reprojected_gdf[export_columns] + +def merge_utm_files_to_wgs84(utm_folder, output_file, ext=".shp", driver="ESRI Shapefile"): + """ + Merge UTM zone shapefiles into a single WGS84 file. + + Args: + utm_folder (str): Directory containing UTM zone shapefiles. + output_file (str): Path to save the merged WGS84 file. + + Returns: + None + """ + all_gdfs = [] + for file in os.listdir(utm_folder): + if file.endswith(ext): + print("processing file: ", file) + gdf = gpd.read_file(os.path.join(utm_folder, file)) + gdf = gdf.to_crs("EPSG:4326") # Reproject to WGS84 + all_gdfs.append(gdf) + + merged_gdf = gpd.GeoDataFrame(pd.concat(all_gdfs, ignore_index=True)) + output_file = ensure_file_extension(output_file, driver) + merged_gdf.to_file(output_file, driver=driver) + + +def ensure_file_extension(output_file, driver): + """ + Ensure the correct file extension for the given driver. + + Args: + output_file (str): The output file path. + driver (str): The file driver (e.g., "ESRI Shapefile", "GeoJSON"). + + Returns: + str: The file path with the correct extension. + """ + extension_map = { + "ESRI Shapefile": ".shp", + "GeoJSON": ".geojson", + "GPKG": ".gpkg", + } + desired_extension = extension_map.get(driver, "") + if not desired_extension: + raise ValueError(f"Unsupported driver: {driver}") + + base, ext = os.path.splitext(output_file) + if ext.lower() != desired_extension: + output_file = f"{base}{desired_extension}" + return output_file + + + def get_utm_zone_from_latlng(latlng): """ Get the UTM zone from a latlng list and return the corresponding EPSG code. @@ -272,6 +504,7 @@ def get_utm_zone_from_latlng(latlng): print(test_lats[:10]) print(test_rows[:10]) print(test_cols[:10]) + # Make line segments from the points to their corresponding grid points lines = [] @@ -282,3 +515,12 @@ def get_utm_zone_from_latlng(latlng): lines.to_file(f'testlines_{dist}km.geojson',driver='GeoJSON') grid.points.to_file(f'testgrid_{dist}km.geojson',driver='GeoJSON') + + # Test 2: Single Grid Cell Outline + grid_cell_name = "1U_17R" # Replace with a valid grid cell name + outline = grid.get_product_outline_for_cell(grid_cell=grid_cell_name) + print(f"Product outline for grid cell {grid_cell_name}: {outline}") + + lat_lon = (45.0, 13.0) # Replace with a valid lat/lon pair + outline = grid.get_product_outline_for_cell(lat_lon=lat_lon) + print(f"Product outline for lat/lon {lat_lon}: {outline}")