{ "cells": [ { "cell_type": "code", "execution_count": 64, "id": "22b944c0-d8e0-478c-a3d8-ea670d90cc2f", "metadata": {}, "outputs": [], "source": [ "import yaml" ] }, { "cell_type": "code", "execution_count": 65, "id": "2607ecaa-758c-4a3e-80f2-050f2c73df8b", "metadata": {}, "outputs": [], "source": [ "# load menu\n", "with open(\"mnt/city-directories/01-user-input/menu.yml\", 'r') as f:\n", " menu = yaml.safe_load(f)" ] }, { "cell_type": "code", "execution_count": 66, "id": "a6992e26-b62f-41fc-85e7-66c95e7380e9", "metadata": {}, "outputs": [], "source": [ "if menu['all_stats']:\n", " import os\n", " import glob\n", " import math\n", " import geopandas as gpd\n", " import pandas as pd\n", " import numpy as np\n", " from io import StringIO\n", " import requests\n", " from sklearn.preprocessing import MinMaxScaler\n", " from shapely.geometry import shape\n", " from shapely.ops import unary_union\n", " import pint\n", " import folium\n", " from pathlib import Path\n", " import matplotlib.pyplot as plt\n", " import requests\n", " import re\n", " import rasterio\n", " from rasterio.mask import mask\n", " from shapely.geometry import Point\n", " from fiona.crs import from_epsg\n", " from nbconvert import MarkdownExporter\n", " import nbformat\n", " import base64\n", " import pickle\n", " import plotly.graph_objects as go\n", " import seaborn as sns\n", " import plotly.express as px" ] }, { "cell_type": "code", "execution_count": 67, "id": "0ca8d4ff-ae23-4850-95f6-4a374d77ab88", "metadata": {}, "outputs": [], "source": [ "url = \"https://raw.githubusercontent.com/compoundrisk/monitor/databricks/src/country-groups.csv\"\n", "country_groups = pd.read_csv(url)\n", "\n", "# Source helper functions\n", "helpers_url = \"https://raw.githubusercontent.com/compoundrisk/monitor/databricks/src/fns/helpers.R\"\n", "helpers_code = requests.get(helpers_url).text\n", "\n", "# Define tolatin function\n", "def tolatin(x):\n", " return stri_trans_general(x, id=\"Latin-ASCII\")\n", "\n", "# Define normalize function\n", "def normalize(x):\n", " x_min = np.min(x)\n", " x_max = np.max(x)\n", " return (x - x_min) / (x_max - x_min)\n", "\n", "def print_text(x, linebreaks=2):\n", " print(x + \"\\n\" + \"
\" * linebreaks)" ] }, { "cell_type": "code", "execution_count": 68, "id": "a5df3a99-85dd-4b6b-93d7-bbd0a68679f7", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAHFCAYAAABW0zbJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSCElEQVR4nO3deVxU5eIG8Gdm2GSZURSUXUUEEUEUF1RMDUlcgux2bVHUNNOsNH/dTNO6pobeytJSs7I0Kylz61ZStuBWmigYLimuoA7gxr4P5/cH16lRUYQZ3jMzz/fzmU85y5nnjDCP55z3vEchSZIEIiIiGVGKDkBERHQjlhMREckOy4mIiGSH5URERLLDciIiItlhORERkeywnIiISHZYTkREJDssJyIikh2WE5EMLVu2DAqFAiEhIXU+59ixYxg3bhx8fX1hZ2eHVq1aYejQodi2bdtNz01JSYFCocBXX31lythERsNyIpKhjz76CABw5MgR7Nu376bHN23ahPDwcPz++++YO3cufvzxR6xcuRIAMHToULzwwgtNmpfI2GxEByAiQ6mpqTh06BCGDRuGb7/9FqtXr0avXr30j586dQpjxoxBly5dkJKSAicnJ/1jDz30EKZMmYLXX38d3bp1w8MPPyxiFYgajVtORDKzevVqAMCiRYvQp08fJCUlobS0VP/4W2+9hdLSUrzzzjsGxXTdm2++iebNm2PhwoVNlpnI2FhORDJSVlaG9evXo0ePHggJCcHjjz+OoqIibNiwQf+c7du3o3Xr1ujdu/ctl+Ho6IiYmBgcPnwYOTk5TRWdyKhYTkQy8tVXX6GgoAATJkwAAIwaNQrOzs76rSkAyMrKQrt27W67nOuPZ2VlmS4skQmxnIhkZPXq1WjWrJn+WJGzszMeeugh7Nq1C5mZmfVezvXLtCkUCpPkJDI1lhORTJw8eRI7d+7EsGHDIEkS8vPzkZ+fj3/84x8A/hrB5+vrizNnztx2WWfPngUA+Pj4mDQzkamwnIhk4qOPPoIkSfjqq6/QokUL/W3YsGEAgLVr10Kn02Hw4MHIzc3F3r17b7mc0tJSbN++HSEhIWjTpk1TrgKR0XAoOZEM6HQ6rF27Fv7+/vjwww9vevybb77Bm2++iW3btuG5557DRx99hGeeeeamoeQA8Pzzz+PatWv6856IzBHLiUgGtm3bhosXL2Lx4sUYMGDATY+HhITg3XffxerVq7F582asW7cOjz32GHr06IEZM2YgMDAQubm5+Oijj7Bt2zY8//zzGDVqVNOvCJGRsJyIZGD16tWws7PD+PHjb/l4q1at8MADD+Crr75Cbm4uHnzwQXTq1An/+c9/MG/ePOTm5sLFxQU9e/bEt99+i6FDhzbxGhAZl0K6PqyHiIhIJjgggoiIZIflREREssNyIiIi2WE5ERGR7LCciIhIdlhOREQkOzzPqYFqampw8eJFuLi4cHJNIrJakiShqKgInp6eUCqNt73DcmqgixcvclJNIqL/yc7Ohre3t9GWx3JqIBcXFwC1fyFqtVpwGiIiMQoLC+Hj46P/TjQWllMDXd+Vp1arWU5EZPWMfXiDAyKIiEh2WE5ERCQ7LCciIpIdlhMREckOy4mIiGSH5URERLLDciIiItlhORERkeywnIiISHZYTkREJDssJyIikh2WExERyQ7LiYiIZIflRHQXqnQ1KCyvEh2DyOLxkhlkkap0NSit0KG0qholFTqUVlajtPJv/63QoeRv95VU6FBWaXhfaaUOJRXV/7u/9r4qnQQA6NnWFU8P6oCogFa8EjKRCQgtp5UrV2LlypU4e/YsAKBz5854+eWXERsbCwDIzc3FzJkz8cMPPyA/Px/9+/fHO++8g4CAgDqXuWbNGowfP/6m+8vKyuDg4HDT/YmJiZg9ezamTZuGt99+2yjrRQ1XpavBxfwyZF8tQ35ZZW3BVFbry+HvxWJYJjqUVtQ+r6xSh0pdjUlz/n72KhI++h2h3hpMHdgBgzu1hlLJkiIyFqHl5O3tjUWLFqFDhw4AgLVr1yIuLg5paWkIDg5GfHw8bG1tsXXrVqjVaixZsgTR0dE4evQonJyc6lyuWq3G8ePHDe67VTHt378f77//PkJDQ427YlQnSZJwqbgC2VfLkH21FNlXS5F1tRTZ10qRfbUM2oIy1EiiU9bfH+cL8OS6Awhs7YKnBvpjeKgnVCwpokYTWk4jRoww+PPChQuxcuVK7N27F7a2tti7dy8OHz6Mzp07AwBWrFgBd3d3rF+/HhMnTqxzuQqFAm3atLntexcXF+Oxxx7DBx98gAULFjR+ZUivuKLaoHjOXyurLaD/lVB5lWm3akQ4nluEaUnpeGv7CTw1oAPiw71gZ8NDukQNJZtjTjqdDhs2bEBJSQkiIyNRUVEBwHCLR6VSwc7ODrt3775tORUXF8PPzw86nQ5du3bF/PnzER4ebvCcqVOnYtiwYYiOjq5XOVVUVOgzAUBhYeHdrqLF+Puut+xrpX8Vz9VSZF8rw9WSStERhTl7pRQvbPwDb/94ApMH+OOfET5wsFWJjkVkdoSXU0ZGBiIjI1FeXg5nZ2ds3rwZwcHBqKqqgp+fH2bNmoVVq1bByckJS5YsQU5ODrRabZ3LCwoKwpo1a9ClSxcUFhZi6dKl6Nu3Lw4dOqQ/VpWUlISDBw9i//799c6ZmJiIefPmNXp9zYEkSbhcXPm/rZ5SZF35a7db1tVSs9v1JsLFgnK8vPUIlv10Ek9EtcNjvf3gbC/8143IbCgkSRL6NVNZWYmsrCzk5+dj48aN+PDDD7Fjxw4EBwfjwIEDmDBhAg4dOgSVSoXo6GgolbW7Sr777rt6Lb+mpgbdunVD//79sWzZMmRnZyMiIgI//PADwsLCAAADBgxA165dbzsg4lZbTj4+PigoKIBarW74ByBISUV17VbPldqtnRt3w5VV6URHtCiaZrZ4vG87jOvTFhpHW9FxiIymsLAQGo3G6N+FwsvpRtHR0fD398eqVav09xUUFKCyshJubm7o1asXIiIisHz58nov84knnsD58+exbds2bNmyBQ888ABUqr92teh0OigUCiiVSlRUVBg8VhdT/YUY0/WyyTIYeFBbRNa8600kZ3sbjIn0w4R+7dDK2V50HKJGM9V3oez2M0iSZLCFAgAajQYAkJmZidTUVMyfP/+ulpeeno4uXboAAO69915kZGQYPGf8+PEICgrCzJkz61VMcnf4QgEWbfsTu09eFh2FblBcUY2VKafw8Z4zeKSnLyb1bw8PTTPRsYhkR2g5zZ49G7GxsfDx8UFRURGSkpKQkpKC5ORkAMCGDRvg5uYGX19fZGRkYNq0aYiPj0dMTIx+GQkJCfDy8kJiYiIAYN68eejduzcCAgJQWFiIZcuWIT09Xb+l5eLigpCQEIMcTk5OaNmy5U33m5vsq6V444fj2Jp+UXQUuoPyqhp8vOcsPt17Dv/o7o3J9/jDr2Xdp0cQWRuh5ZSbm4sxY8ZAq9VCo9EgNDQUycnJGDx4MABAq9VixowZyM3NhYeHBxISEjB37lyDZWRlZemPQwFAfn4+Jk2ahJycHGg0GoSHh2Pnzp3o2bNnk65bU7pWUol3fj6JT/eeM/nJp2RcVToJ63/Pxhf7sxHX1QtPDfBHQGsX0bGIhJPdMSdzIYdjTmWVOny05wzeSzmFoopqIRnIuBQKYEjnNpg6sANCvDSi4xDdkdUcc6I709VI+OpANt7anomcwnLRcciIJAnYdjgH2w7nYECgG2YM7ohQ7+aiYxE1OZaTGZEkCT8dy8Pi5D+RmVcsOg6ZWMrxS9h3+iqSJvVGmE9z0XGImhTnVzETaVnXMGrVXkz8JJXFZEXKqnSYsHY/sq6Uio5C1KRYTjJ3+lIxpnx6AA+s+BW/n70qOg4JcLm4EuM+/h3XeG4aWRGWk0xdKqrAnC0ZGPzWTmw7nCM6Dgl2+nIJJn6SinLO3EFWguUkM8UV1Xhr+wnc8/ov+HRvFnScxI7+58C5a3jui3T+TJBVYDnJRJWuBut+O4sBr/+CpT9lorSS/0Kmm207nIOF3x4THYPI5DhaTzBJkrDtcA5e//44zlwuER2HzMBHe87Aq0UzTOjXTnQUIpNhOQm07/QVJG77E+nZ+aKjkJlZ8O1ReGocENvFQ3QUIpNgOQlQXqXD058fxI/H8kRHITMlScC0L9Lh5mKPiLauouMQGR2POQngYKuCrysn+aTGqayuwcRPUnHqEs97I8vDchJk1tAgdPdrIToGmbn80iqM+/h3XCqquPOTicwIy0kQW5USyx/thpZOdqKjkJnLvlqGCWv3o7SSk/+S5WA5CdRG44Blj4RDqRCdhMzdH+cL8MznaajmJVPIQrCcBOvboRVmDO4oOgZZgJ/+zMMrXx8Br4JDloDlJANPDeiAQUHuomOQBfhsXxZW7jglOgZRo7GcZECpVGDJP8Pg3aKZ6ChkAf6TfBxb0y+IjkHUKCwnmWjuaIeVj3WHnYp/JdR4z284hF9PXRYdg6jB+E0oI128Nfj3/Z1FxyALUKWT8OS6AzieUyQ6ClGDsJxk5pGePhjZzUt0DLIAReXVGP/x78gtLBcdheiusZxkRqFQYGF8FwS1cREdhSzAxYJyjPt4P4rKq0RHIborLCcZamanworHusHZnlMfUuMd0xbiqc8OoornQJEZYTnJVHs3Z7zxUKjoGGQhdmVexqxNGTwHiswGy0nGhoR4YCKv2UNG8tWB83j7x0zRMYjqheUkczNjg9CjLSeIJeNY+lMmvtyfLToG0R2xnGTOVqXEu492Qytne9FRyELM2pyBHScuiY5BdFssJzPQWu2AdzhBLBmJrkbCU58ewOELBaKjENWJ5WQmIv1b4vn7AkXHIAtRUqnD42v240J+megoRLfEcjIjk/v7I7oTJ4gl48grqsC4j35HQSnPgSL5YTmZEaVSgTcf6gofV04QS8aRmVeMSetSUVGtEx2FyADLycxoHG1rJ4i14V8dGce+M1fxrw1/oKaG50CRfPAbzgyFeGnwKieIJSP6+tBF/Of746JjEOmxnMzUqB4+eKi7t+gYZEHe23EK6/aeEx2DCADLyWwpFArMjw9BJw+16ChkQV7Zehjbj+aKjkHEcjJnDrYqrHysG1w4QSwZSY0EPLP+INKz80VHISvHcjJzbVs5YfPUPvhHd2/YqniWLjVeeVUNJqzZj3NXSkRHISvGcrIAHdxd8MZDYdg9cxCmDPCHiwO3pKhxrpRUYtzH+3G1pFJ0FLJSLCcL0lrtgJlDgvDbrHvx8vBgeDXn+VDUcGcul2Di2v0or+I5UNT0WE4WyNneBo/3a4cd/xqAdx4JRxcvjehIZKYOZuVjelI6dDwHipoYy8mC2aiUGBHmia+f7oukSb1xbxCnPqK7l3wkBwu+PSo6BlkZHpywAgqFAr3bt0Tv9i1xMq8IH+46g00HL6CSl+2mevp4z1l4NW+GiVHtRUchK6GQeN3mBiksLIRGo0FBQQHUavM71yivqByf/HoO6/aeQ0EZJ/6kO1MogOWPdsPQLh6io5CMmOq7kOXUQOZeTteVVlZjQ+p5fLj7NLKv8vIJdHt2Nkp8NrEXerR1FR2FZMJU34U85mTlHO1sMLZPW6Q8PxArHuuGMJ/moiORjFVW1+CJT1Jx6lKx6Chk4VhOBABQKRUY2sUDW57qgw2TIzE4uDUUPKeXbiG/tArjPv4dl4oqREchC8ZyIgMKhQI92rrig4QI/DjjHjzay5eX56CbZF8tw4S1+1FaWS06Clkood86K1euRGhoKNRqNdRqNSIjI7Ft2zb947m5uRg3bhw8PT3h6OiIIUOGIDMz87bLXLNmDRQKxU238vLyer8v1fJ3c8ZrD3TBry8OwrR7A9DC0VZ0JJKRP84X4OnP01DNUZ9kAkLLydvbG4sWLUJqaipSU1MxaNAgxMXF4ciRI5AkCfHx8Th9+jS2bt2KtLQ0+Pn5ITo6GiUlt5/zS61WQ6vVGtwcHBzq9b50s1bO9nhucEf8+uK9WBAfgrYtHUVHIpn4+c88vPx17e8rkTHJbrSeq6srXn/9dURFRSEwMBCHDx9G5861F9bT6XRwd3fH4sWLMXHixFu+fs2aNZg+fTry8/Mb9L4TJkyo1/MtZbReQ+hqJGw/mosPdp3GgXPXRMchGfjXfYGYOrCD6BgkgMWP1tPpdEhKSkJJSQkiIyNRUVF7sPXvWzwqlQp2dnbYvXv3bZdVXFwMPz8/eHt7Y/jw4UhLS6v3+9aloqIChYWFBjdrpVIqMCSkDTZO6YONUyIxpHMbDp6wcq9/fxyb086LjkEWRHg5ZWRkwNnZGfb29pg8eTI2b96M4OBgBAUFwc/PD7NmzcK1a9dQWVmJRYsWIScnB1qtts7lBQUFYc2aNfj666+xfv16ODg4oG/fvjcdq6rrfeuSmJgIjUajv/n4+BjtMzBn3f1c8d6Y7vjl/wZgTG8/ONgK/5EiQV746g/8evKy6BhkIYTv1qusrERWVhby8/OxceNGfPjhh9ixYweCg4Nx4MABTJgwAYcOHYJKpUJ0dDSUytovv++++65ey6+pqUG3bt3Qv39/LFu2rF7veysVFRX6rTmgdlPWx8fHKnfr3c7Vkkp8uvcc1v56Fld4uQWr42Jvgw1TIhHUhr8T1sJqZoiIjo6Gv78/Vq1apb+voKAAlZWVcHNzQ69evRAREYHly5fXe5lPPPEEzp8/f9sRebd639ux5mNO9VFepcNXB87jte+OobSSl1ywJh4aB2x6qg88NLxkizWw+GNO10mSZLCFAgAajQZubm7IzMxEamoq4uLi7mp56enp8PC4/Xxgt3pfajgHWxVG9/bDhsmRaKN2uPMLyGJoC8ox/uP9KCrnnI3UcEJnJZ89ezZiY2Ph4+ODoqIiJCUlISUlBcnJyQCADRs2wM3NDb6+vsjIyMC0adMQHx+PmJgY/TISEhLg5eWFxMREAMC8efPQu3dvBAQEoLCwEMuWLUN6errBltad3peMp7OnBlum9sXja/bjqNZ6B5FYmz9zijDl04P4aFwPnsRNDSK0nHJzczFmzBhotVpoNBqEhoYiOTkZgwcPBgBotVrMmDEDubm58PDwQEJCAubOnWuwjKysLP1xKADIz8/HpEmTkJOTA41Gg/DwcOzcuRM9e/as9/uScbXROGDD5Eg8uz4NP/2ZJzoONZHdJy9j/jdHMT8+RHQUMkOyO+ZkLnjM6e7paiTM/+Yo1vx6VnQUakL/+Uco/hnB0a2WymqOOZHlUikV+Pf9nTHv/s5Q8rwoqzFny2Ecys4XHYPMDMuJmtzYPm3x4dgIONmpREehJlBZXYPJnx7A5WIOOKL6YzmREIOCWmPD5D4cyWcltAXlmPrZQVRxkliqJ5YTCRPsqcaWqX3R2ZPH7KzBvjNX8dp3x0THIDPBciKh2mgc8OWTkYju5C46CjWBj/ec5Rx8VC8sJxLOyd4Gq8ZEYHzftqKjUBN4cWMGDl8oEB2DZI7lRLKgUirwyojOeDWOI/ksXUV1DZ5cdwBXOfci3QbLiWQlIbItVo/twZF8Fu5CfhmeWX+QV9GlOrGcSHYGBrlzJJ8V2HPyCv7z/XHRMUimWE4kSxzJZx3e33kaXx+6KDoGyRDLiWSLI/msw8yv/sAxTgpMN2A5kaxdH8n3eN92oqOQiZRV6fDkugPIL+UACfoLy4lkT6VU4OURwZjPkXwWK+tqKZ5NSoeuhvNQUy2WE5mNMZFtsXocR/JZqp0nLuHNHzhAgmqxnMisDAysHcnnoeFIPku0IuUUtmVoRccgGWA5kdm5PpIvxIsj+SzR/204hBO5RaJjkGAsJzJLrdXXR/K1Fh2FjKy0snaAREFZlegoJBDLicyWo50NVo3pjgn9OJLP0py5XILnvkhHDQdIWC2WE5k1lVKBucODMT8+hCP5LMzPf+bh7Z8yRccgQVhOZBHG9PbjSD4LtOynTPxwJEd0DBKA5UQWY2CgO76awpF8lmbGl4dwMq9YdAxqYiwnsiidPDiSz9IUV1TjyXWpKCqX9wAJzrBuXApJknjEsQEKCwuh0WhQUFAAtZpfhHJTWlmNaUnp2H40V3QUMpKY4NZ4b3R3KJv44KIkSSgsq4a2sAw5BeXILSyH9n//zSn46/+LK6qxbkIv9G7fsknziWaq70KWUwOxnORPVyMh8btj+HD3GdFRyEj+b3BHPHNvgNGWp6uRcKmoAjmF5cgpqC2fnMKK2v//X/nkFJajvKp+W0UtHG2xdWo/+LZ0NFpGuWM5yQzLyXys23sO//76COdtswAKBfDR2B4YGHTnmerLKnX6grlpa6ewHLkF5cgrKoexfyw6tnbGxil94OJga9wFyxTLSWZYTuYl5Xgenv48DcUV1aKjUCO5ONhg/RO9oVIq9Fs2OQV/beVcLyKRJ/HeG+SO9xMioLKC8xtYTjLDcjI/x7SFmLBmPy4WlIuOQlZg8j3+eDE2SHQMkzPVdyFH65HVuD6Sr4uXRnQUsgLv7TiFzWnnRccwWywnsiruagd88WRvxARzTj4yvZkbM3Aw65roGGaJ5URWx9HOBitHd8cTUZyTj0yrsroGkz45AG1BmegoZoflRFZJpVTgpWHBWBAfYhUHrUmcy8UVeOKTVJRV6kRHMSssJ7Jqo3v74aNxPeBsbyM6ClmwwxcK8fyGQ+D4s/pjOZHVu6ejG76aEglPzslHJvRthhbLfjopOobZYDkRAQhqUzuSL8DdWXQUsmBv/XiCl6GvJ5YT0f+4qx3wxkNhUPAQFJnQjC8P4cjFAtExZI/lRPQ3YT7NMaa3n+gYZMHKqnR4Ym0qLhVViI4iaywnohs8f18g3FzsRccgC3axoBxPrktFRTVH8NWF5UR0A7WDLV4eHiw6Blm4g1n5mL3pMEfw1YHlRHQLw0M9EBXQSnQMsnAbD57Hh7t4SZdbYTkR3YJCocCC+BDY2fBXhEzrtW3H8MufeaJjyA5/84jq4NfSCc8M7CA6Blk4SQKeWZ/GKY5uwHIiuo1J97SHv5uT6Bhk4YorqnHwXL7oGLLCciK6DXsbFRbEdxEdg6xAZl6R6AiywnIiuoNI/5YY2c1LdAyycJm5xaIjyArLiageXhraCZpmtqJjkAXjlpMhlhNRPbR0tscsK7jkNolz5nIJqnQ1omPIhtByWrlyJUJDQ6FWq6FWqxEZGYlt27bpH8/NzcW4cePg6ekJR0dHDBkyBJmZmbdd5po1a6BQKG66lZeX65+TmJiIHj16wMXFBe7u7oiPj8fx48dNtp5kGf4Z4YMIvxaiY5CFqtJJOHelRHQM2RBaTt7e3li0aBFSU1ORmpqKQYMGIS4uDkeOHIEkSYiPj8fp06exdetWpKWlwc/PD9HR0Sgpuf1foFqthlarNbg5OPx1OYQdO3Zg6tSp2Lt3L7Zv347q6mrExMTccblk3ZRKBRY8EAIbXpyQTOQEjzvpKSSZzZ3h6uqK119/HVFRUQgMDMThw4fRuXNnAIBOp4O7uzsWL16MiRMn3vL1a9aswfTp05Gfn1/v97x06RLc3d2xY8cO9O/fv16vKSwshEajQUFBAdRqdb3fi8xf4rZjWLXjtOgYZIGmRwdgenRH0THuiqm+C2VzzEmn0yEpKQklJSWIjIxERUXtjL1/3+JRqVSws7PD7t27b7us4uJi+Pn5wdvbG8OHD0daWtptn19QUDt9vaurayPXgqzBtHsD4NW8megYZIEy87jldJ3wcsrIyICzszPs7e0xefJkbN68GcHBwQgKCoKfnx9mzZqFa9euobKyEosWLUJOTg602rov1hUUFIQ1a9bg66+/xvr16+Hg4IC+ffvWeaxKkiTMmDED/fr1Q0hISJ3LraioQGFhocGNrJOjnQ1ejessOgZZoMxcjti7Tng5BQYGIj09HXv37sWUKVMwduxYHD16FLa2tti4cSNOnDgBV1dXODo6IiUlBbGxsVCpVHUur3fv3hg9ejTCwsIQFRWFL7/8Eh07dsQ777xzy+c//fTT+OOPP7B+/frb5kxMTIRGo9HffHx8GrXeZN7u7dQa93VuLToGWRiO2PuL7I45RUdHw9/fH6tWrdLfV1BQgMrKSri5uaFXr16IiIjA8uXL673MJ554AufPnzcYCQgAzzzzDLZs2YKdO3eiXbt2t11GRUWFflcjULuf1cfHh8ecrJi2oAzRb+5ASSWvyUPG8+OM/ujg7iI6Rr1Z/DGn6yRJMigBANBoNHBzc0NmZiZSU1MRFxd3V8tLT0+Hh4eHwX1PP/00Nm3ahJ9//vmOxQQA9vb2+iHv129k3Tw0zfDcYPM6eE3yx5kiatmIfPPZs2cjNjYWPj4+KCoqQlJSElJSUpCcnAwA2LBhA9zc3ODr64uMjAxMmzYN8fHxiImJ0S8jISEBXl5eSExMBADMmzcPvXv3RkBAAAoLC7Fs2TKkp6cbbGlNnToVn3/+ObZu3QoXFxfk5OQAqC3BZs14oJvqb1yftth08AKOankMkozjRG4xYjmdo9hyys3NxZgxY6DVaqHRaBAaGork5GQMHjwYAKDVajFjxgzk5ubCw8MDCQkJmDt3rsEysrKyoFT+tQGYn5+PSZMmIScnBxqNBuHh4di5cyd69uypf87KlSsBAAMGDDBY1scff4xx48aZZmXJItmolFj4QAhGrvwV8tpBTubqBKcxAiDDY07mguc50d/N2ZKBT/dmiY5BFiCwtQu+f65+51vKgdUccyIyR/+6LwitnO1FxyALcPpyMUfsgeVEZBSaZraYO7yT6BhkAWrn2CsVHUM4lhORkdwf5ol+HVqJjkEWgCfjspyIjEahUGB+fAjsbPhrRY3DaYxYTkRG1a6VE6YO6CA6Bpm5E9xyYjkRGdvkAe3RvpWT6BhkxngiLsuJyOjsbVRYEF/3JMJEd+Ku5shPlhORCfTp0AoPhHuJjkFmqJWzPd78Z5joGMKxnIhMZPbQTlA7CJ2EhcyMQgEse7gr3F0c7vxkC8dyIjIRNxd7vBjLc5+o/p4dFIA+PB0BAMuJyKQe7uGDcN/momOQGYhs3xLP3hsgOoZssJyITEipVOC1B7pApVSIjkIy1srZHksf6cqfk79hORGZWCcPNSb0u/M1w8g6KRTAUh5nugnLiagJTLs3AJ4afvnQzZ4ZFIC+PM50E5YTURNwsrfBvDie+0SGerd3xTQeZ7ollhNRExkc3BqDgtxFxyCZaOVsh2UPh/M4Ux1YTkRNaPbQIH4ZERQK4K1RXeGu5q7eujS4nE6dOoU5c+bgkUceQV5eHgAgOTkZR44cMVo4IkvTwd0FD/fwER2DBHtmYAdEBbiJjiFrDSqnHTt2oEuXLti3bx82bdqE4uLaSQr/+OMPvPLKK0YNSGRppkd3hJOdSnQMEqRXO1dMi+4oOobsNaicXnzxRSxYsADbt2+HnZ2d/v6BAwfit99+M1o4Ikvk5mKPKQP8RccgAVo62WHZIzzOVB8NKqeMjAw88MADN93v5uaGK1euNDoUkaWb0K89PDi03KpcP87UmseZ6qVB5dS8eXNotdqb7k9LS4OXF2diJrqTZnYqPB8TKDoGNaGO7i6ICuD5TPXVoHJ69NFHMXPmTOTk5EChUKCmpgZ79uzB888/j4SEBGNnJLJID4R7IdhDLToGNZHjuUX4MjVbdAyz0aByWrhwIXx9feHl5YXi4mIEBwejf//+6NOnD+bMmWPsjEQWSalUYM4wzlpuTRZ+ewx5heWiY5gFhSRJUkNffOrUKaSlpaGmpgbh4eEICLCeM50LCwuh0WhQUFAAtZr/+qWGe3zNfvz8Z57oGNREhnZpgxWPdRcdw2hM9V3YqCuh+fv7w9+fo46IGmNWbBBSjuehpsH/TCRz8l1GDr4/koP7OrcRHUXW6l1OM2bMqPdClyxZ0qAwRNYooLULHu7pi8/3ZYmOQk1k7pbD6N2+JTTNbEVHka16l1NaWprBnw8cOACdTofAwNoRRydOnIBKpUL37pazuUrUVKZHB2Br2gWUVOpER6EmkFdUgcXJf+K1B7qIjiJb9S6nX375Rf//S5YsgYuLC9auXYsWLVoAAK5du4bx48cjKirK+CmJLJy7iwMm3+OPN7efEB2Fmsjn+7IQF+aJXu1bio4iSw0aEOHl5YUffvgBnTt3Nrj/8OHDiImJwcWLF40WUK44IIKMraxShwFv/ILcwgrRUaiJtG/lhO+mRcHB1nynszLVd2GDhpIXFhYiNzf3pvvz8vJQVFTU6FBE1ogn5lqf05dL8O7PJ0XHkKUGldMDDzyA8ePH46uvvsL58+dx/vx5fPXVV5gwYQJGjhxp7IxEVmNkN2904om5VuW9HadwTFsoOobsNKic3nvvPQwbNgyjR4+Gn58f/Pz88NhjjyE2NhYrVqwwdkYiq6FSKvDSUJ6Ya02qayS8uPEP6HgugYFGnYRbUlKCU6dOQZIkdOjQAU5OTsbMJms85kSmNO7j35Fy/JLoGNSE5g4PxoR+7UTHuGuyOuZ0nZOTE0JDQxEWFmZVxURkarNiO4FXVbAub3x/HNlXS0XHkI0GzRAxcOBAKBR1/+b8/PPPDQ5EREBgGxeM6uGD9b9zolBrUValw+zNGfjk8Z63/X61Fg3acuratSvCwsL0t+DgYFRWVuLgwYPo0oUnlREZw3PRHeHIK+ZalV2Zl7E57YLoGLLQoC2nt95665b3//vf/9Zfsp2IGsdd7YAn+/vjrR95Yq41efWbo+jf0Q2tnO1FRxGqUcecbjR69Gh89NFHxlwkkVV7on87uLtY95eUtckvrcL8b46KjiGcUcvpt99+g4MDL0FMZCyOdjY8MdcKbU2/iF+s/DIqDdqtd+OJtpIkQavVIjU1FXPnzjVKMCKq9WB3b3y05wz+zOHsK9bkpc0Z+GHGPXC2b9SVjcxWg7ac1Go1NBqN/ubq6ooBAwbgu+++wyuvvGLsjERWTaVUYDZPzLU6FwvK8cb3x0XHEKZRJ+FaM56ES00t4aPfsfMET8y1JgoFsHFKH3TzbSE6Sp1kdRJu+/btceXKlZvuz8/PR/v27RsdiohuNntoEE/MtTKSBLy48Q9UVteIjtLkGlROZ8+ehU5380XRKioqcOECx+gTmUJQGzUe6u4jOgY1sRO5xViZckp0jCZ3V0favv76a/3/f//999BoNPo/63Q6/PTTT2jbtq3RwhGRof+L6YivD11EWRWvmGtN3v0lE0O7tEFAaxfRUZrMXR1zUiprN7QUCgVufJmtrS3atm2LN998E8OHDzduShniMScS5e0fT+DtHzNFx6Am1t2vBTY8GQmlzPbtyuKYU01NDWpqauDr64u8vDz9n2tqalBRUYHjx4/fVTGtXLkSoaGhUKvVUKvViIyMxLZt2/SP5+bmYty4cfD09ISjoyOGDBmCzMzb/1KuWbMGCoXiplt5ebn+OTt37sSIESPg6ekJhUKBLVu23M3HQCTUpP7teWKuFTpw7ho+23dOdIwm06BjTmfOnEGrVq0a/ebe3t5YtGgRUlNTkZqaikGDBiEuLg5HjhyBJEmIj4/H6dOnsXXrVqSlpcHPzw/R0dEoKSm57XLVajW0Wq3B7e8nB5eUlCAsLAzvvvtuo9eBqKk52tng/2I6io5BAixOPo6L+WWiYzSJeu/WW7ZsGSZNmgQHBwcsW7bsts999tlnGxzI1dUVr7/+OqKiohAYGIjDhw+jc+fOAGqPa7m7u2Px4sWYOHHiLV+/Zs0aTJ8+Hfn5+fV6P4VCgc2bNyM+Pv6ucnK3Homkq5EwdOkuHM/libnW5t4gd3w4NkI2M5eb6ruw3gMi3nrrLTz22GNwcHCoc+JXoPbLviHlpNPpsGHDBpSUlCAyMhIVFRUAYLDFo1KpYGdnh927d9dZTgBQXFwMPz8/6HQ6dO3aFfPnz0d4ePhdZ/q7iooKfSag9i+ESBSVUoHZwzph7Ee/i45CTeynP/NwTFuEYE/L/kdxvcvpzJkzt/z/xsrIyEBkZCTKy8vh7OyMzZs3Izg4GFVVVfDz88OsWbOwatUqODk5YcmSJcjJyYFWq61zeUFBQVizZg26dOmCwsJCLF26FH379sWhQ4cQEBDQ4JyJiYmYN29eg19PZGz3dHRDVEAr7Mq8LDoKNbETuZZfTg065vTqq6+itPTmKzaWlZXh1VdfvatlBQYGIj09HXv37sWUKVMwduxYHD16FLa2tti4cSNOnDgBV1dXODo6IiUlBbGxsVCp6r7GTe/evTF69GiEhYUhKioKX375JTp27Ih33nnnrtfz72bNmoWCggL9LTubF4Ej8WYP7QSZ7N2hJpSZZ/m7cxtUTvPmzbvldZtKS0vveuvCzs4OHTp0QEREBBITExEWFoalS5cCALp374709HTk5+dDq9UiOTkZV65cQbt27eq9fKVSiR49etxxlN+d2Nvb60cVXr8RidbJQ42HunuLjkFNLDPX8q+b16BykiTplgfjDh06BFdX10YFkiTJ4NgOAGg0Gri5uSEzMxOpqamIi4u7q+Wlp6fDw8OjUbmI5GrG4EA0s+UVc63JyTzLL6e7miGiRYsW+vOGOnbsaFBQOp0OxcXFmDx5cr2XN3v2bMTGxsLHxwdFRUVISkpCSkoKkpOTAQAbNmyAm5sbfH19kZGRgWnTpiE+Ph4xMTH6ZSQkJMDLywuJiYkAarfqevfujYCAABQWFmLZsmVIT0/H8uXL9a8pLi7GyZMn9X8+c+YM0tPT4erqCl9f37v5SIiEa6NxwBP922PZTzwx11qcvVKCimod7G0s9x8ld1VOb7/9NiRJwuOPP4558+YZTF9kZ2eHtm3bIjIyst7Ly83NxZgxY6DVaqHRaBAaGork5GQMHjwYAKDVajFjxgzk5ubCw8MDCQkJN10vKisrSz9zBVA7+eykSZOQk5MDjUaD8PBw7Ny5Ez179tQ/JzU1FQMHDtT/ecaMGQCAsWPHYs2aNXfzkRDJwpP92+PzfVm4XFxx5yeT2auRgNOXStDJw3IPLzTokhk7duxAnz59YGtra4pMZoHnOZHcrP89C7M2ZYiOQU1k2SPhuD/MU3QMeUxfdN0999yjL6aysjIUFhYa3Iio6T3U3RsdWzuLjkFN5KSFn4DdoHIqLS3F008/DXd3dzg7O6NFixYGNyJqejYqJWbxirlWI9PCB0U0qJz+9a9/4eeff8aKFStgb2+PDz/8EPPmzYOnpyc++eQTY2ckonoa0NEN/To0ft5Lkj+W0y3897//xYoVK/CPf/wDNjY2iIqKwpw5c/Daa6/hs88+M3ZGIqonhUKBWUODeGKuFTh7ucSir5DboHK6evWq/kRYtVqNq1evAgD69euHnTt3Gi8dEd21zp4aPNiNJ+ZauuoaCeeu3P4KDeasQeXUvn17nD17FgAQHByML7/8EkDtFtXfh5cTkRj/F9MRDrYN+vUmM3LCgmeKaNBP7/jx43Ho0CEAtXPOXT/29Nxzz+GFF14wakAiunsemmZ4Iqq96BhkYpY8x95dnYR73XPPPaf//4EDB+LPP/9Eamoq3Nzc8PHHHxstHBE13JP3+GP971m4XFwpOgqZiCUPijDKdr+vry9GjhwJtVqNtWvXGmORRNRIzvY2mDqwg+gYZEInuVuPiMzRPyN8oGlmvTO5WLrTl4tRrbPMEXssJyIL5mRvg9G9OZmxparSSTh39eZr61kClhORhRvbpy3sVPxVt1SZFjqN0V0NiBg5cuRtH8/Pz29MFiIyAXcXB4zs5oWk/bx6syXKzC3GkBDRKYzvrsrpTucwaTQaJCQkNCoQERnfxKj2LCcLZakj9u6qnDhMnMg8dXB3RnQnd/x4LE90FDIySy0n7ogmshKT+vuLjkAmcOpSMXQ1d31ZPtljORFZiR5tW6CrT3PRMcjIKqtrkG2BI/ZYTkRWQqFQYFJ/TmlkiSxx1x7LiciK3Ne5DfxaOoqOQUZ2wgKHk7OciKyISqnAxH7tRMcgIzvJLSciMnf/6O6DFo6c0siSWOLs5CwnIivTzE6FMZFtRccgIzqZV4waCxuxx3IiskJjI/1gb8Nff0tRXlWDC/llomMYFX86iaxQS2d7/KM7L+VuSSxt1x7LichKTYxqD4VCdAoyFku7ZDvLichKtWvlhJjg1qJjkJFkspyIyFJwSiPLcZK79YjIUnT3a4Hufi1ExyAjyMwrhiRZzog9lhORleOURpahtFKHiwXlomMYDcuJyMoN7tQa7Vo5iY5BRmBJV8VlORFZOaVSgYlRnNLIEljSoAiWExHhwW7eaOlkJzoGNdLZKyWiIxgNy4mI4GCrwtg+bUXHoEaK9G8pOoLRsJyICAAwprcfHGz5lWCuWjjaYrAFnbfGn0QiAgC0cLLDPyN8RMegBnqwmzfsbVSiYxgNy4mI9Cb2aw8lpzQySw/3tKx/WLCciEjPt6UjYkM8RMeguxTh1wId3F1ExzAqlhMRGXiCJ+WanVE9LGurCWA5EdENuvo0R892rqJjUD252NtgWKjlbe2ynIjoJk9y68ls3N/VE452NqJjGB3LiYhuMjDQHf5unNLIHDzS01d0BJNgORHRTZRKBSeENQOdPdUI8dKIjmESLCciuqX4cC+4udiLjkG38bCFbjUBLCciqoO9jQpzhnUSHYPq4GCrRFxXT9ExTIblRER1iuvqhacG8Gq5cjSsiyfUDraiY5gMy4mIbuv5mECLmrPNUjxiYTNC3EhoOa1cuRKhoaFQq9VQq9WIjIzEtm3b9I/n5uZi3Lhx8PT0hKOjI4YMGYLMzMzbLnPNmjVQKBQ33crLDa8QuWLFCrRr1w4ODg7o3r07du3aZZJ1JDJ3SqUCb4/qiqA2ljUDgTnr4O6M7n4tRMcwKaHl5O3tjUWLFiE1NRWpqakYNGgQ4uLicOTIEUiShPj4eJw+fRpbt25FWloa/Pz8EB0djZKS21+zRK1WQ6vVGtwcHBz0j3/xxReYPn06XnrpJaSlpSEqKgqxsbHIysoy9SoTmSUnext8ODaC13ySiYd7+EChsOxJEBWSJEmiQ/ydq6srXn/9dURFRSEwMBCHDx9G586dAQA6nQ7u7u5YvHgxJk6ceMvXr1mzBtOnT0d+fn6d79GrVy9069YNK1eu1N/XqVMnxMfHIzExsV45CwsLodFoUFBQALVaXf8VJDJjqWev4pEP9qJKJ6uvDatiq1Jg76x70dJZHiMpTfVdKJtjTjqdDklJSSgpKUFkZCQqKioAwGCLR6VSwc7ODrt3777tsoqLi+Hn5wdvb28MHz4caWlp+scqKytx4MABxMTEGLwmJiYGv/76qxHXiMjyRLR1xWsPdBEdw6rFdG4jm2IyJeHllJGRAWdnZ9jb22Py5MnYvHkzgoODERQUBD8/P8yaNQvXrl1DZWUlFi1ahJycHGi12jqXFxQUhDVr1uDrr7/G+vXr4eDggL59++qPVV2+fBk6nQ6tWxse4G3dujVycnLqXG5FRQUKCwsNbkTW6KEIH56gK9DDFjjJ660IL6fAwECkp6dj7969mDJlCsaOHYujR4/C1tYWGzduxIkTJ+Dq6gpHR0ekpKQgNjYWKlXdF9Tq3bs3Ro8ejbCwMERFReHLL79Ex44d8c477xg878b9tZIk3XYfbmJiIjQajf7m42MdPyBEtzJzSBAGBbmLjmF1vFs0Q1//VqJjNAnh5WRnZ4cOHTogIiICiYmJCAsLw9KlSwEA3bt3R3p6OvLz86HVapGcnIwrV66gXbt29V6+UqlEjx499FtOrVq1gkqlumkrKS8v76atqb+bNWsWCgoK9Lfs7OwGrC2RZVApFVj6cFcEuDuLjmJVRkX4QGklV4MUXk43kiRJf7zpOo1GAzc3N2RmZiI1NRVxcXF3tbz09HR4eNROKW9nZ4fu3btj+/btBs/bvn07+vTpU+dy7O3t9UPer9+IrJmLgy1Wj+2BFo6WeyKonCgVtbtUrYXQedZnz56N2NhY+Pj4oKioCElJSUhJSUFycjIAYMOGDXBzc4Ovry8yMjIwbdo0xMfHGwxmSEhIgJeXl36U3bx589C7d28EBASgsLAQy5YtQ3p6OpYvX65/zYwZMzBmzBhEREQgMjIS77//PrKysjB58uSm/QCIzJxvS0esHN0doz/ch+oajuAzpYGB7mijcbjzEy2E0HLKzc3FmDFjoNVqodFoEBoaiuTkZAwePBgAoNVqMWPGDOTm5sLDwwMJCQmYO3euwTKysrKgVP61AZifn49JkyYhJycHGo0G4eHh2LlzJ3r27Kl/zqhRo3DlyhW8+uqr0Gq1CAkJwXfffQc/P7+mWXEiC9K7fUvMjw/BrE0ZoqNYNEue5PVWZHeek7ngeU5Ehub99wg+3nNWdAyLteuFgfBxdRQd4yYWf54TEZm3l4Z2Qv+ObqJjWKzMvCLREZoUy4mIjMJGpcQ7j4SjPa+gaxLHc4pFR2hSLCciMhpNs9oRfJpmHMFnbCdyueVERNRg7Vo5YcVj3aCykvNxmsrxHJYTEVGj9O3QCv8eESw6hkU5eakY1boa0TGaDMuJiExiTGRbjOnN0zOMpbK6BueuloqO0WRYTkRkMi+PCEbfDi1Fx7AYJ6xo1x7LiYhMxlalxPJHu6FtS/mdn2OOjlvRoAiWExGZVHNHO3w4tgdcHIROSGMRrGnEHsuJiEyug7sz3n20GziAr3GsacQey4mImsQ9Hd0wZxhH8DXG2SulKK/SiY7RJFhORNRkxvdtazVXcjUFXY2E05dKRMdoEiwnImoyCoUCr8aFoGc7V9FRzJa1HHdiORFRk7KzUeK90d3h49pMdBSzZC0j9lhORNTkXJ3ssHpsDzjZqURHMTvWcq4Ty4mIhOjY2gXLHgmHgiP47soJK7l0BsuJiIS5t1NrzIoNEh3DrGRfLUNJRbXoGCbHciIioZ6Iao8Hu3mLjmFWMvMs/9pOLCciEkqhUOC1kSHo7tdCdBSzYQ3HnVhORCScvY0Kq8Z0h1dzjuCrD2sYscdyIiJZaOVsjw8SIuDIEXx3ZA3nOrGciEg2gj3VeGtUV9ExZM8a5thjORGRrNzXuQ3+dV+g6BiylldUgWsllaJjmBTLiYhk56kB/ojv6ik6hqxZ+q49lhMRyY5CocCiB0MR5tNcdBTZ+tPCd+2xnIhIlhxsVfhgTHd4aBxER5Glb/64KDqCSbGciEi23NUO+CAhAg62/Kq60f6z13D0YqHoGCbDv3EikrUQLw2W/LOr6BiytG7vWdERTIblRESyN7SLB56L7ig6huxsSbuIgtIq0TFMguVERGbhmUEd0M23uegYslJWpcOGA9miY5gEy4mIzIJSqUDiyFDYKHmNjb/7dO851NRIomMYHcuJiMxGYBsXTL7HX3QMWTl7pRQ7My+JjmF0LCciMitPD+qAdq2cRMeQlXW/nRMdwehYTkRkVhxsVVgYHyI6hqz8fDwP2VdLRccwKpYTEZmdPh1a4R/deYHC6ySp9tiTJWE5EZFZemloJ7g62YmOIRtfpGajvEonOobRsJyIyCy1cLLDy8ODRceQjfzSKnx9yHKmNGI5EZHZiuvqiaiAVqJjyMYnv52FJFnGsHKWExGZLYVCgYXxXTj33v8cvlCItOx80TGMgn+jRGTWfFs6YjqnNtKzlGHlLCciMnsT+rVDJw+16Biy8O0fWlwurhAdo9FYTkRk9mxVSiwa2QUKzmyESl0Nvthv/vPtsZyIyCKE+TTHuD5tRceQhU/3nkO1rkZ0jEZhORGRxfi/mEBeOReAtqAcPx7LEx2jUVhORGQxnO1tMD+OUxsBtcPKzRnLiYgsSnRwawzt0kZ0DOF+PXUFJ/OKRMdoMJYTEVmcf4/oDBcHG9ExhPvEjIeVCy2nlStXIjQ0FGq1Gmq1GpGRkdi2bZv+8dzcXIwbNw6enp5wdHTEkCFDkJmZWe/lJyUlQaFQID4+3uD+oqIiTJ8+HX5+fmjWrBn69OmD/fv3G2u1iEgwd7UDXowNEh1DuI0HzqOo3Dwv4y60nLy9vbFo0SKkpqYiNTUVgwYNQlxcHI4cOQJJkhAfH4/Tp09j69atSEtLg5+fH6Kjo1FSUnLHZZ87dw7PP/88oqKibnps4sSJ2L59O9atW4eMjAzExMQgOjoaFy5cMMVqEpEAj/TwRYRfC9ExhCqp1GFzmnl+rykkmU3E5Orqitdffx1RUVEIDAzE4cOH0blzZwCATqeDu7s7Fi9ejIkTJ9a5DJ1Oh3vuuQfjx4/Hrl27kJ+fjy1btgAAysrK4OLigq1bt2LYsGH613Tt2hXDhw/HggUL6pWzsLAQGo0GBQUFUKt58h+RHJ3ILcKwZbtQpZPV11yT6uDujO3P9YfCRCeBmeq7UDbHnHQ6HZKSklBSUoLIyEhUVNSe4ezg8NewUJVKBTs7O+zevfu2y3r11Vfh5uaGCRMm3PRYdXU1dDqdwXIBoFmzZrddbkVFBQoLCw1uRCRvHVu7YIqVX9b9ZF4xfjt9RXSMuya8nDIyMuDs7Ax7e3tMnjwZmzdvRnBwMIKCguDn54dZs2bh2rVrqKysxKJFi5CTkwOtVlvn8vbs2YPVq1fjgw8+uOXjLi4uiIyMxPz583Hx4kXodDp8+umn2Ldv322Xm5iYCI1Go7/5+Pg0et2JyPSeGtgB7a38su6f/Gp+AyOEl1NgYCDS09Oxd+9eTJkyBWPHjsXRo0dha2uLjRs34sSJE3B1dYWjoyNSUlIQGxsLlUp1y2UVFRVh9OjR+OCDD9CqVd3T6K9btw6SJMHLywv29vZYtmwZHn300TqXCwCzZs1CQUGB/padbf7TgxBZAwdbFV4b2UV0DKG2H8vFxfwy0THuiuyOOUVHR8Pf3x+rVq3S31dQUIDKykq4ubmhV69eiIiIwPLly296bXp6OsLDww1KpqamdgoPpVKJ48ePw9//r038kpISFBYWwsPDA6NGjUJxcTG+/fbbeuXkMSci8/LCV4fwZep50TGEeXpgBzx/X6DRl2vxx5yukyRJf7zpOo1GAzc3N2RmZiI1NRVxcXG3fG1QUBAyMjKQnp6uv91///0YOHAg0tPTb9oV5+TkBA8PD1y7dg3ff/99ncslIvM3e2gntLTiy7on7c9CRbX5XMZd6Flqs2fPRmxsLHx8fFBUVISkpCSkpKQgOTkZALBhwwa4ubnB19cXGRkZmDZtGuLj4xETE6NfRkJCAry8vJCYmAgHBweEhBhOXdK8eXMAMLj/+++/hyRJCAwMxMmTJ/Gvf/0LgYGBGD9+vOlXmoiEaO5oh5dHBGNaUrroKEJcLq7EtowcxId7iY5SL0LLKTc3F2PGjIFWq4VGo0FoaCiSk5MxePBgAIBWq8WMGTOQm5sLDw8PJCQkYO7cuQbLyMrKglJ5dxuABQUFmDVrFs6fPw9XV1c8+OCDWLhwIWxtbY22bkQkP/eHeWLTwQvYceKS6ChCfPLbWbMpJ9kdczIXPOZEZJ6yr5Yi5q2dKKsyn11cxvTNM/0Q4qUx2vKs5pgTEZEp+bg64rnBAaJjCGMus5WznIjI6jzetx2CrfSy7lvTL+JaSaXoGHfEciIiq2OjUmLRg12gtMLLuldU12DDAfmfp8lyIiKrFOrdHOP7thMdQ4h1e89BVyPv4QYsJyKyWjMGd4RX82aiYzS57Ktl2HFC3pdxZzkRkdVysrfB/PjOomMIIfcLEbKciMiqDQpqjWGhHqJjNLmU45dw9vKdr40nCsuJiKzeKyOCrfKy7p/ule/WE8uJiKyeu4sDZg/tJDpGk/syNRtllfI8GZnlREQEYFSED3q2dRUdo0kVlldja7o8L+POciIiAqBUKvDayBDYqqzr5Ke1v52DHGexYzkREf1PB3cXPDWgg+gYTeqYthAHzl0THeMmLCcior95aqA/2rtZ12Xd18pwWDnLiYjob+xtVEh8wLou656WdU12u/ZYTkREN+jVviUe6elz5ydaABulAksf7gqFQl7H2lhORES38OKQTmjlbC86hsnNHtoJ3f3kN0qR5UREdAsaR1u8MiJYdAyTGhbqgfF924qOcUssJyKiOgwP9cDAQDfRMUyivZsTFj8YKrvdedexnIiI6qBQKDA/PgTNbFWioxhVM1sV3hvdHc728p2yieVERHQb3i0c8X8xHUXHMKpFD3ZBx9YuomPcFsuJiOgOxvVpixAvy7ise0KkH+K6eomOcUcsJyKiO7BRKbFoZKjZX9a9q09zvDTMPCa4ZTkREdVDiJcGE/qZ72XdWzjaYvlj3WBvYx7Hz1hORET19JyZXtZdoQCWPhxuVtlZTkRE9eRoZ4MFD4SIjnHXpt/bEf07mteQeJYTEdFdGBjojhFhnqJj1NuAQDc8M8j8ZlpnORER3aWXhwdDbQaXdfdq3gxv/bMrlGY4koPlRER0l9xc7GU/6s1OpcSKx7qhhZOd6CgNwnIiImqAf0b4oFc7+U2Yet3LI4IR5tNcdIwGYzkRETWAQqHAayO7wE4lv6/RkeFeeKyXr+gYjSL/naZERDLl7+aMqQM74K0fTwjN4WJvgxAvDUJ9NAjzbo6Bge6yndC1vlhORESNMHlAe/z3j4s4mVfcJO9nb6NEZ081Qr2bI8xHg1Dv5mjX0sksBz3cDsuJiKgR7G1USBzZBQ+995vRl22jVCCwjQtCvWtLKNRbg46tXWArw12JxsZyIiJqpB5tXfFoL198vi+rwctQKID2rZwQ9r8SCvVpjmAPNRws7HId9cVyIiIygplDgrD9aC4uFVXU6/neLZr9VUTezRHipYaLg62JU5oPlhMRkRFomtni3yM6Y+rnB296zM3FHmH/K6Eu3hqEemnQ0tleQErzwXIiIjKSoV3a4P4wT1wtqdRvEYX5aNBG7WD2o+eaGsuJiMhIFAoFlj0SLjqGRbD8IR9ERGR2WE5ERCQ7LCciIpIdlhMREckOy4mIiGSH5URERLLDciIiItlhORERkewILaeVK1ciNDQUarUaarUakZGR2LZtm/7x3NxcjBs3Dp6ennB0dMSQIUOQmZlZ7+UnJSVBoVAgPj7e4P7q6mrMmTMH7dq1Q7NmzdC+fXu8+uqrqKmpMdaqERFRIwidIcLb2xuLFi1Chw4dAABr165FXFwc0tLSEBwcjPj4eNja2mLr1q1Qq9VYsmQJoqOjcfToUTg5Od122efOncPzzz+PqKiomx5bvHgx3nvvPaxduxadO3dGamoqxo8fD41Gg2nTpplkXYmIqP4UkiRJokP8naurK15//XVERUUhMDAQhw8fRufOnQEAOp0O7u7uWLx4MSZOnFjnMnQ6He655x6MHz8eu3btQn5+PrZs2aJ/fPjw4WjdujVWr16tv+/BBx+Eo6Mj1q1bV6+chYWF0Gg0KCgogFqtbtjKEhGZOVN9F8rmmJNOp0NSUhJKSkoQGRmJioraaecdHBz0z1GpVLCzs8Pu3btvu6xXX30Vbm5umDBhwi0f79evH3766SecOFF7aeVDhw5h9+7dGDp0qJHWhoiIGkP4xK8ZGRmIjIxEeXk5nJ2dsXnzZgQHB6Oqqgp+fn6YNWsWVq1aBScnJyxZsgQ5OTnQarV1Lm/Pnj1YvXo10tPT63zOzJkzUVBQgKCgIKhUKuh0OixcuBCPPPJIna+pqKjQFyZQ+68FIiIyDeFbToGBgUhPT8fevXsxZcoUjB07FkePHoWtrS02btyIEydOwNXVFY6OjkhJSUFsbCxUqltfGbKoqAijR4/GBx98gFatWtX5nl988QU+/fRTfP755zh48CDWrl2LN954A2vXrq3zNYmJidBoNPqbj49Po9ediIhuTXbHnKKjo+Hv749Vq1bp7ysoKEBlZSXc3NzQq1cvREREYPny5Te9Nj09HeHh4QbldX0EnlKpxPHjx+Hv7w8fHx+8+OKLmDp1qv55CxYswKeffoo///zzlrluteXk4+PDY05EZNVMdcxJ+G69G0mSZFACAKDRaAAAmZmZSE1Nxfz582/52qCgIGRkZBjcN2fOHBQVFWHp0qX6rZ3S0lIolYYbjSqV6rZDye3t7WFv/9eVK693OnfvEZE1u/4daOztHKHlNHv2bMTGxsLHxwdFRUVISkpCSkoKkpOTAQAbNmyAm5sbfH19kZGRgWnTpiE+Ph4xMTH6ZSQkJMDLywuJiYlwcHBASEiIwXs0b94cAAzuHzFiBBYuXAhfX1907twZaWlpWLJkCR5//PF6Zy8qKgIA7t4jIkLtd+L1DQljEFpOubm5GDNmDLRaLTQaDUJDQ5GcnIzBgwcDALRaLWbMmIHc3Fx4eHggISEBc+fONVhGVlbWTVtBd/LOO+9g7ty5eOqpp5CXlwdPT088+eSTePnll+u9DE9PT2RnZ8PFxcVsL798fddkdnY2d02aCD9j0+Lna3p3+owlSUJRURE8PT2N+r6yO+ZETYfnapkeP2PT4udreqI+Y+Gj9YiIiG7EciIiItlhOVkxe3t7vPLKKwajEMm4+BmbFj9f0xP1GfOYExERyQ63nIiISHZYTkREJDssJyIikh2WExERyQ7LyQIkJiZCoVBg+vTpBvcfO3YM999/PzQaDVxcXNC7d29kZWXVuZxNmzYhIiICzZs3h5OTE7p27XrLiy+uWLEC7dq1g4ODA7p3745du3YZe5VkpSk/38TERPTo0QMuLi5wd3dHfHw8jh8/borVkpWm/hm+0/taoqb+jC9cuIDRo0ejZcuWcHR0RNeuXXHgwIF655XdxK90d/bv34/3338foaGhBvefOnUK/fr1w4QJEzBv3jxoNBocO3bM4OKNN3J1dcVLL72EoKAg2NnZ4ZtvvsH48ePh7u6O++67D0Dt5UamT5+OFStWoG/fvli1ahViY2Nx9OhR+Pr6mnRdRWjqz3fHjh2YOnUqevTogerqarz00kuIiYnB0aNH4eTkZNJ1FaWpP+M7va8laurP+Nq1a+jbty8GDhyIbdu2wd3dHadOndLPdVovEpmtoqIiKSAgQNq+fbt0zz33SNOmTdM/NmrUKGn06NGNfo/w8HBpzpw5+j/37NlTmjx5ssFzgoKCpBdffLHR7yU3Ij7fG+Xl5UkApB07djT6veRI1Gd8u/e1NCI+45kzZ0r9+vVr1DK5W8+MTZ06FcOGDUN0dLTB/TU1Nfj222/RsWNH3HfffXB3d0evXr2wZcuWei9bkiT89NNPOH78OPr37w8AqKysxIEDBwxmhQeAmJgY/Prrr41eH7lp6s/3VgoKCgDU/mvVEon6jOt6X0sk4jP++uuvERERgYceegju7u4IDw/HBx98cHfBG1VtJMz69eulkJAQqaysTJIkyeBfRFqtVgIgOTo6SkuWLJHS0tKkxMRESaFQSCkpKbddbn5+vuTk5CTZ2NhI9vb20urVq/WPXbhwQQIg7dmzx+A1CxculDp27GjcFRRMxOd7o5qaGmnEiBGN/heoXIn6jG/3vpZG1Gdsb28v2dvbS7NmzZIOHjwovffee5KDg4O0du3aemdnOZmhrKwsyd3dXUpPT9ff9/cfuusl8sgjjxi8bsSIEdLDDz9822XrdDopMzNTSktLk9544w1Jo9FIv/zyi8Fyf/31V4PXLFiwQAoMDGz8ismEqM/3Rk899ZTk5+cnZWdnN2p95EjUZ3yn97UkIn+ObW1tpcjISIPXPPPMM1Lv3r3rnZ/lZIY2b94sAZBUKpX+BkBSKBSSSqWSysvLJRsbG2n+/PkGr3vhhRekPn363NV7TZgwQYqJiZEkSZIqKioklUolbdq0yeA5zz77rNS/f//GrZSMiPp8/+7pp5+WvL29pdOnTzdqXeRK1Gd8p/etrq422jqKJvLn2NfXV5owYYLBc1asWCF5enrWe5kcrWeG7r333psuRz9+/HgEBQVh5syZsLe3R48ePW4agnzixAn4+fnd1XtJkoSKigoAgJ2dHbp3747t27fjgQce0D9n+/btiIuLa+DayI+oz/f6n5955hls3rwZKSkpaNeuXcNXRMZEfcZ3el+VStWAtZEnkT/Hffv2bfxy76oeSbZu3DWxadMmydbWVnr//felzMxM6Z133pFUKpW0a9cu/XPGjBljMMrutddek3744Qfp1KlT0rFjx6Q333xTsrGxkT744AP9c5KSkiRbW1tp9erV0tGjR6Xp06dLTk5O0tmzZ5tkPUVpqs93ypQpkkajkVJSUiStVqu/lZaWNsl6itRUn/Gd3teSNdVn/Pvvv0s2NjbSwoULpczMTOmzzz6THB0dpU8//bTeWVlOFuJWv2CrV6+WOnToIDk4OEhhYWHSli1bbnrN2LFj9X9+6aWX9M9v0aKFFBkZKSUlJd30XsuXL5f8/PwkOzs7qVu3bhY7zPnvmurzBXDL28cff2yiNZOPpvwZvtP7Wqqm/Iz/+9//SiEhIZK9vb0UFBQkvf/++3eVlZfMICIi2eF5TkREJDssJyIikh2WExERyQ7LiYiIZIflREREssNyIiIi2WE5ERGR7LCciMzU2bNnoVAokJ6ebpLlKxSKu7p8ApExsZyIGmjcuHGIj48X9v4+Pj7QarUICQkBAKSkpEChUCA/P19YJiJj4cSvRGZKpVKhTZs2omMQmQS3nIhMYMeOHejZsyfs7e3h4eGBF198EdXV1frHBwwYgGeffRYvvPACXF1d0aZNG/z73/82WMaff/6Jfv36wcHBAcHBwfjxxx8NdrX9fbfe2bNnMXDgQABAixYtoFAoMG7cOABA27Zt8fbbbxssu2vXrgbvl5mZif79++vfa/v27Tet04ULFzBq1Ci0aNECLVu2RFxcHM6ePdvYj4rollhOREZ24cIFDB06FD169MChQ4ewcuVKrF69GgsWLDB43tq1a+Hk5IR9+/bhP//5D1599VV9KdTU1CA+Ph6Ojo7Yt28f3n//fbz00kt1vqePjw82btwIADh+/Di0Wi2WLl1ar7w1NTUYOXIkVCoV9u7di/feew8zZ840eE5paSkGDhwIZ2dn7Ny5E7t374azszOGDBmCysrKu/l4iOqFu/WIjGzFihXw8fHBu+++C4VCgaCgIFy8eBEzZ87Eyy+/DKWy9t+EoaGheOWVVwAAAQEBePfdd/HTTz9h8ODB+OGHH3Dq1CmkpKTod90tXLgQgwcPvuV7qlQquLq6AgDc3d3RvHnzeuf98ccfcezYMZw9exbe3t4AgNdeew2xsbH65yQlJUGpVOLDDz+EQqEAAHz88cdo3rw5UlJSEBMTc3cfEtEdsJyIjOzYsWOIjIzUf4kDtRdfKy4uxvnz5+Hr6wugtpz+zsPDA3l5eQBqt358fHwMjin17NnTZHl9fX31xQQAkZGRBs85cOAATp48CRcXF4P7y8vLcerUKZPkIuvGciIyMkmSDIrp+n0ADO63tbU1eI5CoUBNTU2dy2gopVKJG6+MU1VVdVO2G7P8XU1NDbp3747PPvvspue6ubkZJSfR37GciIwsODgYGzduNCiYX3/9FS4uLvDy8qrXMoKCgpCVlYXc3Fy0bt0aALB///7bvsbOzg4AoNPpDO53c3ODVqvV/7mwsBBnzpwxyJuVlYWLFy/C09MTAPDbb78ZLKNbt2744osv4O7uDrVaXa91IGoMDoggaoSCggKkp6cb3CZNmoTs7Gw888wz+PPPP7F161a88sormDFjhv54050MHjwY/v7+GDt2LP744w/s2bNHPyCiri0qPz8/KBQKfPPNN7h06RKKi4sBAIMGDcK6deuwa9cuHD58GGPHjoVKpdK/Ljo6GoGBgUhISMChQ4ewa9eumwZfPPbYY2jVqhXi4uKwa9cunDlzBjt27MC0adNw/vz5hnx0RLfFciJqhJSUFISHhxvcXnnlFXz33Xf4/fffERYWhsmTJ2PChAmYM2dOvZerUqmwZcsWFBcXo0ePHpg4caL+9Q4ODrd8jZeXF+bNm4cXX3wRrVu3xtNPPw0AmDVrFvr374/hw4dj6NChiI+Ph7+/v/51SqUSmzdvRkVFBXr27ImJEydi4cKFBst2dHTEzp074evri5EjR6JTp054/PHHUVZWxi0pMglepp3ITOzZswf9+vXDyZMnDcqFyBKxnIhkavPmzXB2dkZAQABOnjyJadOmoUWLFti9e7foaEQmxwERRDJVVFSEF154AdnZ2WjVqhWio6Px5ptvio5F1CS45URERLLDARFERCQ7LCciIpIdlhMREckOy4mIiGSH5URERLLDciIiItlhORERkeywnIiISHZYTkREJDv/D2zNwaEYk+NZAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# SET UP ##############################################\n", "\n", "# load city inputs files, to be updated for each city scan\n", "with open(\"city_inputs.yml\", 'r') as f:\n", " city_inputs = yaml.safe_load(f)\n", "\n", "city = city_inputs['city_name'].replace(' ', '_').lower()\n", "country = city_inputs['country_name'].replace(' ', '_').lower()\n", "# load global inputs, such as data sources that generally remain the same across scans\n", "with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", "# Read AOI shapefile --------\n", "# transform the input shp to correct prj (epsg 4326)\n", "aoi_file = gpd.read_file(city_inputs['AOI_path']).to_crs(epsg = 4326)\n", "features = aoi_file.geometry\n", "\n", "# Define output folder ---------\n", "output_folder = Path('mnt/city-directories/02-process-output')\n", "# Define render folder ---------\n", "render_folder = Path('mnt/city-directories/03-render-output')\n", "multi_scan_folder = Path('multi-scan-materials')\n", "\n", "if not os.path.exists(output_folder):\n", " os.mkdir(output_folder)\n", "\n", "fig, ax = plt.subplots()\n", "features.plot(ax=ax)\n", "plt.title('AOI')\n", "plt.xlabel('Longitude')\n", "plt.ylabel('Latitude')\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": 69, "id": "18954664-49f8-43b5-96c1-7faef1809401", "metadata": {}, "outputs": [], "source": [ "#AOI AREA\n", "def calculate_aoi_area(features):\n", " # Create a unit registry\n", " ureg = pint.UnitRegistry()\n", "\n", " # Combine geometries using unary_union from shapely\n", " combined_geometry = unary_union(features)\n", "\n", " # Calculate the area in square kilometers\n", " area_km2 = combined_geometry.area / 1e6 # Convert square meters to square kilometers\n", "\n", " # Print the result using the pint library for unit formatting\n", " area_quantity = area_km2 * ureg.km**2\n", " return area_quantity.to('km^2')\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 70, "id": "c1cc951b-ac81-470a-b400-dd38767a235b", "metadata": {}, "outputs": [], "source": [ "#Climate\n", "def get_koeppen_classification():\n", " # Load global inputs from YAML file\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " # Calculate centroid of AOI\n", " centroid = features.centroid.values[0]\n", "\n", " # Retrieve centroid coordinates\n", " coords = {' Lon': centroid.x, 'Lat': centroid.y}\n", "\n", " # Read Köppen climate classification file\n", " koeppen_file_path = global_inputs.get('koeppen_source') #use an alternate dataset\n", " koeppen = pd.read_csv(koeppen_file_path)\n", "\n", " # Filter Köppen data for the region around the centroid with a buffer of 0.5 degrees\n", " lon_min, lon_max = coords[' Lon'] - 0.5, coords[' Lon'] + 0.5\n", " lat_min, lat_max = coords['Lat'] - 0.5, coords['Lat'] + 0.5\n", " koeppen_city = koeppen[\n", " (koeppen[' Lon'].between(lon_min, lon_max)) &\n", " (koeppen['Lat'].between(lat_min, lat_max))\n", " ][' Cls'].unique()\n", "\n", " # Print Köppen climate classification\n", " koeppen_text = ', '.join(koeppen_city)\n", " print(f\"Köppen climate classification: {koeppen_text} (See https://en.wikipedia.org/wiki/Köppen_climate_classification for classes)\")\n", "\n", " # Return Köppen climate classification if needed\n", " return koeppen_city\n" ] }, { "cell_type": "code", "execution_count": 71, "id": "39d2f590-64e2-46d8-8fda-62b3eb6019b3", "metadata": {}, "outputs": [], "source": [ "#Age Distribution by World Pop\n", "def age_stats():\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " # Get city inputs\n", " city_inputs = global_inputs.get('city_inputs')\n", "\n", " if menu['demographics']: \n", " age_file = os.path.join(output_folder, f\"{city}_demographics.csv\") \n", " if os.path.exists(age_file): \n", " pass\n", " else:\n", " print(\"Demographics file does not exist.\")\n", " return \n", " pop_dist_group_wp = pd.read_csv(age_file)\n", "\n", " \n", " pop_dist_group_wp = pop_dist_group_wp.rename(columns={\"age_group\": \"Age_Bracket\", \"sex\": \"Sex\"})\n", " pop_dist_group_wp['Count'] = pd.to_numeric(pop_dist_group_wp['population'], errors='coerce')\n", " pop_dist_group_wp['Age_Bracket'] = pop_dist_group_wp['Age_Bracket'].replace({'<1': '0-4', '1-4': '0-4'})\n", " pop_dist_group_wp = pop_dist_group_wp.groupby(['Age_Bracket', 'Sex']).agg(Count=('Count', 'sum')).reset_index()\n", "\n", " \n", " pop_dist_group_wp['Percentage'] = pop_dist_group_wp.groupby('Sex')['Count'].transform(lambda x: x / x.sum())\n", " pop_dist_group_wp['Sexed_Percent'] = pop_dist_group_wp.groupby('Sex')['Count'].transform(lambda x: x / x.sum())\n", " pop_dist_group_wp['Sexed_Percent_cum'] = pop_dist_group_wp.groupby('Sex')['Sexed_Percent'].cumsum()\n", "\n", " # Plot age-sex distribution\n", " sns.barplot(data=pop_dist_group_wp, x='Age_Bracket', y='Percentage', hue='Sex', dodge=True)\n", " plt.title(f\"Population distribution in {city} by sex\")\n", " plt.xlabel(\"Age Bracket\")\n", " plt.ylabel(\"Percentage\")\n", " plt.legend(title=\"Sex\", loc=\"upper right\")\n", " plt.xticks(rotation=45)\n", " plt.tight_layout()\n", " render_path = os.path.join(render_folder, \"age_stats.png\")\n", " plt.savefig(render_path)\n", " plt.close()\n", "\n", " plt.show()\n", "\n", " # Plotly\n", " fig = px.bar(pop_dist_group_wp, x='Age_Bracket', y='Percentage', color='Sex', barmode='group', \n", " title=f\"Population distribution in {city} by sex\", \n", " labels={'Age_Bracket': 'Age Bracket', 'Percentage': 'Percentage', 'Sex': 'Sex'})\n", "\n", " \n", " fig.update_layout(xaxis_title=\"Age Bracket\", yaxis_title=\"Percentage\", legend_title=\"Sex\",plot_bgcolor='white') #swap the colors\n", " fig.update_xaxes(tickangle=45)\n", "\n", " \n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", " under5 = pop_dist_group_wp[pop_dist_group_wp['Age_Bracket'] == '0-4']['Percentage'].sum()\n", " youth = pop_dist_group_wp[pop_dist_group_wp['Age_Bracket'].isin(['15-19', '20-24'])]['Percentage'].sum()\n", " working_age = pop_dist_group_wp[pop_dist_group_wp['Age_Bracket'].isin(['15-64'])]['Percentage'].sum()\n", " elderly = pop_dist_group_wp[pop_dist_group_wp['Age_Bracket'].isin(['60-64', '65-69', '70-74', '75-79', '80+'])]['Percentage'].sum()\n", " female_pct = pop_dist_group_wp[pop_dist_group_wp['Sex'] == 'f']['Percentage'].sum()\n", " sex_ratio = (1 - female_pct) / female_pct * 100\n", " reproductive_age = pop_dist_group_wp[(pop_dist_group_wp['Sex'] == 'f') & (pop_dist_group_wp['Age_Bracket'].isin(['15-19', '20-24', '25-29', '30-34', '35-39', '40-44', '45-49']))]['Sexed_Percent'].sum()\n", "\n", " \n", " print(f\"under5: {under5:.2%}\")\n", " print(f\"youth (15-24): {youth:.2%}\")\n", " print(f\"working_age (15-64): {working_age:.2%}\")\n", " print(f\"elderly (60+): {elderly:.2%}\")\n", " print(f\"reproductive_age, percent of women (15-50): {reproductive_age:.2%}\")\n", " print(f\"sex_ratio: {round(sex_ratio, 2)} males to 100 females\")\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 72, "id": "a83777d1-ef48-4145-a841-ba3c934b99ef", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'\\ncity=\\'Mumbai\\'\\ndef oxford_age_stats():\\n with open(\"global_inputs.yml\", \\'r\\') as f:\\n global_inputs = yaml.safe_load(f)\\n\\n # Get city inputs\\n city_inputs = global_inputs.get(\\'city_inputs\\')\\n \\n# Define the \\'oxford_age_stats\\' function\\ndef oxford_age_stats(city):\\n with open(\"global_inputs.yml\", \\'r\\') as f:\\n global_inputs = yaml.safe_load(f)\\n\\n # Get city inputs\\n city_inputs = global_inputs.get(\\'city_inputs\\')\\n\\n if \\'oxford\\' in menu and menu[\\'oxford\\']: \\n oxford_data_path = os.path.join(multi_scan_folder, \"Oxford Global Cities Data.csv\")\\n if os.path.exists(oxford_data_path): \\n oxford_data = pd.read_csv(oxford_data_path)\\n indicators = oxford_data[\\'Indicator\\'].drop_duplicates()\\n pop_dist_inds = [indicator for indicator in indicators if \"Population\" in indicator and indicator not in [\"Population 0-14\", \"Population 15-64\", \"Population 65+\"]] \\n\\n if city in oxford_data[\\'Location\\'].values:\\n print(f\"{city} exists in the Oxford file.\")\\n pop_dist_structure = oxford_data.loc[(oxford_data[\\'Location\\'] == city) & (oxford_data[\\'Indicator\\'].isin(pop_dist_inds))]\\n print(pop_dist_structure.head(3))\\n pop_dist_structure[\\'Age_Bracket\\'] = pop_dist_structure[\\'Indicator\\'].str[11:19]\\n # Convert to numeric, handling errors by setting them to NaN\\n pop_dist_structure[\\'Age_Bracket\\'] = pd.to_numeric(pop_dist_structure[\\'Age_Bracket\\'], errors=\\'coerce\\')\\n pop_dist_structure[\\'Group\\'] = pd.cut(\\n pop_dist_structure[\\'Age_Bracket\\'],\\n bins=[0, 4, 14, np.inf], # Replace float(\\'inf\\') with np.inf\\n labels=[\\'Young\\', \\'Working\\', \\'65+\\']\\n )\\n \\n pop_dist_structure = pop_dist_structure.groupby([\\'Year\\', \\'Group\\']).agg(Count=(\\'Value\\', \\'sum\\')).reset_index()\\n pop_dist_structure[\\'Percent\\'] = pop_dist_structure.groupby(\\'Year\\')[\\'Count\\'].transform(lambda x: x / x.sum())\\n pop_dist_structure[\\'pct_sum\\'] = pop_dist_structure.groupby(\\'Year\\')[\\'Percent\\'].cumsum()\\n return pop_dist_structure\\n else:\\n print(f\"{city} does not exist in the Oxford file.\")\\n else:\\n print(\"Oxford file does not exist.\")\\n else:\\n print(\"Oxford is not selected in the menu.\")\\n\\n# Example usage\\noxford_age_stats(\\'Mumbai\\')\\n'" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "#Age structure Oxford Economics (If in Oxford)\n", "\n", "'''\n", "city='Mumbai'\n", "def oxford_age_stats():\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " # Get city inputs\n", " city_inputs = global_inputs.get('city_inputs')\n", " \n", "# Define the 'oxford_age_stats' function\n", "def oxford_age_stats(city):\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " # Get city inputs\n", " city_inputs = global_inputs.get('city_inputs')\n", "\n", " if 'oxford' in menu and menu['oxford']: \n", " oxford_data_path = os.path.join(multi_scan_folder, \"Oxford Global Cities Data.csv\")\n", " if os.path.exists(oxford_data_path): \n", " oxford_data = pd.read_csv(oxford_data_path)\n", " indicators = oxford_data['Indicator'].drop_duplicates()\n", " pop_dist_inds = [indicator for indicator in indicators if \"Population\" in indicator and indicator not in [\"Population 0-14\", \"Population 15-64\", \"Population 65+\"]] \n", "\n", " if city in oxford_data['Location'].values:\n", " print(f\"{city} exists in the Oxford file.\")\n", " pop_dist_structure = oxford_data.loc[(oxford_data['Location'] == city) & (oxford_data['Indicator'].isin(pop_dist_inds))]\n", " print(pop_dist_structure.head(3))\n", " pop_dist_structure['Age_Bracket'] = pop_dist_structure['Indicator'].str[11:19]\n", " # Convert to numeric, handling errors by setting them to NaN\n", " pop_dist_structure['Age_Bracket'] = pd.to_numeric(pop_dist_structure['Age_Bracket'], errors='coerce')\n", " pop_dist_structure['Group'] = pd.cut(\n", " pop_dist_structure['Age_Bracket'],\n", " bins=[0, 4, 14, np.inf], # Replace float('inf') with np.inf\n", " labels=['Young', 'Working', '65+']\n", " )\n", " ''''''\n", " pop_dist_structure = pop_dist_structure.groupby(['Year', 'Group']).agg(Count=('Value', 'sum')).reset_index()\n", " pop_dist_structure['Percent'] = pop_dist_structure.groupby('Year')['Count'].transform(lambda x: x / x.sum())\n", " pop_dist_structure['pct_sum'] = pop_dist_structure.groupby('Year')['Percent'].cumsum()\n", " return pop_dist_structure\n", " else:\n", " print(f\"{city} does not exist in the Oxford file.\")\n", " else:\n", " print(\"Oxford file does not exist.\")\n", " else:\n", " print(\"Oxford is not selected in the menu.\")\n", "\n", "# Example usage\n", "oxford_age_stats('Mumbai')\n", "'''\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "07361797-faf5-4edc-9374-3e416c3dff56", "metadata": {}, "outputs": [], "source": [ "# Share of GDP, Emp, Pop (If in Oxford)" ] }, { "cell_type": "code", "execution_count": 12, "id": "bccb56a3-73a4-4c63-ba6a-7554098455f0", "metadata": {}, "outputs": [], "source": [ "# GDP, Pop, Emp Growth (If in Oxford)" ] }, { "cell_type": "code", "execution_count": 13, "id": "4fbd4ebe-58ff-42c7-9bb1-d7bc692ebcfe", "metadata": {}, "outputs": [], "source": [ "# GDP per capita (If in Oxford)" ] }, { "cell_type": "code", "execution_count": 14, "id": "82a89901-b37b-4440-8a70-50c288bd3700", "metadata": {}, "outputs": [], "source": [ "# Share of employment by sector (If in Oxford)" ] }, { "cell_type": "code", "execution_count": 15, "id": "976d9c30-ee1e-42ee-b3cc-881d5830a540", "metadata": {}, "outputs": [], "source": [ "#Economic Inequality (If in Oxford)" ] }, { "cell_type": "code", "execution_count": 73, "id": "7ae3cd22-889c-4859-8555-c955df2d6bfd", "metadata": {}, "outputs": [], "source": [ "def wsf_stats():\n", " # Load global inputs from YAML file\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " # Get city inputs\n", " city_inputs = global_inputs.get('city_inputs')\n", "\n", " if menu['wsf']: \n", " wsf_stats_file = Path(output_folder / f\"{city}_wsf_stats.csv\") \n", " if os.path.exists(wsf_stats_file): \n", " pass\n", " else:\n", " print(\"WSF stats file does not exist.\")\n", " return \n", " wsf = pd.read_csv(wsf_stats_file)\n", " \n", " wsf = wsf.rename(columns={'year': 'Year'}).\\\n", " loc[:, ['Year', 'cumulative sq km']].\\\n", " rename(columns={'cumulative sq km': 'uba_km2'})\n", "\n", " \n", " wsf['growth_pct'] = (wsf['uba_km2'] / wsf['uba_km2'].shift(1) - 1)\n", " wsf['growth_km2'] = wsf['uba_km2'] - wsf['uba_km2'].shift(1)\n", "\n", " # Plot\n", " plt.figure(figsize=(8, 6))\n", " plt.plot(wsf['Year'], wsf['uba_km2'], marker='o', linestyle='-')\n", " plt.title(\"Urban Built-up Area, 1985-2015\")\n", " plt.xlabel(\"Year\")\n", " plt.ylabel(\"Urban built-up area (sq. km)\")\n", " plt.grid(True)\n", " # Save as PNG in output folder\n", " render_path = os.path.join(render_folder, \"urban_built_up_area.png\")\n", " plt.savefig(render_path)\n", " plt.close()\n", "\n", " # Create Plotly\n", " fig = go.Figure()\n", "\n", " # Add trace for urban built-up area\n", " fig.add_trace(go.Scatter(\n", " x=wsf['Year'],\n", " y=wsf['uba_km2'],\n", " mode='lines+markers',\n", " name='Urban built-up area (km^2)'\n", " ))\n", "\n", " \n", " fig.update_layout(\n", " title=\"Urban Built-up Area, 1985-2015\",\n", " xaxis_title=\"Year\",\n", " yaxis_title=\"Urban built-up area (km^2)\",\n", " template='plotly_white',\n", " showlegend=True,\n", " hovermode='x'\n", " )\n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", " first_area = wsf['uba_km2'].iloc[0]\n", " latest_area = wsf['uba_km2'].iloc[-1]\n", " first_year = wsf['Year'].iloc[0]\n", " latest_year = wsf['Year'].iloc[-1]\n", " pct_growth = 100 * (latest_area - first_area) / first_area\n", " print(f\"The city's built-up area grew from {round(first_area, 2)} sq. km in {first_year} to {round(latest_area, 2)} in {latest_year} for {round(pct_growth, 2)}% growth\")\n" ] }, { "cell_type": "code", "execution_count": 74, "id": "0b3913e5-3125-4d45-bca2-9b5558d82040", "metadata": {}, "outputs": [], "source": [ "#Landcover Graph\n", "def lc_stats():\n", " \n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " \n", " city_inputs = global_inputs.get('city_inputs')\n", "\n", " if menu['landcover']: \n", " lc_stats_file = Path(output_folder / f\"{city}_lc.csv\") \n", " if not lc_stats_file.exists(): \n", " print(\"Land Cover stats file does not exist.\")\n", " return \n", " \n", " \n", " lc = pd.read_csv(lc_stats_file)\n", "\n", " # Define colors\n", " lc_colors = {\n", " \"Tree cover\": \"#397e48\",\n", " \"Built-up\": \"#c4281b\",\n", " \"Grassland\": \"#88af52\",\n", " \"Bare / sparse vegetation\": \"#a59b8f\",\n", " \"Cropland\": \"#e49634\",\n", " \"Water bodies\": \"#429bdf\",\n", " \"Permanent water bodies\": \"#00008b\",#change this\n", " \"Mangroves\": \"#90EE90\",#change this\n", " \"Moss and lichen\": \"#013220\",#change this\n", " \"Shrubland\": \"#dfc25a\",\n", " \"Herbaceous wetland\": \"#7d87c4\",\n", " \"Snow and ice\": \"#F5F5F5\"\n", " }\n", " #connect it to the frontend layers.yml\n", "\n", " \n", " total_pixels = lc['Pixel Count'].sum()\n", "\n", " \n", " lc['Percentage'] = lc['Pixel Count'] / total_pixels * 100\n", "\n", " \n", " plt.figure(figsize=(8, 8))\n", " plt.pie(lc['Pixel Count'], labels=None, colors=[lc_colors[lc] for lc in lc['Land Cover Type']], startangle=90)\n", " plt.title(\"Land Cover Distribution\")\n", " \n", "\n", " \n", " render_path = os.path.join(render_folder, \"landcover.png\")\n", " plt.savefig(render_path)\n", " plt.close()\n", " \n", " \n", " fig = px.pie(lc, values='Pixel Count', names='Land Cover Type', color='Land Cover Type', color_discrete_map=lc_colors)\n", " fig.update_traces(textposition='inside', textinfo='percent+label')\n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", "\n", " # Dictionary to convert integers to ordinal strings\n", " ordinal_dict = {1: \"first\", 2: \"second\", 3: \"third\"}\n", "\n", " # Print the percentage of the first three highest values\n", " highest_values = lc.sort_values(by='Percentage', ascending=False).head(3)\n", " for i, (index, row) in enumerate(highest_values.iterrows()):\n", " ordinal_str = ordinal_dict.get(i + 1, str(i + 1) + \"th\")\n", " print(f\"The {ordinal_str} highest landcover value is {row['Land Cover Type']} with {row['Percentage']:.2f}% of the total land area\")\n", " \n", " \n", "\n" ] }, { "cell_type": "code", "execution_count": 75, "id": "41ad0554-e5ec-4ed7-b374-b9b3a48c6428", "metadata": {}, "outputs": [], "source": [ "#Elevation Graph\n", "def elev_stats():\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", " city_inputs = global_inputs.get('city_inputs')\n", " if menu['elevation']: \n", " elev_stats_file = Path(output_folder / f\"{city}_elevation.csv\") \n", " if not elev_stats_file.exists(): \n", " print(\"Elevation stats file does not exist.\")\n", " return \n", " elev = pd.read_csv(elev_stats_file)\n", " elev['percent'] = pd.to_numeric(elev['percent'], errors='coerce')\n", "\n", " \n", " elevation = elev.dropna(subset=['legend']).copy() # Create a copy to avoid the warning\n", "\n", " # Preprocess the data\n", " elevation['Percent'] = elevation['percent'].apply(lambda x: f\"{x:.0%}\") \n", " elevation['Elevation'] = pd.to_numeric(elevation['legend'].str.split('-').str[0], errors='coerce') \n", " elevation = elevation.dropna(subset=['Elevation']) \n", " elevation['legend'] = pd.Categorical(elevation['legend'], categories=elevation['legend'].unique()) \n", " print(elevation)\n", " # Define colors\n", " elevation_colors = {\n", " \"0-2\": \"#f5c4c0\",\n", " \"2-5\": \"#f19bb4\",\n", " \"5-10\": \"#ec5fa1\",\n", " \"10-20\": \"#c20b8a\",\n", " \"20+\": \"#762175\"\n", " }\n", " # Create the pie chart\n", " plt.figure(figsize=(8, 5))\n", " wedges, texts, autotexts = plt.pie(elevation['percent'], labels=None, colors=elevation_colors.values(), startangle=90, autopct='%1.0f%%')\n", " plt.title(\"Elevation Distribution\")\n", " plt.legend(elevation['Elevation'], loc=\"center left\", bbox_to_anchor=(1, 0, 0.5, 1))\n", "\n", " # Add percentage labels\n", " for autotext in autotexts:\n", " autotext.set_color('white')\n", " autotext.set_size(10)\n", " \n", " # Save the plot\n", " render_path = os.path.join(render_folder, \"elevation.png\")\n", " plt.savefig(render_path, bbox_inches='tight')\n", " plt.close()\n", "\n", " # Create a pie chart with percentage labels\n", " fig = px.pie(elevation, values='percent', names='Elevation', hole=0.3,\n", " color_discrete_sequence=list(elevation_colors.values()), \n", " title=\"Elevation Distribution\", \n", " labels={'percent': 'Percentage', 'Elevation': 'Elevation Range'},\n", " template='plotly_white')\n", "\n", " # Update layout\n", " fig.update_layout(legend=dict(orientation=\"h\", x=0.5, y=1.1),\n", " margin=dict(l=0, r=0, t=50, b=0))\n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", "\n", " # Print the percentage of land at an elevation of 20+\n", " max_percent_index = elevation['percent'].idxmax()\n", " highest_percent_row = elevation.loc[max_percent_index]\n", " print(\"Highest percentage entry for Elevation:\")\n", " print(highest_percent_row)\n", "\n", "\n", "\n", "\n", "\n", "#why are they different?" ] }, { "cell_type": "code", "execution_count": 76, "id": "be08dae0-5584-42e6-b59d-9eaa5cd4d6a7", "metadata": {}, "outputs": [], "source": [ "#Slope Graph\n", "def slope_stats():\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", " city_inputs = global_inputs.get('city_inputs')\n", " if menu['elevation']: \n", " slope_stats_file = Path(output_folder / f\"{city}_slope.csv\") \n", " if not slope_stats_file.exists(): \n", " print(\"Elevation stats file does not exist.\")\n", " return \n", " slope = pd.read_csv(slope_stats_file)\n", " slope['percent'] = pd.to_numeric(slope['percent'], errors='coerce')\n", "\n", " # Drop rows with NaN values in 'legend'\n", " slope = slope.dropna(subset=['legend']).copy()\n", "\n", " # Preprocess the data\n", " slope['Percent'] = slope['percent'].apply(lambda x: f\"{x:.0%}\")\n", " slope['Slope'] = slope['legend'].str.extract(r\"(\\d+)\").astype(float)\n", " slope['legend'] = pd.Categorical(slope['legend'], categories=slope['legend'].unique())\n", "\n", " # Define colors\n", " slope_colors = {\n", " \"0-2\": \"#ffffd4\",\n", " \"2-5\": \"#fed98e\",\n", " \"5-10\": \"#fe9929\",\n", " \"10-20\": \"#d95f0e\",\n", " \"20+\": \"#993404\"\n", " }\n", "\n", " # Plot the donut chart\n", " plt.figure(figsize=(8, 5))\n", " plt.pie(slope['percent'], labels=None, colors=slope_colors.values(), startangle=90)\n", " plt.title(\"Slope Distribution\")\n", " plt.legend(slope['Slope'], loc=\"center left\", bbox_to_anchor=(1, 0, 0.5, 1))\n", " # Save the plot\n", " render_path = os.path.join(render_folder, \"slope.png\")\n", " plt.savefig(render_path, bbox_inches='tight')\n", " plt.close()\n", "\n", " # Plotly\n", " fig = px.pie(slope, values='percent', names='Slope', hole=0.3,\n", " color_discrete_sequence=list(slope_colors.values()), \n", " title=\"Slope Distribution\", \n", " labels={'percent': 'Percentage', 'Slope': 'Slope Range'},\n", " template='plotly_white')\n", "\n", " # Update layout\n", " fig.update_layout(legend=dict(orientation=\"h\", x=0.5, y=1.1),\n", " margin=dict(l=0, r=0, t=50, b=0))\n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", "\n", " # Print the highest percentage value and consequent class\n", " max_percent_index = slope['percent'].idxmax()\n", " highest_percent_row = slope.loc[max_percent_index]\n", " print(\"Highest percentage entry for Slope:\")\n", " print(highest_percent_row)\n" ] }, { "cell_type": "code", "execution_count": 31, "id": "ab2851a1-315a-445d-8982-bec89bbe1bc2", "metadata": {}, "outputs": [], "source": [ "#Cyclone Graph" ] }, { "cell_type": "code", "execution_count": null, "id": "505a5d75", "metadata": {}, "outputs": [], "source": [ "#earthquake timeline\n" ] }, { "cell_type": "code", "execution_count": 77, "id": "e56dba45", "metadata": {}, "outputs": [], "source": [ "import yaml\n", "import geopandas as gpd\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from shapely.geometry import Point, LineString\n", "\n", "def flood_timeline():\n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " city_inputs = global_inputs.get('city_inputs')\n", "\n", " if menu['flood']: \n", " flood_archive_path = global_inputs.get('flood_timeline_source')\n", " flood_archive = gpd.read_file(flood_archive_path)\n", " flood_archive = flood_archive[flood_archive.is_valid]\n", " aoi = features.to_crs(flood_archive.crs)\n", " flood_archive = flood_archive[flood_archive.intersects(aoi.unary_union)]\n", " fig, ax = plt.subplots()\n", " flood_archive.plot(column='DEAD', ax=ax, legend=True)\n", "\n", " floods = flood_archive[['BEGAN', 'ENDED', 'DEAD', 'DISPLACED', 'MAINCAUSE', 'SEVERITY']]\n", "\n", " # Tally of flood events\n", " print(\"Tally of flood events\")\n", " print(floods.agg({'DEAD': 'sum', 'DISPLACED': 'sum', 'BEGAN': 'count'}))\n", "\n", " duration = (pd.to_datetime(floods['ENDED']) - pd.to_datetime(floods['BEGAN'])).dt.days\n", "\n", " # Prepare text for plotting\n", " flood_text = floods.copy()\n", " flood_text['severity'] = np.select([flood_text['SEVERITY'] == 1, flood_text['SEVERITY'] == 1.5, flood_text['SEVERITY'] == 2],\n", " ['Large event', 'Very large event', 'Extreme event'])\n", " flood_text['duration'] = duration\n", " flood_text['text'] = flood_text.apply(lambda row: f\"{row['BEGAN']}, {row['severity'].lower()} {row['MAINCAUSE']}, \"\n", " f\"{row['duration']} days, {row['DEAD']:,} fatalities, {row['DISPLACED']:,} displaced\",\n", " axis=1)\n", "\n", " # Calculate coordinates for text plotting\n", " flood_text['node_x'] = pd.to_datetime(flood_text['BEGAN']) + pd.to_timedelta(1460 * (2 * (flood_text.index % 2) - 1), unit='D')\n", " flood_text['node_y'] = 0\n", "\n", " # Prepare lines for text plotting\n", " flood_lines = pd.concat([flood_text[['BEGAN', 'node_x', 'node_y']],\n", " pd.DataFrame({'BEGAN': flood_text['BEGAN'], 'node_x': flood_text['BEGAN'], 'node_y': 0}),\n", " pd.DataFrame({'BEGAN': flood_text['BEGAN'], 'node_x': flood_text['BEGAN'], 'node_y': 0})])\n", "\n", " # Plotting\n", " fig, ax = plt.subplots(figsize=(20, 5.833))\n", " ax.plot(flood_lines['node_x'], flood_lines['node_y'], color='blue', linestyle='-')\n", " ax.scatter(flood_text['node_x'], flood_text['node_y'], color='red')\n", " for idx, row in flood_text.iterrows():\n", " ax.text(row['node_x'] + pd.to_timedelta(30, unit='D'), row['node_y'] - 50, row['text'], rotation=30, ha='center')\n", " ax.set_xlim(pd.Timestamp('1984-01-01'), pd.Timestamp('2020-12-31'))\n", " ax.set_ylim(-1800, 1800)\n", " ax.set_yticks([1825, 3650])\n", " ax.set_yticklabels([1825, 3650])\n", " plt.show()\n", "\n", "\n", "\n", "#workshop this one \n" ] }, { "cell_type": "code", "execution_count": 78, "id": "5e8fe622-51c5-4159-906f-74405ebe4101", "metadata": {}, "outputs": [], "source": [ "def extract_monthly_stats():\n", " \n", " with open(\"global_inputs.yml\", 'r') as f:\n", " global_inputs = yaml.safe_load(f)\n", "\n", " pv_directory = global_inputs.get('solar_graph_source')\n", "\n", " \n", " pv_files = [f for f in os.listdir(pv_directory) if f.endswith('.tif')]\n", "\n", " \n", " monthly_pv = []\n", "\n", " aoi = aoi_file.geometry\n", "\n", " for f in pv_files:\n", " pattern = re.compile(r'PVOUT_(\\d{2})')\n", " match = pattern.search(f)\n", " if match:\n", " m = int(match.group(1)) # Extract the month\n", "\n", " file_path = os.path.join(pv_directory, f)\n", "\n", " with rasterio.open(file_path) as src:\n", " \n", " raster_data, raster_transform = mask(src, aoi.geometry, crop=True)\n", "\n", " stats = {\n", " 'month': m,\n", " 'max': np.nanmax(raster_data),\n", " 'min': np.nanmin(raster_data),\n", " 'mean': np.nanmean(raster_data),\n", " 'sum': np.nansum(raster_data)\n", " }\n", "\n", " \n", " monthly_pv.append(stats)\n", " else:\n", " print(f\"No match found for filename: {f}\")\n", "\n", " monthly_pv_df = pd.DataFrame(monthly_pv)\n", "\n", " # Sort by 'month'\n", " monthly_pv_df.sort_values(by='month', inplace=True)\n", "\n", " # Calculate daily PV energy yield\n", " monthly_pv_df['daily'] = monthly_pv_df['mean'] \n", "\n", " \n", " highest_value = monthly_pv_df['daily'].max()\n", " lowest_value = monthly_pv_df['daily'].min()\n", "\n", " \n", " ratio = highest_value / lowest_value\n", "\n", " # Check if the ratio is greater than 2.5\n", " if ratio > 2.5:\n", " print(\"Seasonality is high, making solar energy available throughout the year\")\n", " else:\n", " print(\"Seasonality is low to moderate, making solar energy available in only some of the months\")\n", "\n", " # Plot\n", " plt.figure(figsize=(10, 6))\n", " plt.plot(monthly_pv_df['month'], monthly_pv_df['daily'], marker='o', linestyle='-')\n", " plt.text(1, 4.6, 'Excellent Conditions', color='darkgrey', verticalalignment='bottom', horizontalalignment='left')\n", " plt.text(1, 3.6, 'Favorable Conditions', color='darkgrey', verticalalignment='bottom', horizontalalignment='left')\n", " plt.xlabel('Month')\n", " plt.ylabel('Daily PV energy yield (kWh/kWp)')\n", " plt.title('Seasonal availability of solar energy')\n", " plt.xticks(np.arange(1, 13), ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])\n", " plt.grid(True)\n", " plt.axhline(y=3.5, linestyle='--', color='black')\n", " plt.axhline(y=4.5, linestyle='--', color='black')\n", " plt.tight_layout()\n", " # Save the plot\n", " render_path = os.path.join(render_folder, \"_PV_graph.png\")\n", " plt.savefig(render_path, bbox_inches='tight')\n", " plt.close()\n", " #Plotly\n", " fig = px.line(monthly_pv_df, x='month', y='daily', markers=True)\n", " fig.add_annotation(x=1, y=4.6, text='Excellent Conditions', showarrow=False, font=dict(color='darkgrey'), xshift=5)\n", " fig.add_annotation(x=1, y=3.6, text='Favorable Conditions', showarrow=False, font=dict(color='darkgrey'), xshift=5)\n", " fig.update_xaxes(title='Month', tickvals=list(range(1, 13)), ticktext=['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])\n", " fig.update_yaxes(title='Daily PV energy yield (kWh/kWp)', range=[0, 5])\n", " fig.add_shape(type=\"line\", x0=1, y0=3.5, x1=12, y1=3.5, line=dict(color=\"black\", width=1, dash='dash'))\n", " fig.add_shape(type=\"line\", x0=1, y0=4.5, x1=12, y1=4.5, line=dict(color=\"black\", width=1, dash='dash'))\n", " fig.update_layout(title='Seasonal availability of solar energy', xaxis=dict(showgrid=True, zeroline=False),plot_bgcolor='white')\n", "\n", " fig.show()\n", " fig.write_html(render_path.replace('.png', '.html'))\n", "\n" ] }, { "cell_type": "code", "execution_count": 79, "id": "5e48c01f-f509-4137-8d5b-78a7f336dbf6", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/ipshitakarmakar/mambaforge/envs/geo/share/jupyter/nbconvert/templates/base/display_priority.j2:32: UserWarning:\n", "\n", "Your element with mimetype(s) dict_keys(['application/vnd.plotly.v1+json']) is not able to be represented.\n", "\n" ] } ], "source": [ "#Save output\n", "def export_outputs_to_markdown(notebook_path, output_path):\n", " # Load the notebook\n", " with open(notebook_path, 'r', encoding='utf-8') as f:\n", " notebook_content = nbformat.read(f, as_version=4)\n", " \n", " # Initialize the Markdown exporter\n", " markdown_exporter = MarkdownExporter()\n", " markdown_exporter.exclude_input = True # Exclude input cells from the Markdown\n", " \n", " # Convert the notebook to Markdown format\n", " markdown_output, resources = markdown_exporter.from_notebook_node(notebook_content)\n", " \n", " # Remove the folders for plots and pickles from the resources\n", " resources.pop('outputs', None)\n", " resources.pop('output_files', None)\n", " \n", " # Write the Markdown content to a file\n", " with open(output_path, 'w', encoding='utf-8') as f:\n", " f.write(markdown_output)\n", "\n", "# Path to the input notebook\n", "input_notebook_path = \"/Users/ipshitakarmakar/Documents/GitHub/city-scan-automation/scan_assembly.ipynb\"\n", "\n", "# Path for the output Markdown file\n", "output_markdown_path = \"/Users/ipshitakarmakar/Documents/GitHub/city-scan-automation/output_notebook.md\"\n", " \n", "# Export the outputs to Markdown\n", "export_outputs_to_markdown(input_notebook_path, output_markdown_path)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 5 }