plot_data_sideview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8147690-cf50-452d-a62e-aa53ecd389fc",
   "metadata": {},
   "source": [
    "# Plotting function test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "34d81763-e3d9-4b7f-ae74-adf6ae6e698f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "import numpy as np\n",
    "from osgeo import gdal\n",
    "from pathlib import Path\n",
    "import rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e40ca61-3dc5-4b1f-a9ce-fa0ce9f37cec",
   "metadata": {},
   "source": [
    "### Data extraction from GeoTIFFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "40acdb7e-cc11-4152-b4d8-178c428c9ac4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the location as lon,lat\n",
    "# Bounding box: [-179.5,-50,5,85]\n",
    "lon = -110\n",
    "lat =   55"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "7a270634-db16-4e67-89bf-ea42dfd40e63",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Files\n",
    "base_path = Path('/Users/wmk934/data/perceptual_models/data/geotiff_same_resolution/res_0.01/')\n",
    "file_paths = [['climate_fs', base_path / 'worldclim/wc2.1_30s_climate_index_fs_resampled.tif'],\n",
    "              ['climate_im', base_path / 'worldclim/wc2.1_30s_climate_index_im_resampled.tif'],\n",
    "              ['climate_imr', base_path / 'worldclim/wc2.1_30s_climate_index_imr_resampled.tif'],\n",
    "              ['lai_01', base_path / 'lai/2013_2023_01_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_02', base_path / 'lai/2013_2023_02_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_03', base_path / 'lai/2013_2023_03_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_04', base_path / 'lai/2013_2023_04_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_05', base_path / 'lai/2013_2023_05_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_06', base_path / 'lai/2013_2023_06_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_07', base_path / 'lai/2013_2023_07_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_08', base_path / 'lai/2013_2023_08_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_09', base_path / 'lai/2013_2023_09_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_10', base_path / 'lai/2013_2023_10_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_11', base_path / 'lai/2013_2023_11_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['lai_12', base_path / 'lai/2013_2023_12_MOD_Grid_MOD15A2H_Lai_500m_resampled.tif'],\n",
    "              ['veg_height', base_path / 'forest_height/forest_height_2020_resampled.tif'],\n",
    "              ['land_glclu', base_path / 'glclu2019/glclu2019_map_resampled.tif'],\n",
    "              ['land_igbp', base_path / 'modis_land/2013_2022_mode_MCD12Q1_LC_Type1_resampled.tif'],\n",
    "              ['agriculture', base_path / 'lgrip30/lgrip30_agriculture_resampled.tif'],\n",
    "              ['aspect', base_path / 'merit/merit_hydro_aspect_resampled.tif'],\n",
    "              ['slope', base_path / 'merit/merit_hydro_slope_resampled.tif'],\n",
    "              ['elev', base_path / 'merit/merit_hydro_elv_resampled.tif'],\n",
    "              ['clay1', base_path / 'soilgrids/clay_0-5cm_mean_resampled.tif'],\n",
    "              ['clay2', base_path / 'soilgrids/clay_5-15cm_mean_resampled.tif'],\n",
    "              ['clay3', base_path / 'soilgrids/clay_15-30cm_mean_resampled.tif'],\n",
    "              ['clay4', base_path / 'soilgrids/clay_30-60cm_mean_resampled.tif'],\n",
    "              ['clay5', base_path / 'soilgrids/clay_60-100cm_mean_resampled.tif'],\n",
    "              ['clay6', base_path / 'soilgrids/clay_100-200cm_mean_resampled.tif'],\n",
    "              ['sand1', base_path / 'soilgrids/sand_0-5cm_mean_resampled.tif'],\n",
    "              ['sand2', base_path / 'soilgrids/sand_5-15cm_mean_resampled.tif'],\n",
    "              ['sand3', base_path / 'soilgrids/sand_15-30cm_mean_resampled.tif'],\n",
    "              ['sand4', base_path / 'soilgrids/sand_30-60cm_mean_resampled.tif'],\n",
    "              ['sand5', base_path / 'soilgrids/sand_60-100cm_mean_resampled.tif'],\n",
    "              ['sand6', base_path / 'soilgrids/sand_100-200cm_mean_resampled.tif'],\n",
    "              ['silt1', base_path / 'soilgrids/silt_0-5cm_mean_resampled.tif'],\n",
    "              ['silt2', base_path / 'soilgrids/silt_5-15cm_mean_resampled.tif'],\n",
    "              ['silt3', base_path / 'soilgrids/silt_15-30cm_mean_resampled.tif'],\n",
    "              ['silt4', base_path / 'soilgrids/silt_30-60cm_mean_resampled.tif'],\n",
    "              ['silt5', base_path / 'soilgrids/silt_60-100cm_mean_resampled.tif'],\n",
    "              ['silt6', base_path / 'soilgrids/silt_100-200cm_mean_resampled.tif'],\n",
    "              ['soil_depth', base_path / 'pelletier/average_soil_and_sedimentary-deposit_thickness_resampled.tif']\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "42a10b8d-5c03-49c3-8dc7-badcb6c6059c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_values = get_location_data_dict_from_geotiffs(file_paths,lon,lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d7be07b-5175-436c-ab44-697b1cf108a2",
   "metadata": {},
   "source": [
    "### Define variables of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "8d896839-26a4-4e93-b6fd-8b3184189a82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forest height (GLCLU2020)\n",
    "vegetation_height    = data_values['veg_height'] # [m]\n",
    "\n",
    "# GLCLU2019\n",
    "land_type1           = data_values['land_glclu'] # Categorical\n",
    "\n",
    "# MODIS-IGBP\n",
    "land_type2           = data_values['land_igbp'] # Categorical\n",
    "\n",
    "# LGRIP30\n",
    "land_agriculture     = 1 # Categorical. 0: water, 1: non-cropland, 2: irrigated, 3: rainfed\n",
    "\n",
    "# Modis\n",
    "vegetation_lai       = np.array([data_values['lai_01'],data_values['lai_02'],data_values['lai_03'],\n",
    "                                 data_values['lai_04'],data_values['lai_05'],data_values['lai_06'],\n",
    "                                 data_values['lai_07'],data_values['lai_08'],data_values['lai_09'],\n",
    "                                 data_values['lai_10'],data_values['lai_11'],data_values['lai_12']])\n",
    "\n",
    "# Merit\n",
    "topography_elevation = data_values['elev'] # [m.a.s.l.]\n",
    "topography_slope     = math.tan(math.radians(data_values['slope'])) # [degrees] to [m m-1], derived\n",
    "topography_aspect    = data_values['aspect'] # [degrees], derived\n",
    "\n",
    "# Pelletier\n",
    "soil_depth           = data_values['soil_depth'] # [m]\n",
    "\n",
    "# Soilgrids\n",
    "soil_layer_top       = [0.00, 0.05, 0.15, 0.30, 0.60, 1.00] # Top of layer depths [m]\n",
    "soil_layer_bot       = [0.05, 0.15, 0.30, 0.60, 1.00, 2.00] # Bottom of layer depths in [m]\n",
    "soil_sand            = [data_values['sand1']/1000, data_values['sand2']/1000, data_values['sand3']/1000,\n",
    "                        data_values['sand4']/1000, data_values['sand5']/1000, data_values['sand6']/1000]\n",
    "soil_silt            = [data_values['silt1']/1000, data_values['silt2']/1000, data_values['silt3']/1000,\n",
    "                        data_values['silt4']/1000, data_values['silt5']/1000, data_values['silt6']/1000]\n",
    "soil_clay            = [data_values['clay1']/1000, data_values['clay2']/1000, data_values['clay3']/1000,\n",
    "                        data_values['clay4']/1000, data_values['clay5']/1000, data_values['clay6']/1000]\n",
    "\n",
    "# WorldCLim\n",
    "climate_aridity      = data_values['climate_im']\n",
    "climate_seasonality  = data_values['climate_imr']\n",
    "climate_snow         = data_values['climate_fs']\n",
    "\n",
    "# Currently unused\n",
    "geology = 0 # Glhymps \n",
    "lake = 0 # # Hydrolakes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "660163a9-6c4e-45d4-babb-e42b4f3fee1d",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "4e5eb5a7-700b-49f2-930e-f17b44dc6cd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pm_plot_color(object):\n",
    "\n",
    "    '''Defines standard colors for objects in perceptual model plots'''\n",
    "\n",
    "    if object == 'bedrock boundary': color = [0.3, 0.3, 0.3] # Dark-ish grey\n",
    "    elif object == 'bedrock fracture': color = [0.7, 0.7, 0.7] # Light-ish grey\n",
    "    elif object == 'surface': color = 'saddlebrown' # Brown\n",
    "    elif object == 'soil':\n",
    "        print('WARNING: soil is defined based on sand/silt/clay percentages and requires get_pm_soil_color()')\n",
    "        color = None\n",
    "\n",
    "    return color"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "98a51e69-634c-4f93-9743-1f9f1c737f5e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALEAAAGRCAYAAAAq8lJHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9WUlEQVR4nO2deZgU1bn/v9X7Mr3MzsAAgyibMOworle4wnUjGLmgD0SiUa8ETfKo13jBJQmJGn/R5zHRh1wNos51u4oGVxblRtEEiDAiywwDw84w+0zvW3Wd3x/NKbqnq2bprpnprj6f55kHput0LVPfOvWe97znfTlCCAGDkcVoBvsEGIx0YSJmZD1MxIysh4mYkfUwETOyHiZiRtbDRMzIepiIGVkPEzEj69H1tqEgCPjyyy/xj3/8A8ePH0c4HIZOp4PJZILRaEQgEEBHRweWLl2KG2+8Men7Z86cwWuvvYZ9+/ahra0NhYWFmDt3LhYvXgy73a7oRTFyDNJL7rrrLgKgx59169YlfXf79u3EbrdLth89ejQ5e/Zsb0+DwUiCI6R3sRObN2/Gpk2bYDQaUVRUhLKyMpSVlSEajeKWW26Bx+OB0+nEyZMnYbPZxO+1traisrISZ8+ehdlsxl133YX8/Hzs2rUL27ZtQzgcxrx587Bp0yZwHKfg48nIGdJ5AiKRCLn99tsJAGI0GsnWrVuT2jz++OMEANFqteTLL79M2LZnzx6i0WgIALJ58+Z0ToWRw6QsYrfbTa677joCgFgsFvLFF18ktYlGo2TEiBEEAHnwwQcl97Nw4UICgNx2222pngojx+n1wC6eQCCABQsW4G9/+xsA4IMPPsCcOXOS2u3YsQMnT54Ex3F48MEHJfd17bXX4q9//SsOHz4suX316tUAALPZDLvdDpPJBL1eD6vVitmzZ8PpdMJgMIDjOOj1emaS5CApifihhx4SBazT6bBjxw7885//REVFBX74wx/CbDYDAL7//nsAwPTp01FWVia5r4qKCgCA3++X3P70009DEATJbTt37kR5eTn8fj8EQYBer0/lchhZTkoi3rdvn/h/nufxxBNPiL9XVFRgw4YNmDZtmti7Tp48WXZfPM/HTkSXfCqEEFkBAzG3n8/nEweSPp8PFoulbxfDyHpSEvErr7yC1157Da2trTAYDNBoNKivr8enn36K48eP484778SePXtw9uxZAEBxcbHsvlpbWwEARUVFSdui0Wi350EISTAfemrPUCcpifjCCy/EmjVrkj5/9dVXcccdd2Dv3r3YtWsXIpEIACAUCsnua+fOnQDOmxVdWbZsGXieh9/vh8fjQTAYBM/zcLlc0Gq10Gq1Ylsm4tykzyIOhULiQKory5cvx6OPPoozZ87gH//4h/hqP3PmjOS+BEHApk2bAACzZs1KPjmdDlVVVd2eT7xwqWnCyC36FDuxa9cuWK1W/PSnP5XcTggRbVuTySQO5rZs2YJAIJDU/ssvv8TJkyeh1WqxYMGCvp47ACT0xEzEuUmfRHzw4EFEo1G89NJLOHHiRNL2P/3pTzhx4gQ0Gg3mzZuHGTNmAAA6Ozuxdu3ahLaCIOC3v/0tAOCWW26R9V70BvrgUPOFkWP0xal85swZYrVaCQAyZ84ccvToUSIIAjl69ChZuXIl4TiOACArVqwghMRm9GbOnEkAEI1GQ/74xz+SYDBIeJ4nDz/8sDiTd/DgwbSc3WazmQAg9fX1ae2HkZ30ecbunXfeIQaDQQzg0ev1CQE9119/PQkGg2L7M2fOkEmTJonbjUajKDoA5Lnnnkv7IgoKCggAcuDAgbT3xcg+Upp23rt3L1m4cCHR6XSiGJ1OJ/n1r39NIpFIUnuv10seeughkp+fL7YvLi4mb7zxRtoXQAghw4cPJwDIrl27FNkfI7vodRSbFD6fD/X19dBqtbjwwgthNBq7bR8KhVBTUwOe51FZWQmDwZDqoRO46KKLcOTIEWzfvh1XXHGFIvtkZA8p+YkpVqsVlZWVvW5vNBoxZcqUdA4pex5A7KFi5B6qWJ5E/dFy8RcMdaMKEZtMJgBAMBgc5DNhDAaqEDHriXMbVYiYLjR1uVyDfCaMwUAVIi4pKQEANDY2DvKZMAYDVYi4sLAQQGx6m5F7qELE+fn5AIC2trZBPhPGYKAKEZeXlwOQD/lkqBtViHjYsGEAgNOnTw/ymTAGA1WImA7s6FInRm6hChHT9XmhUIi52XIQVYjYarWKvmLmZss9VCFiACgoKADA3Gy5iGpE7HA4ALBZu1xENSKm5oTb7R7kM2EMNKoRMc0C5PF4BvlMGAONakTMeuLcRTUizsvLAwB4vd5BPhPGQKMaEbP4idxFNSKmEx5MxLmHakTMbOLcRTUiZjZx7qIaEbNl+7mLakTMVjznLqoRMTMnchfViNjpdAJgsRO5iGpETCs2sdwTuYdqREyTGXZXH4ShTlQn4kgkgjQSfTKyENWIOL4OHqvdkVuoRsSsAE3uokoRd1eFlKE+VCNiRu6SVqb49vZ2PP3009ixYwd0Oh1GjRqFyspKzJ8/H+PGjRPbffrpp3j55ZcRiURgMpnE1cmCIIAQAovFgqFDh+KBBx5I+VziizJqNOzZzClSLfbxzTffkJKSkoTKSfE/Dz/8sNj2Zz/7mWw7+jNlypS0io+0traK+5IqfsNQLyn1xPv27cO8efPg8/ngdDqxcuVKjB07Fs3NzXj11Vexf/9+PPPMM1i6dCkqKyvFIom33norVq5ciY6ODhw/fhx1dXU4duwYTp06hbvuuiuthzG+EGO8p4KhflK627/73e/g8/mQn5+P7du34+KLLxa3LVq0SCw2fvDgwYTCNDfeeGO/VTeiItbr9f2yf0bmkrLxWFBQgE8//TRBwECil4CWuqXxDHQJUX9Ao9doNBsjd0ipJ37ppZeg0WjEyLF41q1bByC2hH769OkAIBYnb2xshNfrxebNm1FTU4NRo0Zh0aJF3da/W716NYBYbITdbofJZIJerxcHh06nEzNnzkQ4HAbAeuKcREkD+5133iFarZYAIL/5zW/EzxctWiSWxI2vQgqAjBgxgvz1r3+V3adGo+lxUNjS0kKqq6sJAFJaWqrkJTGyAEV8UZFIBA8//DCWLFmCaDSK66+/Hv/1X/8lbqd2MSEECxYswPPPP49HHnkExcXFOHnyJG677TacOnVK6gHr1cRFXl6eGL1GV3gwcoh0nwKfz0fmz58v9orLli1LKFBOCCGBQIDs2rWLuN3uhM8bGxtJRUUFAUDWrFmTtO9IJNJjL2w0GgkhhHzyyScEAJk2bVq6l8TIMtLyRfl8Plx//fX46quvAAAPPfQQfv/73ydNNphMJsycOTPp+6WlpXjkkUdw7733SvbEALBs2TLwPA+/3w+Px4NgMAie5+FyueDz+cSBJF2qT7NjMnKHlEUcDAbxgx/8QBTwo48+it/85jfgOK5P+7n00ksBxGb/kk5Op0NVVVW336fBPgPhAWFkJinZxNFoFEuXLsUXX3wBAHjqqaewZs0aSQETQrpN8tfc3AwAMBgMqZyKOLFBHwK6TImRO6Qk4l/84hd4//33AQC/+tWv8Mgjj8i2/fDDDzFs2DDs2rVLcvtnn30GAEn+5r5Ck6YwEecefRbxt99+ixdeeAEAcO2112LVqlUQBAEejwcnTpzAjh07sG3bNjQ0NACIeQ48Hg9uuOEGbNmyRdxPIBDA22+/jbVr14LjONxwww1pXQhd5cy8E7lHn23i2tpa8f9bt26VNQP0ej127NiBuXPn4p577sFLL72E+fPno7S0FDabDcePHxft2YceegiTJ09O8RJi0KQpUhMwDHXT55745ptvxoMPPijWjpMjvkdcu3YtXn75ZZSVlaGpqQlHjhwBx3G45JJL8Nprr+GZZ57p+5l3oaWlBQAb2OUiHCGpr6oMBoPo6OiAVquFwWCARqNBMBiEIAgoLi5OiKMAYisuTpw4ASAWV6FknMP06dOxZ88efPzxx2mbJozsIi0/sclkEoN8KDQ7pRQajQajRo1K55CysBm73EU1SyA6OjoAnK+ixMgdVCFiQog4Y1dcXDzIZ8MYaFQh4o6ODtHTQTPGM3IHVYiYeiYcDgcLis9BVCFiOq1Na9kxcgtViJjO1rGJjtxEFSKmEWzMM5GbqELENIKNzdblJqoQcWdnJwAm4lxFFSKmYZjdzRYy1IsqREx7YmYT5yaqEDFbmpTbqELEzJzIbVQhYmpOMBHnJqoQMVtfl9uoQsRsxi63UYWIaU/MYidyE1WImAUA5TaqEDFNHUtL4zJyi6wXMSGEZYnPcbJexPG1OrpL1s1QL1kv4vjqoazgTG6S9SKOT5vB6tflJuyuM7IeVYmY1XTOTbJexPF2cHxpXEbukPUijs/KScuAMXKLrBcxx3HigC7e3cbIHbJexMB5/3AoFBrkM2EMBqoQMc36w0Scm6hCxHS6mZkTuYkqREx7YlqknJFbqELENASTBsczcouURfzNN99g/vz5GDNmDCZPnowlS5bgqaeewpEjR2Tb33rrrbj44otRWlqKCRMm4P77708oZJMqVMTd1ctjqJi+1tEVBIE88cQTsrWWtVotWbduXcJ33njjDcJxnGR7m81G9u7dm1Zt33/9138lAEhVVVVa+2FkJ30W8bPPPisKcPbs2eSVV14hW7ZsIS+++CIpLy8nAEh+fj4JBAKEEELq6+uJxWIhAMjw4cPJyy+/TDZs2EB++ctfkhEjRhAAZMKECSQcDqd8EQsXLiQAyNq1a1PeByN76bOICwoKCAAyd+5cEgwGE7a9+eabosDPnj1LCCHk/vvvJwBIXl4eqa2tTWjf0NBArFYrAUBeffXVlC/iRz/6EQFAnnnmmZT3wche+mwTW61WzJkzBx988EFSEDr105pMJhQVFSEcDuPNN98EADzxxBMYO3ZsQvuysjIsWbIEAPDee+/19VREmE2c2/Q5iry+vl5yGRDP81i7di0AYM6cOdDpdNi2bRva2tqg0+nwk5/8RHJ/V111FV555RWcOnVKcvvq1asBxNbP2e12mEwm6PV6WK1WXHDBBZgxYwYsFguA82XAGLlFn0UsJeBQKIQ777wTu3btglarxZo1awAABw8eBADMnj1bNk8arYMn5+N9+umnZUMsL7/8cnz99ddivgnmYstN0l7Pc/ToUSxZsgTffvstgFgJ3GnTpgEADh8+DAAYP3687PfpSmWp9XGEkG5jhGk9Z1qAkf7OyC3SEvF3332HefPmoaWlBQaDAa+//rpo4wJAc3MzgO6zVdI2paWlSdt6ig/uWkWUmRO5Scoi3r17N+bOnQuXy4W8vDx8+OGHuOaaaxLakHPr37obcG3fvh0AZMvlLlu2DDzPw+/3w+PxIBgMgud5uFwuMc8EfUhaW1tTvRxGFpOSiPfv34958+bB5XLB4XDgk08+weWXX57UjnoNaFHyrkQiEWzatAkAcPXVVyefnE6Hqqqqbs+FEIKhQ4cCABoaGvp0HQx10GcX2/HjxzF//ny0t7cjPz8f27ZtkxQwAJSXlwMAtm7dKpatjeeTTz5BS0sLTCYTrr/++r6eCoBYUHxhYSGA8/WdGTlGX5zKHR0dZNy4ceLkxa5du7pt/8UXX4iTH/feey8RBEHc5vf7ycUXX0wAkBUrVqTg4j5PQ0MDAUA0Gg2JRCJp7YuRffRJxGvWrBFF+dhjj5GGhgZy+vRpcvDgQbJ9+3by7rvvkrfffps0NDQQQmJxFosXLxa/c8cdd5C9e/eS2tpacvPNN4uxE2fOnEnrInieF2MzGhsb09oXI/vok4jvvvtu2cCf+B+n00n8fj8hJNbjLlq0SLKdTqcjGzduVORCSkpKCACyZ88eRfbHyB76ZBP/7ne/w8qVKzFu3DgxEF2j0cBkMom/63Q6TJw4EVqtFkBspu1///d/UVVVhcsuuwwcx4HjOFx55ZXYuXMnFixYkJIZ1BU6uGtqalJkf4zsgSMkLg9UHyCxXjwhdVQ0Gk1YfSxFNBpFJBIRRa8U1113HTZt2oR169bhzjvvVHTfjMwmZT8x7VHjob1vd2i12l616ytDhgwBwHriXEQVy5MAiG42KVceQ92oRsS0miit38HIHVQjYhaOmbuoRsQski13YSJmZD2qETEzJ3IXJmJG1qMaEdPYYrpShJE7qEbEzCbOXVQjYpaPLXdRjYipOcFs4txDdSIOhUKsilKOoRoRxxegYcm2cwvViDg+Mo6VAsstVCPirmGhjNxBNSKOj+1PMc6fkaWoUsSsUHluoZq7He+RYKZFbqEaEccP5uLrPTPUj2pETBN8azSaflnDx8hcVCNimt/YaDQycyLHUI2I6QSHVBJwhrpRnYjjZ+4YuYFqRBwOhwEwEeciqhExtYlpIBAjd1CNiGluYrvdPshnwhhoVCNimjTF6XQO7okwBhzViLizsxPA+UxAjNxBNSKmRWeKiooG+UwYA41qRExt4oKCgkE+E8ZAoxoR02yYTMS5h2pETG1iJuLcQzUipvkmaCYgRu6geMyiIAiYO3cuNBoNPv/8c3Ach/Xr12PNmjVimQOr1Qq73Q5CCKLRKKxWK0aMGIF169alfFyaIZ4N7HIPxUVcV1eHv/3tbwBiBcrHjBmD77//HseOHev2e7NmzUrruFTEZWVlae2HkX0oLuL4TO1NTU0YM2aMuOri5z//OdasWYOOjg4cP34cdXV1OHbsGE6dOoVFixalfMxwOCy62KQKnTPUjeIijq+KxPM8gPOrLmbNmgWbzQabzYYRI0bgqquuUuSYbW1tIISA4zgUFxcrsk9G9qD4wC4+dsFoNAI47znIy8tT+nAAgObmZgCx4jNskWjuoXhPHL80iIZF0jDJ6upqDB06FB999BFqamowcuRIrFy5EhUVFbL7W716NYBYdJrdbofJZIJer4fVaoXD4cD8+fPR0tICgJkSOYvSJUoPHDgglr09fPgwIYSQpUuXypbQ1ev15Mknn5QvearRyH63rKyMEELI//zP/xAA5JprrlH6chhZgOLvXtorAkB+fj4A4F/+5V8AAGPHjsXjjz+O//u//8Pbb7+NmTNnIhKJYNWqVfjuu++kHrBukwNS+5eaE6wnzk0UNycOHz4MIJb7gYZF3nXXXfjRj34k2siUm2++GVdddRV27tyJN998E1OmTEnY3lNONZqTmHpEWARbbtJvInY6nQn2cVcBAzGb+T//8z+xaNEisTftyrJly8DzPPx+PzweD4LBIHieh8vlwujRowEALpcLAAuIz1UUF3FtbS2A82Vqe+Liiy8GIF3OVqfToaqqqsd9eDweAKwnzlUUt4nr6+sBnHenBQIBVFdX99g+ncAd2hNT84KRWyguYmoW0AWbn3zyCaZNm4Y333xTsv37778PAJg0aVLKx6QiZkuTchPFzYmuq47HjBkDjUaD5cuX49SpU1i+fDnsdjtqa2uxYcMGvPLKKzCZTPjhD3+Y8jGpiJk5kaMo7bObMGECAUBuv/128bMXXniBcBwn6evlOI6sXbs2rWNOnDiRACBbt25N9/QZWYji5sSqVatw1VVX4ec//7n42cqVK7Fjxw7Mnz9fjJ2YNGkS7r77buzcuRP33ntvWsdkiVNyG46Q7E+rXlFRgRMnTmDHjh245JJLBvt0GAOMKhL50lUdtKpoPM2NZ+B2dcJZUAhnfiF0OpZwUG2oQsR0MBkfBkopKimD2WKFz+tBXc338Lg7EfD7EA6FoDcYoNFoYTSaYMmzoWzocDjzi6BnZklWoQpzQq/Xg+d5nDp1CuXl5SntQxAEuF0daDp7Gq7OdoSCAQhCFOFwGFqtFnq9EQaDEda8PDjyC1FSOhRGE8v7lglkvYh5nhdzEre2tibNFLo62+HzupFnc8DuyE/7eJFIGF6PC25XJ9pbm2K9ejiMSDgEnU4PrV4HLaeF3miE1ZoHq82BPJsdeXl2mMyDv4hVEAQxCXk0ykOj0YLjuKxOTJ71Ig4EAuIKZ5fLlRQ/8egDP8a+6p0AAIPRhILCYhQUlSA/vxjOgiKUDR2OotIyFBaVIr+wGEXFQxS7oZFwGF6vCz6vB16PGz6vG+FQCOFICHwkgmg0CkGIQohGY4tmhShAAEIEaLU6cBwHAgIOsX+jER7gAEJiZR04jgPPR8BxQDQqwGAwxuqVnDt/IggwGE3QaLQwm83QaLUghCSMCzQaDQQh5vEkRDi3by0IESCcO7+oIJw7RwFCVAA4ABwHrVYLo8EEncEAs8mMPJsDzoIi2OzOAX0osl7EPp9PnOL2eDwJq0dOHKvDz+78QUJ7g8EIi9UGk9kMa54NeoMRer0hpgyOg4bTIM/hhNVqg93uhDO/EDaHEzZ7PuwOBywWGyxWK6w2BwyG5KCmXCYSDsPj7oTL1Y6Azwev14VIOAyO00Cj1UCn08NssaKopAylQ4YpdlxVidjr9SZ4KAghiETC4CMRCESA2Wztt6I0PB9BJBwGz0fEmnocx0Gr1UGn00NvMGT1K1tJBEGAq7MN+QXKrIdUhXeC0vV55DgOBoNxQHpMnU7fJ/ddwO9DZ2cb/D4vgn4/3O4O+L0euF0d8HjcCAb80Ol00OsN0Op00Ol04DgNIpEwojyf8LCEwyGEQyH4fG5EwiEIAkE0ykMQBAR8XoRCQfj93pgpACAaPV/AXYgK4DQacBzAcbG5L0IEaLSxY9NzMBhNMJstMFvzYLXmIc/mhM0eG2c484uQZ3fAZnPA5nDC4SyIvd1k0Gg0igkYUIGI4xeGSq0Ccbs60d7RBq1Gg0AggEDAD5/Xg1AoCK/Xg3AohEgkLNqRHMeBEALtuf1qNFoIQlQ0NbRaLXR6PQwGI0xmMyxmC8wWK8xmS8xzYc+HqRfZ6s0WK8yWZL82JeD3obmpAY0NJ9F09gw6O1rR1tIU+7e1Ge2tzfB53bLfH2xsdidKhgxDYVEpCoqKUVg8BAWFJXA481FcMhQjLxij2Jsp682JYDAoBht1dnYmBQHt/e5b/PLBewb0nHQ6Hax5NthsDuTl2VBQUISi4pJz/5aidEgZiopKUVxSmtbkSygUhMfdCa/HDXdnO9yuDng9Lrg62+Fxu+D3eRAI+BEKBhAKBhAOhxAKBREJhxGJhMV/ozwPQYiCP/dvX9BoNDCZLdBqY/0hOZc6QaPVQqvRQKvTQ6PRgAgChHODSrPFimkzL8ft9zyY8rXHo6qeWGo509ChqfmN04Hnebg6O+Dq7Oi2HcdxsNkcKCgswrDyESgtHYry4SNRXDwEZUPLUVI6pFtTyGg0wVg8BEXFQ5S+BPGtRv++1JMCAmi05//mVLx9JRwOpX+S58h6EceXwKXJWuKxWPsn14USEELgdnfC7e7E8WNHkrZzHIei4lIUF5civ6AQRUWlcObno6ioFIWFxSgoLIr19jYbjMbk2cp06Jq/Q6vVKjooVnKckvUipv5SQoikiKWmorMFQghamhvR0tzYY1uTyQyHMx95eXY48/NhszlgsVhhMplhMplhNJliPbfRCL3BAIPeAL3eAL1eD41WB61GA51eL4o3vgeO/Y014t9aEATRbOD5CPx+X9JbUIhGERWi50wV4dz3OfB8bIJl5qzLYM1TZiVO1osYiE07h8NhSRFrz42yI5HwIJzZwBEMBhBsDKAJDYN9Kr2iuGQIqt76VJF9qSLnE512plVFu2LM4t5Yrfi8HsX2pQoR03QAoZD0YMGksL3ISJ9AwI9oNPnNmQqqEjENyeyKUrYXQzkIIXC5OhXZlypE3DVxYVdYqdzMxO12KbIfVYiYRrEFAgHJ7ZkQAslIxu/zKrIfVYnY7/dLbjebmIgzEbmBeF/JDRGzikoZiVJuT1WImKavojnZumJm5kRGEpEZw/QVVYg4Pp5YCrbCOTNhPXEc1PsgN7CLj69gZA5SM6ypoAoR0/gIOT+xQSI3MmPwEYh8FYC+oAoR99wTM3MiE6ErTdJFFSKmNjHNBNQVZk5kKsqsx1CFiLvW7uiKLsXAbUb/Qtf0pYsqRExXOMv1xPErERiZA6dRZo2dKu4uFbGci42RmWg1yqwUUYWIi4qKAEC2AlOWr4VVLcyciGPIkNhCyfhCkPEQgYk4E1FqzZ4qRNzTtHNPRR0Zg4NOr4zrU3ER9/TqDofDaGpqkp2YSAVqE8sFAKl9fV22kpE98aFDh+BwOLBixYqEzwkhePnllzFhwgSYTCYMGTIEVqsVV155JT788MO0j0uj2CKRiGR4X4RXJuSPoSz6TOyJP//8c3g8HmzYsCHh86eeegr33HMPampqxJ5aEAR8/fXXWLhwIT7++OO0jhtfhFHKpAgFpWfyGIOLUrkyFBUx9dO2traKn1VXV+Oxxx4DAMyaNQs7duzAkSNH8NZbb2HOnDkghODOO+8Ua9Glgl6vF6eepSY85MwMxuDSXdLBvqCoiKmQCCFij/v8889DEAQMGTIEH3/8MS655BKMHj0at956Kz766COUlZWhpaUFf/7zn9M6Ns3B1tnZmbQtEJCeBGEMLlJF61NB0flYKiS9Xg+O4+Dz+fDee+8BiJkUxcWJ6TwtFguWLl2KP/zhD9i8eTN++ctfJu1z9erVAGIPiN1uh8lkgl6vh9Vqhd1ux5QpUzBkyBA4nU40NjbKiJiZE5mI1ZrBGYBoaOT27dvh8/lgMpmwePFiyba07lxjo3SqpqeffloyZSvlqaeewiOPPNJtYLySiToYyqFUijFFzQk6qKKFwg8dOgQAuPzyy0UPQlcKCgoASC+3J4R0K2DgfC9LTRkp151PoVW1DOXQaLSK5QNRVMRNTU0AgPz8WJWi+vp6AMCFF14o+x3ac0rlhujNJEXXQoxSPXHAz2ziTEOuU0sFRc2Jw4cPAzjfu7a3twNAQjGYrpw5cwYAMGyYdCGSZcuWged5+P1+eDweBIPBWP5flws+n0/0RtDSX1JTz15mTmQcSmZlUlTEdXV1AM4LiqYHlRpsUbZt2wYAGD16dPLJ6XSoqqrq9pi0t6bH7OhITGwtCAICAeZiyzRMChayVMycEAQBtbW1AM73xHQSgn7eFZ/Ph82bNwMArr322pSOS6cu5QZ2Ab+PRbFlIJZu6pX0FcVE3NraKk4qUEHR3vWbb77BkSPJmdDfffddeDweOJ1OzJ8/P63jy63u8HgytzhLLpNfUNhzo16imIjjY3npIG3evHlihZyf/vSnCR6IlpYWrFq1CgBw3333pZ30r6sdTmH2cGZiszl6btRLFBNxfOANFeTEiRPx+OOPAwC2bt2KK664AlVVVXj33Xfxb//2bzh79izKy8vx0EMPpX38kpISAMmB8cy9lplYFaylotjAbsSIEefqBAsYNWqU+PkTTzwBg8GARx99FP/85z9x++23i9uKiorwwQcfJJXtSgU6sOvaEwdZ8E9GkpEiLiwsxJNPPokjR45g4cKF4uccx2HVqlW48cYb8d///d/Yvn07wuEw5s6di8cffxylpaWKHF/OT8xm6zITZ36BYvtS1MUmFftAqaysxIsvvqjk4RKQS7Td0dEu1ZwxyGSki22wofPwXet2eNydg3A2jJ7oTeng3qI6EXeNnejsoaonY3BQ0iZWjYipOcGfK/5HYbN1mYnDka/YvlQj4vj1WvEpQ1nwT2aiZOJz1Yg4vhZxfE/c3t4q1ZwxyDBzQoL45d/xIvaznjgjUbKOimpELEeALRLNODQajaI5o1Uj4vjel8ZrAEAopFySFoYy6A3KrHKmqFLE1D7m+QgzJzIQpesKqkbE8UuZqH3sl8lXzBhc7HblItgAlYqY9sQsljgzMSq0ypmiGhHTUFCNRiOK2OtlIs5ElFzVAahQxPGTHj6WOT4jUSppCkU1IqazdPEiDodDcs0Zg4iSwT+AikRMA3/is8rwUWUqVjKURcmlSYCKRNw1ExAAyVzFjMFHydk6QEUipnHE8ZkWlarizlAWpTLEU1QjYrosKT7bEKvVkZnoFS5TrDoRx2eN51mZg4yE9cQy0EzzrCfOfLQK19pWnYhpEhUAIESZKu4MZVEygg1QkYhp0kK73S5+xkGZ2sEMZWHmhAw0B1t8IhZWmDwz0ShU01ncn6J7G0S6JtsGAK22X6o5MNJEo1H2DakaEbe1tQE4n6UeUP6JZ2QmqhExTSQYnxZLb1B2AMHITFQjYpohPr4nNihU7I+R2ahGxHSyI947YTYrG7fKUAZBUDZzvypETAgRC87E98TWbgreMAYPQVB2EkoVIg6FQmIUW3zVUqUKYDOUJT5DkxKoQsS+uAWh8S42g0K1gxnKovRiBVWJ2Gg0Qhc3L69kDlyGciidvV8VIm5tjeVbi7eHAcCiYNI6hnIoLeIBm9Lyer0Ih8MwmUwwm80JWXrShcZNxAf/AIDZYgXHcayOXYbBR7LQJn7qqadgs9lQWFgIq9UKrVYLp9MJu92OvLw8lJaWYvbs2SnvnxZGj48lBmLL95VO1MFIH6XjvAdExLRwucFgwKhRo2C32+FyueDxeODz+dDc3JxWHTsahul0OpO2WSzMzZZpKL32cUAjZB577DE8+uijAGLCO3ToEI4dO4bTp0/jmmuuSXm/dGAnVbld6UWJjPRRuiceEBHTikZjx44VP3M4HJg1axZmzZqV9v6llutTlM42w0ifrBzY0de9LoVlKatXrwYQW4pvt9thMpmg1+thtVoxfPhwzJ49u1sRm5mIMw6lc+QNiIipDfTRRx+ho6MDH330EQ4ePIiKigo88MADCTWgu/L0008npG2NZ968edi8ebMY/CNVmVTJtPoMZfC4XYrub0BETGfR1q9fj/Xr14uf19XVYcuWLbjtttvwxhtvJAmZECIrYOB8sA9d1dHVTxw7NhNxpqF0ve0B8U4sXrwYY8eOxb//+79jw4YNaG9vx8GDB7F8+XIAwFtvvYWtW7cmfa+n1cp0ZbPUqg6KjbnYMg6ly7INSE983XXX4brrrkv4LD8/H+vXr0dnZyc2btyIjz/+GPPmzUv67rJly8DzPPx+PzweD4LBIHieh8vlwsiRIwGcn+yQMicKC4uTPmMMLsFsFLEcHMfhvvvuw8aNG8XlRfHodDpUVVX1uB+5yQ5A+azkjPTheR4+n1cxU2/QYyeGDRsG4HxvmgrU+yEl4vyCopT3y+g/XJ3KFY7vdxG3tbXh9ddfl41f2LdvH4DzYk4F6p3oGjsBAEXFJSnvl9F/tLa2KLavfhfx1q1bsXz5cjz++ONJ2wghePXVVwEAU6dOTfkY/nO16qRm7JhNnJl4PMq52frdJr7ssstgs9nw29/+FvX19fjxj38Mu92OmpoavP/++/jss89QXFyMRYsWpXwMqVIHFKPRBJ1Op/hqAkZ6KFmKot9FPGLECLz33ntYsmQJ3nrrLbz11lsJ2x0OB6qqqhKWFfWV7kQcO0Y+2tqUe30x0kfJ+oID4p2YN28eDh06hGeffRZ///vfQQhBRUUFLrnkEtx2220oKkpv8CVVryMeh5OJONPweT2K7WvAXGwlJSX4/e9/3y/77qknNrMVHhlHe3urYvsadBdbukSjUXFq2iBTM5ites48OrPJxdbfxA/Y5FKGyvXQjMGjvY31xCJSNZ27YmLmRMahZPxE1os4HjkRG1n+iYyDiTiO3qxk1rPEghmHkt6JrBdxb1C6DCsjfZRcoqQqEcv1ykrXTWOkD8/ziq16znoRx9vBckH0OuadyEj8fmWmnrNexPGLT+XiI5Su1sNQhtC5Bb7pkvUi1mq14tq8sEwtZ62GFaDJREIhJmIAsdUhNHsQzVHcFdYTZyZB1hOfhy4Qjc9THA+rZ5eZhBTyUKji7tJgeBocz8gO/ApNeKhCxDTwR84mZqldM5NIRPp+9RVViJgG+MiKWOFqPQxliCq02kYVIu6pJ44onIWRoQxssiMOGuAjJ2KlU4kylKG7FGV9QRUipvkmaE62rvSUDosxOCg1VlGFiGm+ifZ26dUCSpecYmQWOSFipQYQDGVRyn+vChHTWh1y5gTLOZGZaLXKhAOoQsQ9zdgpNUfPUBal1j6qQsQ02bZcUkK5mArG4KLUihtViJgmX5FKDwuwnjhTUWrtoypEXFpaCgBoaGiQ3B5h3omMRKna26oQ8dChQwEAZ8+eldyuVMgfQ1l0Ci0bU4WIqYvN7XZLzgIpXXKKoQzMJo6DDuwIIZIeCiUzMDKUg9nEcZhMJnGtnZSv2K9wySmGMjCbOA6O40STorm5OWFbNMoz70SGYlFL4RmlKCsrA5As4g6ZqWjG4GI2W2SzmPYV1Yi4pCRWYKaxsTHhc5erYzBOh9EDSlZ6VY2Iqa+4pSUxI7zSFd0ZymBRsHD8gCZk8Hq9qKurQ1tbG/Lz8zF58mTF5s9piVxvl4Imfpl4CsbgkqdgkcwB6YmDwSDuv/9+FBYWYvr06Zg3bx5mzpyJIUOG4Mknn1QkOFouMF7JUlMM5bAomDN6QES8cuVKvPDCCwiHwzAajWKv2d7ejtWrV+OPf/xj2segdZ27itirYApRhnIYFMwZPSDFGF955RUAwOLFi3H27Fm43W4cOHAAK1asAAA88sgjOH36dFrHkQvHVCpBB0NZLOYssolffPFFAMCECRPw2muvwWQyib+/8MIL+PLLL3Hw4EH86U9/Squ6EjUnaLFySnu7dGRbV7RaHaLR3gXP97Ztb4tAdtdOq9XBbnfA7nDCZrNDr9fDmmeDzeaAzWaH2WyBXm+ATq+DTqeHTqc7574ywmA0Qq/TQ6vTQaNJ7K8IIRAEAdFoFHwkAj7KIxrlEeV58DyPcDgEj8cNv8+LQCCAQMAHn9cLv98Hj8eNSCSMcDiEUDAIn88rZn7v7TXbFLSJ+1XEra2t+OSTTwAAzzzzjChgikajwZIlS/DEE09g+/btkvtYvXo1AMBsNsNut8NkMsVupNUKu92OoqIizJgxQ1bEP7j5Vlxx5RyEwiH4fT643Z0IBPzweT0IhUIIBgPweT245l+vh9lswbqXn0fj2TMIBgMIh0JJ9voNNy3C/H/7AZ763SqcbTgle+0TJ03DQ7/8Nf7f7x/HkboaGIxGWCxWmM0WmExmWK15MFssKCkdituW/gS1B7+Hz++DLc8Oi9UKi8UKm80OZ35hkgAzEUEQEAwGYDAY0dLciPb2Vni9Hrg6OxDw++DxeuB2dcDtdkGv1+NHy+9V7Ngc6cf0OBs3bsTChQuRn5+PpqYmSU/ERx99hAULFuCiiy5CXV1d0natVtvt0m6z2Qy/3y/uZ+bMmdi1a5di1xCJRMDzEYRCIUT5CAoKiwBwIERAU+NZdHZ2QKPhoNFoodPpoNPpoDcYUFIyBBwXEx/N2tlXCCEIBAJoa2uD2+2GxxN78Px+P7xeL7xeL3w+H0KhEPhzPWg0GkU4HEYwGEQoFDp3/jwEQUh4IDmOg0ajAcdx4nlrtVrxR6/Xw2KxIC8vTxzHmM1mmEwm2Gw2GI1GGAwGsZ3VaoXZbIbVah3wBI792hNTUV5++eWyrjQ6yJMSKn3ldUcgEIAgCLI28QMPPICOjg7o9XqYTCY4HA6YTCZYrVaYTCYYDAbxZtGbYbfbxRtmMBig0+lgNJpgtDvihKlF2dBylA0tlzwvQgiCwWCCoEKhEHw+H4LBIILBoChAr9eL9vZ2uN1utLe3w+/3IxAIwO/3K5ZgZKDgOA4mk0n8+zmdTvHvmpeXB6fTCYPBAK1WC4vFgjlz5qT9pulXER87dgxArL6zHC5XzAVGzYF4epsvIhAIyCYVbG9vl10FnQp6vV78oT0ZEBNtNBoVe8RIJDLoOeAMBoMomPg8zsD5DkK0i8/14ukuqqVvD7okTG6hAhALoZ01axZsNlvKbyugn0VM3V3dhdwdPXoUAFBeLt2jLVu2DDzPw+/3w+PxIBgMgud5uFwu+Hw+NDU1wev1yk52FBQUQKPRnLN/g2n3bJFIpF97R61WC7PZDIvFAovFApvNhqKiItjtdthsNrGXy8vLE3+MRiP0er0o1HgToes4BIgJjfbyGo3m3GDQLAqJitnn88Hr9Ypvi0AggGAwKJo1kUgE4XAYfr9ffHu4XC74/X6EQiEEAgHZxbv0Wk0mU1oCBvpZxDQ8Um7tGwBs2bIFADB27FjJ71dVVXV7DJ7nwXGc+OR3FfFzzz0n/p/2EvRVTl/zgUBAtC/9fj/cbjcCgQBCoRDC4bDYswYCAVHEkUgkwc6kNiYVk16vF1+pRqNR/LFYLDCZTKJJQwVbUFAAm80Gh8OR8HolhODUqVPo7OxEJBKBw+HA6NGjJW/8gQMHxLefzWbDFVdcIfk3O336NPbu3Qsg9ma5+uqrE/ZHe26/34/a2lpotVoYDAZMnjxZjN3ueg/OnDkjPjjUZANiqcWomN1uN1wul/i2om+KdOlXEdN8EHv27AEhJOkP39raim3btgEArrvuupSOQR8UanN310tyHCf2cAUFBXC5XKiuroZOp4PFYknwcsQTjUbR3NwsCtVms4nZ6eMhhKCjo0McTJnNZjFEtCv19fVoamoSb+jw4cMlhXn06FHU1NQAiMVNT5o0SbJdc3OzKGCNRoOpU6dKDrB8Ph/2798v/l5ZWSnZW/M8j+rqajG/XXl5uaSACSHYvXu3GLNSWlqKGTNmiNupUB0Oh9j7cxwHu92OyspKyb9NX+lXEdPedf/+/di5cycuvfTShO3r169HJBLB8OHDcfXVV6d1LPpE04LlPQ0WIpEIdu/eLdrQF154oaSACSH49ttvxZtUUlKCmTNnSu6zpqZGNI+sViuuvPJKyXanT58WhanT6TBjxgxJYXZ2dqK2tlb8fcqUKZI9VzAYxHfffSf+Pn78eEnBRaNR7NmzRxxrDB8+XAxh7XrN+/btE/82TqcTF110keS1HD16VPzbGI1GTJ48WfYho9ei0WgwceJExbwY/eqAvPHGG8Ue64477kBr6/mi1DU1NfjVr34FIDZjF18FKRXivR892ayEEOzduzfhJo0ZM0ay7eHDh3t1k86ePSsKmOM4TJs2TfKa3G439u3bJ/4+adIkcVAaTzgcxu7du0VzZfTo0WJqgngEQcCePXvEHrO4uBgVFRWS13Lw4EFxIG21WnHxxRdLtjt58iTOnDkDIPaQTZ06VbJTaGtr69VD5vf7Ex6ycePGiWECStCvPfHw4cOxbt06LF26FLW1tZgwYQJuv/12mM1m/OUvf4Hf78e0adNw1113pX2s3tSzoxw7dkyMO+7uJrW0tCT4rqdOnSo5SPX5fKKNCcR6QqmbxPM8du/endATDhs2LKkdIQTV1dWine90OiXHDABw6NAh0ftiNBoxZcoUyYfs9OnTOHHiBIBYTyj3kHV2duLAgQPi75MmTRLdl/GEQiHRTARib7Li4uKkdl0fspKSEowaNUryWlKl36edb7vtNuTl5WH58uVoaWnBs88+K26bOnUqNm7cqIhx39sRbkdHh/gqp+cgdZMCgQCqq6vF38eNGyfZE9JXNHVNDR06VPImEULw/fffi6N1h8OBiRMnSp5jfX292PsbDAZMnz5d8iFrampCfX09gNj1z5gxQ/Ih83q9Cb3/xIkTJR+ySCSCPXv2iL75UaNGdfuQhUKxfB5FRUWyb7KamhoxM5PFYpF9yNJhQOKJb7rpJhw7dgzvvPMOvvrqK0QiEcyZMwd33HGHYktUepOwWeoVTYPpu+6ra+8xevRoyX0eOHAg4RVdWVkpeZNOnDgh+kx1Oh2mTZsmaRN2fUVPnTpVchDZ9RU9fvx45OfnJ7Xr2vuXl5dL+u0JIfjuu+8STKzx48dLXnNdXZ1oGhqNRtk3WUNDQ8Jgc/r06Yrd73gGLCje4XDgnnvuwT333NMv+48XsZQ4qB1ME6kUFBTIvqLr6urQ0RFb1mQymWR7j4aGBpw8eRLA+ZskZwcfPHhQ/H3y5MmSvX84HMaePXvE3y+66KJuX9HU9i8tLZV9RR84cECMJ7HZbLK9//Hjx9HU1AQgNr6YNm2apDBbW1tx+PBhAOdtf6ne3+/34/vvvxd/nzBhgqJ2cDyZH1nSS+IHc1JT3PE3yWAwdHuTjhw5AiB2k+R6j643adKkSbI+1K6vaDmPwN69e3v1iq6rq+vVK7qhoQGnTsWClLRabbcPWbyJNWXKFNnBZryJNXbsWBQWFia1EwQB1dXVCSbWyJEjJa9FCVQjYnrzdRJhh52dnUk9oZRvNBgMJvSEY8eOlXxFC4KA3bt3J9wkuRnHffv2iRMwdrsd48aNk2x39OjRhIdMTpjNzc0JD9nUqVMlH1qv15sw2Jw4caI4qxlPJBLBt99+m/CQSZlYUnawnIlVW1srvsksFousb1spVCfirq82+oqmdvAFF1wgawd/9913Ca4quZsU76rq7iZ1dVXJ2cHt7e1JriqphywQCCS5qqQesq7+4GHDhkk+ZF39wQ6HQ/YhO3LkSIKrUe4ha2pqSnA1yj1kSqIaEdNRf7ytKeUPlrtJXQcr3b2ijx8/DuC8uSF1k9xud8LM2KRJkyR7QqnBJk0/EA/t/eMHmxdccIHktezfv1+MW7FarbIP2fHjxxMGm9OnT5d8yFpbW3Ho0CHx96lTp0o+ZH6/P8HckBtsKo1qREx9qvEj+RMnTiS8ouVcVVJ2cG8GK3KuKtoT0lf0yJEjZV1V8XZwYWGh7GDz8OHDoh1sNpsH1Q6Wm3SJt4PLysoU9wfLoRoRdzUnIpFIkh0s5arieT7hFT127FjJeAcquHg7WC7EtLa2NsEOnjBhgmS706dPJzxkcq6qjo6OJI+A1GAzFAol+YPl4h2qq6sT7OAhQ4ZInuP+/fsT7OALL7xQst3Ro0cT7GA5V2N/oFoR63Q6jBkzBhzHoaKiQtIOpu0mTJgAnU6HwsJCWTuY4ziMHTsWZrMZZrO528FKRUUFnE6nODMmFyNQWloqeirkBpsAxMg1ABgzZozsK9pgMGDcuHHQaDQoKyuTHWxyHIfx48fDaDTCZrPJmljA+ZgSvV7f7URFeXk5ioqKBswOToCohI0bNxIAZNasWQmfd3Z2Ep7ne/y+3+8nfr+/x3bhcJi43e4e20WjUdLR0dFjO0EQSHt7e4/tCCGkvb2dCILQYzuPx0NCoVCP7YLBIPF6vT2243medHZ29tiuL9eiJAOaAag/kRrYAei1g13K1JCCruroCY1GI4aidgfHcb0e/PS2ndQAUgoa49wTWq22V3/HvlyLkqjGnKBBMIPxR2QMLqoRMV09IjWDxFA3qhEx9Yv21/w8I3NRjYipD1XKpcRQN6oTMbOJcw/ViJhOLkitk2OoG9WIWGramZEbqEbENKqM2cS5h2pETFcvMBHnHqoRMZ3s6M0sGUNdqELEhBDRnJDLuMNQL6oQcTgcFoPK6cCOEIJwONyrVdCMgaG/7okqRExXMAPnRdzU1ASj0ShmitTr9d0mNmT0P/11T1QRxWa328XUrXKJBXmeTztVFiM9+uueqOKuchyXlCaULtCk6HQ6yVwPjIGjv+6JKswJKeg0NKW0tJT1xINMf90T1Yq4a1FG1gsPPv11T1Qr4q5p9lmI5uDTX/dEtSLuOohgpsTg01/3RLUijne7AZBdScwYOPrrnqhWxF3tLxaiOfj01z1RrYhpaCZFKrsNY2Dpr3uiWhEDiXmKmU2cGfTHPenX2s6ZADlXOZMQwoScISh9T1QvYob6UbU5wcgNmIgZWQ8TMSPrUe1I5+uvv8bWrVtx+vRpEEIwceJE3HLLLf1aAIWRSHV1NTZs2IAHHnggacVNTU0NNm7ciGPHjiEQCGD06NFYuHAhJk+e3PcDDXgezn5GEATyk5/8hABI+tHpdOTNN98c7FPMGW644QYCgDzzzDMJnz/33HNEo9FI3qNVq1b1+TiqE/Ef/vAH8Q9y0UUXkV/84hfkpptuIg6HgwAgRqORHDhwYLBPMyeYOHEiAUDuvvtu8bMtW7YQjuMIAFJUVETuu+8+cuutt5KhQ4eK923Dhg19Oo6qROxyuYjNZiMAyNVXX03C4bC4ze/3k6uuuooAIDfddNMgnmXuUFlZSQCQZcuWEUJib8nJkycTAGTkyJGkublZbBuNRsmKFSsIAFJRUUEikUivj6MqEf/lL38hAIjFYiEnT55M2r5t2zbxaW9oaBiEM8wtrrjiCgKA/PjHPyaEEPLtt9+Kf//PP/88qX1zczPRarUEAPnss896fRxVeSfeeustAMDy5csxfPjwpO1XXnmlGIgdX9KK0T/Q2AiajZ7en0svvRRz585Nal9cXIyZM2cCQEJlp55QlYhpqaoFCxZIbtfpdKJ3giYgZPQftH4gjVajVark7g8AscBOX+6PakTc1tYmZgGaPn16j+3lKhoxlKNrut26ujoAyt8f1YiYrqTlOE6yWCCFlUUYGAghOHbsGICYiAkhYvVSqYqplFTuj2pETMvFEkLEgoldOXPmjFj8kE169C+NjY3imrrCwkJEo1Gx1rTc/REEQSwQ35f7oxoRx+clpkLtyqZNmwDECnZ31xsw0ie+4HphYSG0Wq2Y2Obs2bOS36murkZzczM4jsOUKVN6fSzViDhelJ999plkm9dffx0AsHDhwgEr2ZqrUPsXiK1q5jhOvEdy96eqqgoAMHv2bNkyvVKoRsTFxcViEfDnn38+aVHi3//+d3z11VcAgP/4j/8Y8PPLNeJ7W/qWnDp1KoCYWLv2xg0NDXjttdcAACtWrOjbwZRxa2cGf/7zn0Vn+oIFC0hTUxMhhJCTJ0+SkSNHEgBk8eLFg3yWucHPfvYz8V4cOXKEEELI119/LU5mVFZWktraWkIIIV6vV5xNHT9+fJ9m6whR2YydIAjkwQcfFP94RqORjB07lhiNRgKADBs2jDQ2Ng72aeYEjz76qHgfWlpaxM/Xr18vCpnjODJmzBgxrsVkMpFdu3b1+ViqEjHljTfeIKNGjUqIjrr66qvJsWPHBvvUcoavv/6aACCzZ89O2rZ9+3YyY8aMhPszZswYsmPHjpSOpdo1doIgYPfu3WhoaMDIkSMxefJkNpgbYE6ePImSkhLZJCk1NTU4cuQInE4nLrvsspQnoFQrYkbuoBrvBCN3YSJmZD1MxIysh4mYkfUwETOyHiZiRtbDRMzIepiIGVkPEzEj62EiZmQ9TMSMrIeJmJH1MBEzsp7/D5y57S6M8dKoAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.xkcd():\n",
    "\n",
    "    # Determine the spacing\n",
    "    total_width  = 10 # m\n",
    "    total_height = vegetation_height + soil_depth + 1 # [m], +1 to have some space for bedrock\n",
    "    slope_l = soil_depth + 0.5*(topography_slope * total_width)\n",
    "    slope_r = soil_depth - 0.5*(topography_slope * total_width)\n",
    "    \n",
    "    # General figure components\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_axes((0.1, 0.2, 0.8, 0.7))\n",
    "    ax.set_aspect('equal') \n",
    "    ax.spines[['top', 'right', 'bottom']].set_visible(False)\n",
    "    #ax.set_xticks([])\n",
    "    #ax.set_yticks([])\n",
    "    plt.ylim([-total_height*0.1, total_height])\n",
    "    \n",
    "    # Bedrock\n",
    "    n_bedrock_shading = 10\n",
    "    for x in np.linspace(0,total_width,n_bedrock_shading, endpoint=False):\n",
    "        plt.plot([x,x+total_width/n_bedrock_shading],\n",
    "                 [-total_height*0.05,0], \n",
    "                 color=get_pm_plot_color('bedrock fracture'))\n",
    "    plt.plot([0,total_width],[0,0], color=get_pm_plot_color('bedrock boundary'))\n",
    "    \n",
    "    # Surface\n",
    "    #plt.plot([0,total_width],[slope_l,slope_r], color=get_pm_plot_color('surface'))\n",
    "    \n",
    "    # Soil\n",
    "    for ix,(soil_t,soil_b) in enumerate(zip(soil_layer_top,soil_layer_bot)): # SOILGRIDS layer loop\n",
    "        if soil_depth >= soil_t: # the top of this soil layer is within estimated soil depth\n",
    "    \n",
    "            # Get the soil type in this layer\n",
    "            _,soil_name = find_usda_soilclass(soil_sand[ix]*100,soil_silt[ix]*100,soil_clay[ix]*100)\n",
    "            soil_color = get_pm_soil_color(soil_name)\n",
    "            \n",
    "            # Check where the bottom of the plotting area is\n",
    "            if soil_depth < soil_b: # Estimated depth ends before this soilgrids layer does\n",
    "                soil_b = soil_depth # Therefore end the plotting layer at estimated soil depth\n",
    "    \n",
    "            # Handle the special case where estimated soil depth is deeper than deepest soil grid layer\n",
    "            if soil_b == 2: \n",
    "                soil_b = soil_depth # extend this layer all the way to bedrock\n",
    "    \n",
    "            # Calculate at which fraction of total soil depth the layer interfaces are\n",
    "            soil_t_frac = soil_t / soil_depth\n",
    "            soil_b_frac = soil_b / soil_depth\n",
    "    \n",
    "            # Create a Polygon patch with the specified characteristics\n",
    "            x = [total_width, total_width, 0, 0, total_width]\n",
    "            y = [(1-soil_b_frac) * slope_r, (1-soil_t_frac) * slope_r, \n",
    "                 (1-soil_t_frac) * slope_l, (1-soil_b_frac) * slope_l,\n",
    "                 (1-soil_b_frac) * slope_r]\n",
    "            \n",
    "            soil_patch = patches.Polygon(\n",
    "                xy=list(zip(x, y)),  # Combine x and y coordinates\n",
    "                closed=True,  # Connect the last point to the first to close the shape\n",
    "                facecolor=soil_color,\n",
    "                edgecolor='none'#,\n",
    "                #linestyle='dotted',\n",
    "                #linewidth=1.0  # Width of the outline\n",
    "            )\n",
    "    \n",
    "            ax.add_patch(soil_patch)\n",
    "    \n",
    "    # Vegetation\n",
    "    \n",
    "    \n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ce3007e-b514-49a9-abce-538bff0184af",
   "metadata": {},
   "source": [
    "### Plotting functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3c832568-717d-4fa0-a167-8fad4b892b0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_usda_soilclass(sand,silt,clay,no_data_value=255):\n",
    "\n",
    "    # Adapted from: Knoben (2021)\n",
    "    # Knoben, W. J. M. (2021). Global USDA-NRCS soil texture class map, \n",
    "    # HydroShare, https://doi.org/10.4211/hs.1361509511e44adfba814f6950c6e742\n",
    "    \n",
    "    # Based on Benham et al., 2009 and matching the following soil class table:\n",
    "    # SUMMA-ROSETTA-STAS-RUC soil parameter table\n",
    "    # 1  'CLAY' \n",
    "    # 2  'CLAY LOAM'\n",
    "    # 3  'LOAM' \n",
    "    # 4  'LOAMY SAND'\n",
    "    # 5  'SAND'\n",
    "    # 6  'SANDY CLAY'\n",
    "    # 7  'SANDY CLAY LOAM'\n",
    "    # 8  'SANDY LOAM'\n",
    "    # 9  'SILT'\n",
    "    # 10 'SILTY CLAY'\n",
    "    # 11 'SILTY CLAY LOAM'\n",
    "    # 12 'SILT LOAM'\n",
    "\n",
    "    # Check inputs\n",
    "    if (sand<=1) & (silt<=1) & (clay<=1):\n",
    "        print(f'WARNING: find_usda_soilclass(): sand, silt and clay all < 1. Please provide inputs in % [0,100].')\n",
    "        return\n",
    "    \n",
    "    # Legend\n",
    "    soiltype = ['clay','clay loam','loam','loamy sand','sand','sandy clay','sandy clay loam',\n",
    "                'sandy loam','silt','silty clay','silty clay loam','silt loam']\n",
    "    \n",
    "    # Classify\n",
    "    if (clay >= 40) & (sand <= 45) & (silt < 40):\n",
    "        soilclass = 1\n",
    "    elif (clay >= 27) & (clay < 40) & (sand > 20) & (sand <= 45):\n",
    "        soilclass = 2\n",
    "    elif (clay >= 7) & (clay < 27) & (silt >= 28) & (silt < 50) & (sand < 52):\n",
    "        soilclass = 3\n",
    "    elif ((silt + 1.5 * clay) >= 15) & ((silt + 2* clay) < 30):\n",
    "        soilclass = 4\n",
    "    elif ((silt + 1.5 * clay) < 15):\n",
    "        soilclass = 5\n",
    "    elif (clay >= 35 )& (sand > 45):\n",
    "        soilclass = 6\n",
    "    elif (clay >= 20) & (clay < 35) & (silt < 28) & (sand > 45):\n",
    "        soilclass = 7\n",
    "    elif ((clay >= 7) & (clay < 20) & (sand > 52) & ((silt+2*clay) >= 30)) | ((clay < 7) & (silt < 50) & ((silt+2*clay >= 30))):\n",
    "        soilclass = 8\n",
    "    elif (silt >= 80) & (clay < 12):\n",
    "        soilclass = 9\n",
    "    elif (clay >= 40) & (silt >= 40):\n",
    "        soilclass = 10\n",
    "    elif (clay >= 27) & (clay < 40) & (sand <= 20):\n",
    "        soilclass = 11\n",
    "    elif ((silt >= 50) & (clay >= 12) & (clay < 27)) | ((silt >= 50) & (silt < 80) & (clay < 12)):\n",
    "        soilclass = 12\n",
    "    elif (sand == no_data_value) | (silt == no_data_value) | (clay == no_data_value):\n",
    "        soilclass = 0\n",
    "\n",
    "    soilname = 'undefined'\n",
    "    if soilclass > 0:\n",
    "        soilname = soiltype[soilclass-1] # -1 to match 0-based indexing\n",
    "    \n",
    "    return soilclass, soilname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f6a41a12-759d-467c-a75e-78ed17e08486",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pm_soil_color(soil_name):\n",
    "\n",
    "    '''Creates a soil color for a give sand/silt/clay percentage'''\n",
    "\n",
    "    # Color source: https://kremesti.com/water/silt_clay_mud.htm (accessed 2023-11-27)\n",
    "    if soil_name == 'clay':\n",
    "        color = [97,62,38]\n",
    "    elif soil_name == 'clay loam':\n",
    "        color = [68,45,30]\n",
    "    elif soil_name == 'loam':\n",
    "        color = [71,66,49]\n",
    "    elif soil_name == 'loamy sand':\n",
    "        color = [66,61,41]\n",
    "    elif soil_name == 'sand':\n",
    "        color = [91,87,67]\n",
    "    elif soil_name == 'sandy clay':\n",
    "        color = [85,82,71]\n",
    "    elif soil_name == 'sandy clay loam':\n",
    "        color = [78,64,40]\n",
    "    elif soil_name == 'sandy loam':\n",
    "        color = [74,71,42]\n",
    "    elif soil_name == 'silt':\n",
    "        color = [48,48,48]\n",
    "    elif soil_name == 'silty clay':\n",
    "        color = [73,52,39]\n",
    "    elif soil_name == 'silty clay loam':\n",
    "        color = [61,52,47]\n",
    "    elif soil_name == 'silt loam':\n",
    "        color = [88,83,77]\n",
    "    else:\n",
    "        print(f'WARNING: get_pm_soil_color(): soil {soil_name} not defined.')\n",
    "        color = None\n",
    "    return [C/255 for C in color]    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac5fdb35-fbf9-4be0-8e93-1002f6ca2c63",
   "metadata": {},
   "source": [
    "### Data extraction functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b2b4c049-27be-4952-b266-9396969b9c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_geotiff_properties(file):\n",
    "\n",
    "    with rasterio.open(file) as src:\n",
    "        geoTransform = src.transform\n",
    "        pixelSize = [geoTransform[0],geoTransform[4]]\n",
    "        rasterSize = [src.width,src.height]\n",
    "        boundingBox = src.bounds\n",
    "    \n",
    "    return pixelSize, rasterSize, boundingBox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e23061aa-483e-473c-98e8-f14311e84bed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_geotif_data_as_array(file, band=1):\n",
    "    ds = gdal.Open(file) # open the file\n",
    "    band = ds.GetRasterBand(band) # get the data band\n",
    "    data = band.ReadAsArray() # convert to numpy array for further manipulation   \n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c56fc5f0-d257-409f-a4c9-7fa686af37f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_geotif_noData(src_file, band=1):\n",
    "    src_ds = gdal.Open(src_file)\n",
    "    src_band = src_ds.GetRasterBand(band)\n",
    "    no_data_value = src_band.GetNoDataValue()\n",
    "    src_ds = None\n",
    "    return no_data_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1d7d6345-5dfd-4510-9bd4-c2850eefb01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def confirm_dimensions_match(files):\n",
    "    pixel_sizes = []\n",
    "    raster_sizes = []\n",
    "    bounding_boxes = []\n",
    "    \n",
    "    for sub_list in files:\n",
    "        file = sub_list[1]\n",
    "        pixel, raster, bb = get_geotiff_properties(file)\n",
    "        pixel_sizes.append(pixel)\n",
    "        raster_sizes.append(raster)\n",
    "        bounding_boxes.append(bb)\n",
    "\n",
    "    assert all([sub_list == pixel_sizes[0] for sub_list in pixel_sizes]), 'pixel sizes don''t match'\n",
    "    assert all([sub_list == raster_sizes[0] for sub_list in raster_sizes]), 'raster sizes don''t match'\n",
    "    assert all([sub_list == bounding_boxes[0] for sub_list in bounding_boxes]), 'bounding boxes don''t match'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "8c78434f-0fa3-4cc8-a727-4a0fc6ed203d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_pixel_value(tiff_path, pixel_x, pixel_y):\n",
    "\n",
    "    file = str(tiff_path)\n",
    "    \n",
    "    # Get data\n",
    "    with rasterio.open(file) as src:\n",
    "        pixel_value = src.read(1, window=((pixel_y-1, pixel_y), (pixel_x-1, pixel_x)))\n",
    "        pixel_value = pixel_value[0][0]\n",
    "\n",
    "    # Get scaling and offset\n",
    "    src = gdal.Open(file)\n",
    "    scale_factor = src.GetRasterBand(1).GetScale()\n",
    "    offset = src.GetRasterBand(1).GetOffset()\n",
    "    if scale_factor: pixel_value = pixel_value * scale_factor\n",
    "    if offset: pixel_value = pixel_value + offset\n",
    "    src = None\n",
    "    \n",
    "    return pixel_value # Extract the actual value from the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "d3f0d2e6-dc21-4360-8c58-5ca3e64b4d42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main: called from get_location_data_dict_from_geotiffs()\n",
    "def convert_lat_lon_to_raster_x_y(file,lat,lon):\n",
    "    pixel, raster, bb = get_geotiff_properties(file)\n",
    "    assert is_point_inside_bbox((lon,lat),bb), f'({lon},{lat}) not in {bb}' # sub 1\n",
    "    lon_lat_to_pixel(lon,lat, pixel[0], raster, bb) # sub 2\n",
    "    return x,y\n",
    "\n",
    "# Sub 1\n",
    "def is_point_inside_bbox(point, bbox):\n",
    "    x, y = point\n",
    "    min_x, min_y, max_x, max_y = bbox.left, bbox.bottom, bbox.right, bbox.top\n",
    "    return min_x <= x <= max_x and min_y <= y <= max_y\n",
    "\n",
    "# Sub 2\n",
    "def lon_lat_to_pixel(lon, lat, pixel_size, raster_size, bounding_box):\n",
    "    # Calculate pixel coordinates\n",
    "    x = int((lon - bounding_box.left) / pixel_size)\n",
    "    y = int((bounding_box.top - lat) / pixel_size)\n",
    "\n",
    "    # Ensure the coordinates are within the raster bounds\n",
    "    x = max(0, min(x, raster_size[0] - 1))\n",
    "    y = max(0, min(y, raster_size[1] - 1))\n",
    "\n",
    "    return x,y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "a1b1c1e2-9c0e-4b1a-bb1a-6b3ce6a7b3ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_location_data_dict_from_geotiffs(files,lon,lat):\n",
    "\n",
    "    # Confirm the all rasters have the same dimension (function will error out if not)\n",
    "    confirm_dimensions_match(files)\n",
    "    \n",
    "    # Get the raster indices for the provided lat/lon coordinates\n",
    "    x,y = convert_lat_lon_to_raster_x_y(files[0][1],lat,lon)\n",
    "\n",
    "    # Get the data values we want\n",
    "    data_dict = {}\n",
    "    for sublist in files:\n",
    "        data_name,data_file = sublist\n",
    "        data_value = extract_pixel_value(data_file,x,y)\n",
    "        data_dict[data_name] = data_value    \n",
    "    \n",
    "    return data_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7407e158-c313-41f3-b456-f986e9557716",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "perc-model-env",
   "language": "python",
   "name": "perc-model-env"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}