Heat-scan-workflow-city-level.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><center> <font color='cyan'> Urban Heat Scan Workflow </font> </center></h2>\n",
    "<h3><center> <font color='cyan'>Table of Contents</font>   </center></h3>\n",
    "\n",
    "[Preface: Import dependencies and set paths](#section0)\n",
    "<br>\n",
    "[Section 1: Big picture: City map and context](#section1)\n",
    "<br>\n",
    "[Section 2: Population density and tree cover city level](#section2) \n",
    "<br>\n",
    "[Section 3: Visualize Vito heat rasters at city level](#section3) \n",
    "<br>\n",
    "[Section 4: Create the Heat Scan PPT](#section4) \n",
    "<br>\n",
    "[Section 5: Population density and tree cover Sub-city level. NA ](#section5) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section0'></a>\n",
    "<h5><center> <font color='cyan'> Preface: Import dependencies and set paths</font>   </center></h5>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 328,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Aziz\\AppData\\Local\\Temp\\ipykernel_12132\\783695111.py:13: UserWarning: The `utils.config` function is deprecated and will be removed in a future release. Instead, use the `settings` module directly to configure a global setting's value. For example, `ox.settings.log_console=True`.\n",
      "  ox.config(use_cache=True, log_console=True)\n"
     ]
    }
   ],
   "source": [
    "# Import packages\n",
    "import geopandas as gpd , os , time , math, glob\n",
    "from collections import defaultdict\n",
    "\n",
    "from urllib.request import urlopen\n",
    "from pyproj import CRS\n",
    "from os.path import exists\n",
    "from geopandas.tools import overlay\n",
    "from rasterio.warp import calculate_default_transform, reproject, Resampling\n",
    "from osgeo import ogr\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "ox.config(use_cache=True, log_console=True)\n",
    "\n",
    "import  seaborn as sn\n",
    "import geemap \n",
    "import ee \n",
    "import rasterio\n",
    "from rasterio.mask import mask\n",
    "from rasterio.features import shapes\n",
    "import osmnx as ox\n",
    "from shapely.geometry import box\n",
    "from rasterio.plot import show\n",
    "import pylab as plt\n",
    "from geopandas.tools import sjoin\n",
    "import contextily as cx\n",
    "from h3 import h3\n",
    "import h3pandas\n",
    "from rasterstats import zonal_stats\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from rasterio.crs import CRS\n",
    "import matplotlib.colors as colors\n",
    "from matplotlib.lines import Line2D\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "InteractiveShell.ast_node_interactivity = \"last_expr\" #Supress matplotlib output\n",
    "\n",
    "# Stats heavy-lifting\n",
    "from esda.moran import Moran\n",
    "from libpysal.weights import Queen, KNN\n",
    "import  libpysal\n",
    "\n",
    "from sklearn import metrics\n",
    "from sklearn.preprocessing import robust_scale\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.cluster import AgglomerativeClustering\n",
    "from pysal.lib import weights\n",
    "from pysal.explore import esda\n",
    "from pysal.model import spreg\n",
    "\n",
    "from pptx import Presentation\n",
    "from pptx.util import Inches\n",
    "from pptx.util import Inches, Pt\n",
    "from pptx.enum.text import PP_ALIGN\n",
    "\n",
    "\n",
    "from pptx.oxml.xmlchemy import OxmlElement\n",
    "from pptx.util import Cm\n",
    "from pptx.enum.shapes import MSO_SHAPE\n",
    "from pptx.dml.color import RGBColor\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "plt.rcParams['hatch.color'] = 'white'\n",
    "plt.rcParams['hatch.linewidth'] = 4.2\n",
    "plt.rcParams['font.family'] = 'Yu Gothic'\n",
    "\n",
    "# plt.rcParams['figure.figsize'] = 7.5 , 7.5\n",
    "width=10\n",
    "height=7.6\n",
    "plt.rcParams['figure.figsize'] = width , height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_paths(base_dir):\n",
    "    data = os.path.join(base_dir, 'data') \n",
    "    output = os.path.join(base_dir, 'output') \n",
    "    shapefiles = os.path.join(output, 'shapefiles') \n",
    "    maps = os.path.join(output, 'maps') \n",
    "    rasters = os.path.join(output, 'rasters') \n",
    "    tables = os.path.join(output, 'tables') \n",
    "    \n",
    "    dirs_list= [data, output, base_dir, shapefiles , maps , rasters , tables]\n",
    "    for dir in dirs_list:\n",
    "        if not os.path.exists(dir):\n",
    "            os.mkdir(dir)\n",
    "\n",
    "    return data, shapefiles , maps , rasters , output ,tables\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Initialize ee and import shapefile and select city\n",
    "ee.Initialize()\n",
    "run_vito_analysis=True #\n",
    "run_tree_and_pop_analysis=True  #False\n",
    "create_the_ppt=True\n",
    "WGS84=4326 #uses degrees\n",
    "WGS84_meters=3857 #uses meters\n",
    "EPSG_str= 'EPSG:4326'\n",
    "crs = 4326\n",
    "\n",
    "city= \"Nis City\"\n",
    "\n",
    "base_dir = f'C:/Users/Aziz/Dropbox/CRP/UHI/{city}'\n",
    "data, shapefiles , maps , rasters , output ,tables = set_paths(base_dir)\n",
    "\n",
    "# shapefile_path=f'{shapefiles}/geoBoundaries-SRB-ADM1_simplified.shp'\n",
    "\n",
    "# Sub city donot exist somehow. City and municipailty are the same\n",
    "# Serbia is divided into 145 municipalities and 29 cities,[2] which form the basic units of local government. \n",
    "# Each municipality has its own assembly (elected every four years in local elections), a municipal president, \n",
    "# public service property and a budget. Municipalities usually have more than 10,000 inhabitants.[2]\n",
    "\n",
    "shapefile_path=f'{shapefiles}/geoBoundaries-SRB-ADM2.shp'\n",
    "gdf_city = gpd.read_file(shapefile_path).to_crs(WGS84)\n",
    "gdf_city=gdf_city[gdf_city[\"shapeName\"] == city]\n",
    "country= gdf_city.iloc[0]['shapeGroup']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section1'></a>\n",
    "<h5><center> <font color='cyan'> Section 1: Big picture: City map</font>   </center></h5>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "def map_city(maps, vector_file, crs, legend_title,  visualize_column, title, map_output):\n",
    "    ax = vector_file.to_crs(crs).plot(figsize=(10, 10), \\\n",
    "                                   column= visualize_column, \\\n",
    "                                   alpha=0.6,  \\\n",
    "                                   facecolor='none',edgecolor='yellow' , \\\n",
    "                                   linewidth=3 , \n",
    "                                   legend=True,  \n",
    "                                   legend_kwds={'loc':'upper right', \n",
    "                                                'bbox_to_anchor':(1, 1), \n",
    "                                                'markerscale':1.01, \n",
    "                                                'title_fontsize':'small', \n",
    "                                                'fontsize':'x-small'\n",
    "                                                }  \n",
    "                                        )\n",
    "\n",
    "    ax.set_aspect('equal')\n",
    "    cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery,  crs=crs) \n",
    "    cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels,   crs=crs, zoom=11) # zoom=13\n",
    "\n",
    "    vector_file_buffer = vector_file.to_crs(crs).buffer(0.016) #for zoom out\n",
    "    minx, miny, maxx, maxy = vector_file_buffer.total_bounds\n",
    "    ax.set_xlim(minx, maxx)\n",
    "    ax.set_ylim(miny, maxy)\n",
    "    # Legends\n",
    "    LegendElement = [\n",
    "                    Line2D([0],[0],color='yellow',lw=3,label=f'{legend_title}')\n",
    "                    ]\n",
    "    ax.legend(handles=LegendElement,loc='upper right')\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])   \n",
    "    ax.title.set_text(f'{title}')\n",
    "    plt.tight_layout()\n",
    "    ax.figure.savefig(map_output,  bbox_inches='tight',   dpi = 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section2'></a>\n",
    "<h5><center> <font color='cyan'> Section 2: Population density and tree cover city level</font>   </center></h5>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_fc_collection(gdf_city):\n",
    "    poly = gdf_city.geometry.unary_union\n",
    "    gdf_boundary = gpd.GeoDataFrame(geometry=[poly],crs=gdf_city.crs)\n",
    "    geom = gdf_boundary['geometry']\n",
    "    jsonDict = eval(geom.to_json())\n",
    "    for index, row in gdf_boundary.iterrows(): \n",
    "        polygon_list= []\n",
    "        for x in jsonDict['features'][index]['geometry']['coordinates']:\n",
    "            polygon_list.append(x)\n",
    "            region = ee.Geometry.Polygon(polygon_list)\n",
    "            fc_filtered = ee.FeatureCollection(region)\n",
    "    # gdf_boundary.plot()\n",
    "    return fc_filtered  , region \n",
    "\n",
    "\n",
    "def clipToCol(image):\n",
    "  \"\"\"clip gee collection\n",
    "  args:\n",
    "      image: Image collection\n",
    "  returns:\n",
    "    ee-collection: Clipped image collection\n",
    "   \n",
    "   \"\"\"\n",
    "  return image.clip(fc_filtered)\n",
    "\n",
    "def download_tree_cover(ee , EPSG_str, region):\n",
    "    tree_cover=ee.ImageCollection(\"projects/sat-io/open-datasets/GFCC30TC\") \\\n",
    "        .limit(1, 'system:time_start', False).first().clip(fc_filtered)\n",
    "    tree_cover_tif = os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "    geemap.ee_export_image(\n",
    "        tree_cover, filename=tree_cover_tif,  \n",
    "        crs=EPSG_str,\n",
    "        # crs_transform=crs_transform,\n",
    "        scale=30, region=region, file_per_band=False\n",
    "    )\n",
    "    return tree_cover_tif\n",
    "\n",
    "def download_tree_cover_esa_vito(ee , EPSG_str, region):\n",
    "    tree_cover=ee.ImageCollection(\"ESA/WorldCover/v100\") \\\n",
    "        .limit(1, 'system:time_start', False).first().clip(fc_filtered).eq(10) #trees only\n",
    "        # .limit(1, 'system:time_start', False).select(10).first().clip(fc_filtered) #trees only\n",
    "    tree_cover_tif_esa = os.path.join(rasters, 'tree_cover_projected_esa.tif')\n",
    "    geemap.ee_export_image(\n",
    "        tree_cover, filename=tree_cover_tif_esa,  \n",
    "        crs=EPSG_str,\n",
    "        # crs_transform=crs_transform,\n",
    "        scale=10, region=region, file_per_band=False\n",
    "    )\n",
    "    return tree_cover_tif_esa\n",
    "\n",
    "\n",
    "def dowload_pop_density(ee , EPSG_str, region):\n",
    "    #  Define WorldPop & clip using function \n",
    "    pop_density = ee.ImageCollection(\"WorldPop/GP/100m/pop\").map(clipToCol). \\\n",
    "                    filterDate('2020').select('population').mosaic()\n",
    "\n",
    "    pop_density_tif = os.path.join(rasters, 'pop_density_projected.tif')\n",
    "    geemap.ee_export_image( pop_density, \n",
    "                           filename=pop_density_tif,  \n",
    "                           crs=EPSG_str,\n",
    "                           # crs_transform=crs_transform,\n",
    "                           scale=100, \n",
    "                           region=region, \n",
    "                           file_per_band=False\n",
    "                           )\n",
    "\n",
    "    return pop_density_tif\n",
    "\n",
    "# Takes longer\n",
    "def download_google_dynamic_worldcover(ee , EPSG_str, region):\n",
    "    # dynamic_world_cover=ee.ImageCollection(\"GOOGLE/DYNAMICWORLD/V1\") \\\n",
    "    #     .limit(1, 'system:time_start', False).map(clipToCol).mosaic() #all only\n",
    "    dynamic_world_cover=ee.ImageCollection(\"GOOGLE/DYNAMICWORLD/V1\") \\\n",
    "        .limit(1, 'system:time_start', False).mosaic().clip(fc_filtered)\n",
    "    google_dynamic_world_cover_projected = os.path.join(rasters, 'google_dynamic_world_cover_projected.tif')\n",
    "    geemap.ee_export_image(\n",
    "        dynamic_world_cover, filename=google_dynamic_world_cover_projected,  \n",
    "        crs=EPSG_str,\n",
    "        # crs_transform=crs_transform,\n",
    "        scale=10, region=region, file_per_band=False\n",
    "    )\n",
    "    return google_dynamic_world_cover_projected\n",
    "\n",
    "def download_esri_dynamic_worldcover(ee , EPSG_str, region): \n",
    "    dynamic_world_cover=ee.ImageCollection(\"projects/sat-io/open-datasets/landcover/ESRI_Global-LULC_10m_TS\")\\\n",
    "        .filterDate('2022-01-01','2022-12-31').mosaic().clip(fc_filtered)\n",
    "    esri_dynamic_world_cover_projected = os.path.join(rasters, 'esri_dynamic_world_cover_projected.tif')\n",
    "    geemap.ee_export_image(\n",
    "        dynamic_world_cover, filename=esri_dynamic_world_cover_projected,  \n",
    "        crs=EPSG_str,\n",
    "        # crs_transform=crs_transform,\n",
    "        scale=10, region=region, file_per_band=False\n",
    "    )\n",
    "    return esri_dynamic_world_cover_projected\n",
    "\n",
    "\n",
    "def reproject_rasters(crs ,output_dir, unprojected_raster, projected_raster):\n",
    "    if not exists(projected_raster):\n",
    "        with rasterio.open(unprojected_raster) as src:\n",
    "            dst_crs = 'EPSG:' + str(crs)\n",
    "            transform, width, height = calculate_default_transform(\n",
    "                src.crs, dst_crs, src.width, src.height, *src.bounds)\n",
    "            kwargs = src.meta.copy()\n",
    "            kwargs.update({\n",
    "                'crs': dst_crs,\n",
    "                'transform': transform,\n",
    "                'width': width,\n",
    "                'height': height\n",
    "                # ,'dtype': np.float32\n",
    "            })\n",
    "\n",
    "            with rasterio.open(projected_raster, 'w', **kwargs) as dst:\n",
    "                for i in range(1, src.count + 1):\n",
    "                    reproject(\n",
    "                        source=rasterio.band(src, i),\n",
    "                        destination=rasterio.band(dst, i),\n",
    "                        src_transform=src.transform,\n",
    "                        src_crs=src.crs,\n",
    "                        dst_transform=transform,\n",
    "                        dst_crs=dst_crs,\n",
    "                        # dtype= np.float32,\n",
    "                        resampling=Resampling.nearest)\n",
    "\n",
    "def polygonize(raster_file ,output_shapefile, crs, mask_value, raster_val):\n",
    "    # mask = None\n",
    "    with rasterio.Env():\n",
    "        with rasterio.open(raster_file ) as src: #, dtype= np.float32\n",
    "            image = src.read(1) # first band\n",
    "            no_data_value=0\n",
    "            image = np.where(image<0,no_data_value,image) #replaced negative values with 0\n",
    "            results = (\n",
    "            {'properties': {raster_val : v}, 'geometry': s}\n",
    "            for i, (s, v) \n",
    "            in enumerate(\n",
    "                shapes(image, mask=mask_value, transform=src.transform)))\n",
    "\n",
    "    # The result is a generator of GeoJSON features\n",
    "    geoms = list(results)\n",
    "    # first feature\n",
    "    # print(geoms[0])\n",
    "    gpd_polygonized_raster  = gpd.GeoDataFrame.from_features(geoms)\n",
    "    crs_utm = CRS.from_user_input(crs)\n",
    "    driver= 'ESRI Shapefile'\n",
    "    gpd_polygonized_raster.to_file(output_shapefile, driver, crs=crs_utm)\n",
    "    return gpd_polygonized_raster\n",
    "\n",
    "\n",
    "def calcBearing (lat1, long1, lat2, long2):\n",
    "    dLon = (long2 - long1)\n",
    "    x = math.cos(math.radians(lat2)) * math.sin(math.radians(dLon))\n",
    "    y = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.cos(math.radians(dLon))\n",
    "    bearing = math.atan2(x,y)   # use atan2 to determine the quadrant\n",
    "    bearing = math.degrees(bearing)\n",
    "    return bearing\n",
    "\n",
    "def calcNSEW(lat1, long1, lat2, long2):\n",
    "    points = [\"north\", \"north east\", \"east\", \"south east\", \"south\", \"south west\", \"west\", \"north west\"]\n",
    "    bearing = calcBearing(lat1, long1, lat2, long2)\n",
    "    bearing += 22.5\n",
    "    bearing = bearing % 360\n",
    "    bearing = int(bearing / 45) # values 0 to 7\n",
    "    NSEW = points [bearing]\n",
    "    return NSEW\n",
    "\n",
    "def ordinal(n: int):\n",
    "    if 11 <= (n % 100) <= 13:\n",
    "        suffix = 'th'\n",
    "    else:\n",
    "        suffix = ['th', 'st', 'nd', 'rd', 'th'][min(n % 10, 4)]\n",
    "    return str(n) + suffix\n",
    "\n",
    "def get_osm_tags_inside_clusters(cluster_boundaries):   \n",
    "    for index, poi in cluster_boundaries.iterrows():\n",
    "        polygon = cluster_boundaries.iloc[index]['geometry']\n",
    "        park , education, industry, building=get_services_infrastructure_inside_poly(polygon)\n",
    "    return park , education, industry, building\n",
    "\n",
    "def delete_shapefile(output_shapefile):\n",
    "    shapefile = ogr.Open(output_shapefile,1 )\n",
    "    layer=shapefile.GetLayerByIndex(0)\n",
    "    count=layer.GetFeatureCount()\n",
    "    for feature in range(count):\n",
    "        layer.DeleteFeature(feature)\n",
    "\n",
    "\n",
    "def map_clusters(maps, hexagons, legend_title,legends_format, crs ,cmap,  visualize_column, title):\n",
    "    hexagons=hexagons.to_crs(crs)\n",
    "    hexagons = hexagons.sort_values('cluster_area_deg', ascending = False).groupby('ward5wknn').head(2)\n",
    "    hexagons=hexagons.reset_index()\n",
    "    for index, poi in hexagons.iterrows():\n",
    "        try:\n",
    "            cluster_class=poi[\"ward5wknn\"]\n",
    "            polygon = hexagons[hexagons[\"ward5wknn\"]==index]\n",
    "            # polygon = hexagons.iloc[index]['geometry']  \n",
    "            map_output= f'{maps}/{city}_{visualize_column}_{index}_{cluster_class}_cluster.png'\n",
    "            ax = polygon.plot(figsize=(10, 10), \\\n",
    "                                        column= visualize_column, \\\n",
    "                                        facecolor='none',edgecolor='cyan' , \\\n",
    "                                        linewidth=1.5 , \\\n",
    "                                        alpha=0.6, categorical=True , cmap=cmap, \\\n",
    "                                        legend=True,  # Add legend \n",
    "                                        legend_kwds={'loc':'upper right', \n",
    "                                                        'bbox_to_anchor':(1, 1), \n",
    "                                                        'fmt':legends_format,\n",
    "                                                        'markerscale':1.01, \n",
    "                                                        'title_fontsize':'small', \n",
    "                                                        'fontsize':'x-small'\n",
    "                                                        # title_fontsizeint or {'xx-small', 'x-small', 'small'\n",
    "                                                        } ,\n",
    "                                        aspect=1\n",
    "                                                )\n",
    "            # vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "            polygon_buffer = polygon.to_crs(crs).buffer(0.016) #for zoom out\n",
    "            minx, miny, maxx, maxy = polygon_buffer.total_bounds\n",
    "            ax.set_xlim(minx, maxx)\n",
    "            ax.set_ylim(miny, maxy)\n",
    "\n",
    "            cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery,  crs=crs) \n",
    "            # vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "            # cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite,  crs=crs , zoom=12) \n",
    "            cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels,   crs=crs, zoom=11) # zoom=13\n",
    "\n",
    "            # plt.title(f'{title}',fontsize=16)\n",
    "            plt.xticks([])\n",
    "            plt.yticks([]) \n",
    "            leg1 = ax.get_legend()\n",
    "            # Set markers to square shape\n",
    "            for ea in leg1.legendHandles:\n",
    "                ea.set_marker('s')\n",
    "            leg1.set_title(f'{legend_title}')\n",
    "            ax.title.set_text(f'{title}')\n",
    "            plt.tight_layout()\n",
    "            # ax.figure.savefig(map_output)\n",
    "            ax.figure.savefig(map_output,  bbox_inches='tight' ,   dpi = 400)\n",
    "\n",
    "        except:\n",
    "            pass\n",
    "        \n",
    "def map_clusters_land_cover(maps, hexagons, legend_title,legends_format, crs ,cmap,  visualize_column, title):\n",
    "    hexagons=hexagons.to_crs(crs)\n",
    "    hexagons=hexagons.reset_index()\n",
    "    hexagons = hexagons.sort_values('cluster_area_deg', ascending = False).groupby('ward5wknn').head(2)\n",
    "    for index, poi in hexagons.iterrows():\n",
    "        try:\n",
    "            cluster_class=poi[\"ward5wknn\"]\n",
    "            polygon = hexagons[hexagons[\"ward5wknn\"]==index]\n",
    "            # polygon = hexagons.iloc[index]['geometry']  \n",
    "            map_output= f'{maps}/{city}_{visualize_column}_{index}_{cluster_class}_cluster.png'\n",
    "            ax = polygon.plot(figsize=(10, 10), \\\n",
    "                                        column= visualize_column, \\\n",
    "                                        facecolor='none',edgecolor='cyan' , \\\n",
    "                                        linewidth=1.5 , \\\n",
    "                                        alpha=0.6, categorical=True , cmap=cmap, \\\n",
    "                                        legend=True,  # Add legend \n",
    "                                        legend_kwds={'loc':'upper right', \n",
    "                                                        'bbox_to_anchor':(1, 1), \n",
    "                                                        'fmt':legends_format,\n",
    "                                                        'markerscale':1.01, \n",
    "                                                        'title_fontsize':'small', \n",
    "                                                        'fontsize':'x-small'\n",
    "                                                        # title_fontsizeint or {'xx-small', 'x-small', 'small'\n",
    "                                                        } ,\n",
    "                                        aspect=1\n",
    "                                                )\n",
    "            # vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "            polygon_buffer = polygon.to_crs(crs).buffer(0.016) #for zoom out\n",
    "            minx, miny, maxx, maxy = polygon_buffer.total_bounds\n",
    "            print(\"Pre-non bounds\", minx, miny, maxx, maxy)\n",
    "            if not [x for x in (minx, miny, maxx, maxy) if x is None]:\n",
    "                print(minx, miny, maxx, maxy)\n",
    "                ax.set_xlim(minx, maxx)\n",
    "                ax.set_ylim(miny, maxy)\n",
    "\n",
    "                cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery,  crs=crs) \n",
    "                # vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "                # cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite,  crs=crs , zoom=12) \n",
    "                cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels,   crs=crs, zoom=11) # zoom=13\n",
    "\n",
    "                plt.title(f'{title}',fontsize=16)\n",
    "                plt.xticks([])\n",
    "                plt.yticks([]) \n",
    "                leg1 = ax.get_legend()\n",
    "                # Set markers to square shape\n",
    "                for ea in leg1.legendHandles:\n",
    "                    ea.set_marker('s')\n",
    "                leg1.set_title(f'{legend_title}')\n",
    "                ax.title.set_text(f'{title}')\n",
    "                plt.tight_layout()\n",
    "                # ax.figure.savefig(map_output)\n",
    "                ax.figure.savefig(map_output,  bbox_inches='tight' ,   dpi = 400)\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "\n",
    "def map_func(maps, vector_file,  hexagons, legend_title,legends_format, crs ,cmap,  visualize_column, title, map_output, scheme, categorical):\n",
    "    if categorical==True:\n",
    "        ax = hexagons.to_crs(crs).plot(figsize=(10, 10), \\\n",
    "                                    column= visualize_column, \\\n",
    "                                    alpha=0.6, categorical=True , cmap=cmap, \\\n",
    "                                    legend=True,  # Add legend \n",
    "                                    legend_kwds={'loc':'upper right', \n",
    "                                                    'bbox_to_anchor':(1, 1), \n",
    "                                                    'fmt':legends_format,\n",
    "                                                    'markerscale':1.01, \n",
    "                                                    'title_fontsize':'small', \n",
    "                                                    'fontsize':'x-small'\n",
    "                                                    # title_fontsizeint or {'xx-small', 'x-small', 'small'\n",
    "                                                    } ,\n",
    "                                    aspect=1\n",
    "                                            )\n",
    "    elif categorical==False:\n",
    "        ax = hexagons.to_crs(crs).plot(figsize=(10, 10), \\\n",
    "                                    column= visualize_column, \\\n",
    "                                    alpha=0.6, scheme=scheme, cmap=cmap, \\\n",
    "                                    legend=True,  # Add legend \n",
    "                                    legend_kwds={'loc':'upper right', \n",
    "                                                    'bbox_to_anchor':(1, 1), \n",
    "                                                    'fmt':legends_format,\n",
    "                                                    'markerscale':1.01, \n",
    "                                                    'title_fontsize':'small', \n",
    "                                                    'fontsize':'x-small'\n",
    "                                                    # title_fontsizeint or {'xx-small', 'x-small', 'small'\n",
    "                                                    } ,\n",
    "                                    aspect=1\n",
    "                                            )\n",
    "    # vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "    # cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery,  crs=crs) \n",
    "    vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "    cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite,  crs=crs , zoom=12) \n",
    "    cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels,   crs=crs, zoom=11) # zoom=13\n",
    "\n",
    "    # plt.title(f'{title}',fontsize=16)\n",
    "    plt.xticks([])\n",
    "    plt.yticks([]) \n",
    "    leg1 = ax.get_legend()\n",
    "    # Set markers to square shape\n",
    "    for ea in leg1.legendHandles:\n",
    "        ea.set_marker('s')\n",
    "    leg1.set_title(f'{legend_title}')\n",
    "    ax.title.set_text(f'{title}')\n",
    "    plt.tight_layout()\n",
    "    # ax.figure.savefig(map_output)\n",
    "    ax.figure.savefig(map_output,  bbox_inches='tight' ,   dpi = 400)\n",
    "\n",
    "\n",
    "def map_land_cover(vector_file,  hexagons,title, legend_title, crs,  visualize_column, map_output):\n",
    "    class_labels = {1: \"Water , blue\", 2: \"Trees, green\", 4:\"Flooded Vegetation , orange\", 5:\"Crops, seagreen\", 7:\"Built Area, tan\", 8:\"Bare Ground, goldenrod\", 9:\"Snow/Ice, snow\", 10:\"Clouds, ghostwhite\", 11:\"Rangeland, greenyellow\" }\n",
    "    hexagons[['class', 'color']] = hexagons['esri_dyna'].map(class_labels).str.split(', ',expand=True).values\n",
    "    color_map= dict(zip(hexagons[\"class\"], hexagons[\"color\"]))\n",
    "    ax = hexagons.to_crs(crs).plot(figsize=(10, 10), \\\n",
    "                                column= visualize_column, \\\n",
    "                                alpha=0.6, categorical=True , color=hexagons.color, \\\n",
    "                                legend=True,  # Add legend \n",
    "                                legend_kwds={'loc':'upper right', \n",
    "                                                'bbox_to_anchor':(1, 1), \n",
    "                                                'fmt':legends_format,\n",
    "                                                'markerscale':1.01, \n",
    "                                                'title_fontsize':'small', \n",
    "                                                'fontsize':'x-small'\n",
    "                                                # title_fontsizeint or {'xx-small', 'x-small', 'small'\n",
    "                                                } ,\n",
    "                                aspect=1\n",
    "                                        )\n",
    "    \n",
    "    vector_file.to_crs(crs).plot(ax=ax,facecolor='none',edgecolor='k',alpha=0.2)\n",
    "    cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite,  crs=crs , zoom=12) \n",
    "    cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels,   crs=crs, zoom=11) # zoom=13\n",
    "    plt.title(f'{title}',fontsize=16)\n",
    "    plt.xticks([])\n",
    "    plt.yticks([]) \n",
    "    handles = [Line2D([0], [0], marker='s', color='w', markerfacecolor=k, label=v, markersize=8) for v, k in color_map.items()]\n",
    "    ax.legend(title=legend_title, handles=handles, bbox_to_anchor=(1, 1),markerscale=1.01,title_fontsize='small',fontsize='x-small', loc='upper right')\n",
    "    plt.tight_layout()\n",
    "    # ax.figure.savefig(map_output)\n",
    "    ax.figure.savefig(map_output,  bbox_inches='tight' ,   dpi = 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**OSM Network: Get Travel Time and Shortest path from a given point or an edge of a hotspot**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_network_stats(polygon):\n",
    "     try:\n",
    "          # polygon = cluster_boundaries.iloc[index]['geometry'] \n",
    "          graph = ox.graph.graph_from_polygon(polygon, network_type='drive')\n",
    "          # Retrieve only edges from the graph\n",
    "          edges = ox.graph_to_gdfs(graph, nodes=False, edges=True)\n",
    "          graph_proj = ox.project_graph(graph)\n",
    "          # fig, ax = ox.plot_graph(graph_proj) \n",
    "          # Get Edges and Nodes\n",
    "          nodes_proj, edges_proj = ox.graph_to_gdfs(graph_proj, nodes=True, edges=True)\n",
    "          # print(\"Coordinate system:\", edges_proj.crs)\n",
    "          # edges_proj.head()\n",
    "          # Get Edges and Nodes\n",
    "          nodes_proj, edges_proj = ox.graph_to_gdfs(graph_proj, nodes=True, edges=True)\n",
    "          # print(\"Network Coordinate system: projected\", edges_proj.crs)\n",
    "          # edges_proj.head()\n",
    "          # Get the Convex Hull of the network\n",
    "          convex_hull = edges_proj.unary_union.convex_hull\n",
    "          # Calculate the area\n",
    "          area = convex_hull.area\n",
    "          # Calculate statistics with density information\n",
    "          stats = ox.stats.basic_stats(graph_proj, area=area) #replace area by city polygon\n",
    "          stats =pd.Series(stats) #edge_density_km(road density) - edge_length_total per sq km\n",
    "     except:\n",
    "          stats=0\n",
    "          stats =pd.Series(stats) \n",
    "          pass     \n",
    "     return stats\n",
    "\n",
    "def get_osm_network():\n",
    "    G = ox.graph_from_place(place, network_type='drive')\n",
    "    # impute speed on all edges missing data\n",
    "    G = ox.add_edge_speeds(G)\n",
    "    # calculate travel time (seconds) for all edges\n",
    "    G = ox.add_edge_travel_times(G)\n",
    "    # see mean speed/time values by road type\n",
    "    edges = ox.graph_to_gdfs(G, nodes=False)\n",
    "    edges['highway'] = edges['highway'].astype(str)\n",
    "    edges.groupby('highway')[['length', 'speed_kph', 'travel_time']].mean().round(1)\n",
    "\n",
    "    # same thing again, but this time pass in a few default speed values (km/hour)\n",
    "    # to fill in edges with missing `maxspeed` from OSM\n",
    "    hwy_speeds = {'residential': 35,\n",
    "                'secondary': 50,\n",
    "                'tertiary': 60}\n",
    "    G = ox.add_edge_speeds(G, hwy_speeds)\n",
    "    G = ox.add_edge_travel_times(G)\n",
    "    # calculate two routes by minimizing travel distance vs travel time\n",
    "    orig = list(G)[1]\n",
    "    dest = list(G)[-1]\n",
    "    route1 = nx.shortest_path(G, orig, dest, weight='length')\n",
    "    route2 = nx.shortest_path(G, orig, dest, weight='travel_time')\n",
    "\n",
    "    # compare the two routes\n",
    "    route1_length = int(sum(ox.utils_graph.get_route_edge_attributes(G, route1, 'length')))\n",
    "    route2_length = int(sum(ox.utils_graph.get_route_edge_attributes(G, route2, 'length')))\n",
    "    route1_time = int(sum(ox.utils_graph.get_route_edge_attributes(G, route1, 'travel_time')))\n",
    "    route2_time = int(sum(ox.utils_graph.get_route_edge_attributes(G, route2, 'travel_time')))\n",
    "    print('Route 1 is', route1_length, 'meters and takes', route1_time, 'seconds.')\n",
    "    print('Route 2 is', route2_length, 'meters and takes', route2_time, 'seconds.')\n",
    "\n",
    "    # pick route colors\n",
    "    c1 = 'r' #length\n",
    "    c2 = 'b' #travel_time\n",
    "    rc1 = [c1] * (len(route1) - 1)\n",
    "    rc2 = [c2] * (len(route2) - 1)\n",
    "    rc = rc1 + rc2\n",
    "    nc = [c1, c1, c2, c2]\n",
    "\n",
    "    # plot the routes\n",
    "    fig, ax = ox.plot_graph_routes(G, [route1, route2], route_color=rc, route_linewidth=6,\n",
    "                                node_size=0, bgcolor='k')\n",
    "    location_point = (33.299896, -111.831638)\n",
    "    G2 = ox.graph_from_point(location_point, dist=400, truncate_by_edge=True)\n",
    "    # impute speed on all edges missing data\n",
    "    G2 = ox.add_edge_speeds(G2)\n",
    "    # calculate travel time (seconds) for all edges\n",
    "    G2 = ox.add_edge_travel_times(G2)\n",
    "    origin = (33.301821, -111.829871)\n",
    "    destination = (33.301402, -111.833108)\n",
    "    origin_node = ox.distance.nearest_nodes(G2, origin[1], origin[0])\n",
    "    destination_node = ox.distance.nearest_nodes(G2, destination[1], destination[0])\n",
    "    route = ox.shortest_path(G2, origin_node, destination_node)\n",
    "    # compare the two routes\n",
    "    route1_length = int(sum(ox.utils_graph.get_route_edge_attributes(G2, route, 'length')))\n",
    "    route1_time = int(sum(ox.utils_graph.get_route_edge_attributes(G2, route, 'travel_time')))\n",
    "    print('Route 1 is', route1_length, 'meters and takes', route1_time, 'seconds.')\n",
    "    fig, ax = ox.plot_graph_route(G2, route, route_color=\"c\", node_size=0)\n",
    "\n",
    "def get_services_infrastructure(center_point, dist):\n",
    "\n",
    "     tags={'leisure':['park', 'nature_reserve', 'protected_area', 'garden']}\n",
    "     park= ox.features.features_from_point(center_point=center_point, tags=tags, dist=dist)\n",
    "     \n",
    "     tags={'amenity':['college', 'dancing_school', 'driving_school', 'kindergarten', \\\n",
    "                    'language_school', 'library', 'surf_school', 'toy_library', 'research_institute', \\\n",
    "                         'training', 'music_school', 'school', 'university'],\n",
    "                         'landuse': 'education'\n",
    "                         }\n",
    "     education= ox.features.features_from_point(center_point=center_point, tags=tags, dist=dist)\n",
    "     \n",
    "     tags={'landuse':['commercial', 'industrial', 'construction', 'retail']                   }\n",
    "     industry= ox.features.features_from_point(center_point=center_point, tags=tags, dist=dist)\n",
    "     \n",
    "     # tags={'building': True}\n",
    "     tags={'building': 'government'}\n",
    "     building= ox.features.features_from_point(center_point=center_point, tags=tags, dist=dist)\n",
    "     return park , education, industry, building\n",
    "\n",
    "def get_services_infrastructure_inside_poly(polygon):\n",
    "     try:\n",
    "          tags={'leisure':['park', 'nature_reserve', 'protected_area', 'garden']}\n",
    "          park=ox.features.features_from_polygon(polygon, tags)\n",
    "     except:\n",
    "          park=0\n",
    "          pass\n",
    "     try:\n",
    "          tags={'amenity':['college', 'dancing_school', 'driving_school', 'kindergarten', \\\n",
    "                         'language_school', 'library', 'surf_school', 'toy_library', 'research_institute', \\\n",
    "                              'training', 'music_school', 'school', 'university'],\n",
    "                              'landuse': 'education'\n",
    "                              }\n",
    "          education=ox.features.features_from_polygon(polygon, tags)\n",
    "     except:\n",
    "          education=0\n",
    "          pass\n",
    "     try:\n",
    "          tags={'landuse':['commercial', 'industrial', 'construction', 'retail']                   }\n",
    "          industry=ox.features.features_from_polygon(polygon, tags)\n",
    "     except:\n",
    "          industry=0\n",
    "          pass\n",
    "     try:\n",
    "          tags={'building': 'government'}\n",
    "          building=ox.features.features_from_polygon(polygon, tags)\n",
    "     except:\n",
    "          building=0\n",
    "          pass\n",
    "     # print(\"park , education, industry, building:\", park , education, industry, building)\n",
    "     return park , education, industry, building\n",
    "\n",
    "\n",
    "def multi2single(gpdf):\n",
    "    gpdf_singlepoly = gpdf[gpdf.geometry.type == 'Polygon']\n",
    "    gpdf_multipoly = gpdf[gpdf.geometry.type == 'MultiPolygon']\n",
    "\n",
    "    for i, row in gpdf_multipoly.iterrows():\n",
    "        Series_geometries = pd.Series(row.geometry)\n",
    "        df = pd.concat([gp.GeoDataFrame(row, crs=gpdf_multipoly.crs).T]*len(Series_geometries), ignore_index=True)\n",
    "        df['geometry']  = Series_geometries\n",
    "        gpdf_singlepoly = pd.concat([gpdf_singlepoly, df])\n",
    "\n",
    "    gpdf_singlepoly.reset_index(inplace=True, drop=True)\n",
    "    return gpdf_singlepoly\n",
    "\n",
    "\n",
    "def create_ahc_knn_clusters(db  , raster_val , WGS84_meters):\n",
    "    db =db.to_crs(WGS84_meters)\n",
    "    db['area_sqm'] = db.geometry.area\n",
    "    cluster_variables = [raster_val]\n",
    "    db_scaled = robust_scale(db[cluster_variables])\n",
    "    w = KNN.from_dataframe(db, k=4) #k for four nearest neighbors\n",
    "    # Specify cluster model with spatial constraint\n",
    "    model = AgglomerativeClustering(\n",
    "        linkage=\"ward\", connectivity=w.sparse, n_clusters=5\n",
    "    )\n",
    "    # Fit algorithm to the data\n",
    "    model.fit(db_scaled)\n",
    "    db[\"ward5wknn\"] = model.labels_\n",
    "    db['ward5wknn_sum'] = db[cluster_variables].groupby(db['ward5wknn']).transform('sum')\n",
    "    db['area_sum'] = db['area_sqm'].groupby(db['ward5wknn']).transform('sum')\n",
    "    db['area_sum'] = db['area_sum'].apply(lambda x: x/10000) #Hectares\n",
    "    visualize_var= 'ward5wknn_sum_per_meter'\n",
    "    db[visualize_var] = db['ward5wknn_sum'].div(db['area_sum']).round(2)\n",
    "    # select the columns that you with to use for the dissolve and that will be retained\n",
    "    cluster_boundaries = db[[raster_val , visualize_var,\"ward5wknn\", \"geometry\", \"area_sqm\", \"area_sum\", \"ward5wknn_sum\"]]\n",
    "    # dissolve the state boundary by region \n",
    "    cluster_boundaries = cluster_boundaries.dissolve(by=\"ward5wknn\" , aggfunc={raster_val: \"sum\", visualize_var: \"sum\" , \"area_sqm\":\"sum\" , \"area_sum\": \"sum\", \"ward5wknn_sum\": \"sum\",} , )\n",
    "    cluster_boundaries=cluster_boundaries.explode().reset_index() #convert multipolygons to single part polygons. Multis cause errors\n",
    "    cluster_boundaries[\"area_for_dupes\"]=cluster_boundaries.geometry.area\n",
    "#     cluster_boundaries=cluster_boundaries.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"ward5wknn\"])\n",
    "    cluster_boundaries=multi2single(cluster_boundaries)\n",
    "    cluster_boundaries = cluster_boundaries.to_crs(crs)\n",
    "    cluster_boundaries[\"centroid\"] = cluster_boundaries[\"geometry\"].centroid\n",
    "    cluster_boundaries['lon'] = cluster_boundaries['centroid'].x\n",
    "    cluster_boundaries['lat'] = cluster_boundaries['centroid'].y\n",
    "    cluster_boundaries['cluster_area_deg'] = cluster_boundaries.geometry.area\n",
    "    return db , cluster_boundaries, visualize_var\n",
    "\n",
    "\n",
    "def create_ahc_knn_clusters_vito(db  , raster_val , WGS84_meters):\n",
    "     db =db.to_crs(WGS84_meters)\n",
    "     db['area_sqm'] = db.geometry.area\n",
    "     cluster_variables = [raster_val]\n",
    "     db_scaled = robust_scale(db[cluster_variables])\n",
    "     w = KNN.from_dataframe(db, k=4) #k for four nearest neighbors\n",
    "     # Specify cluster model with spatial constraint\n",
    "     model = AgglomerativeClustering(\n",
    "          linkage=\"ward\", connectivity=w.sparse, n_clusters=5\n",
    "     )\n",
    "     # Fit algorithm to the data\n",
    "     model.fit(db_scaled)\n",
    "     db[\"ward5wknn\"] = model.labels_\n",
    "     db['ward5wknn_sum'] = db[cluster_variables].groupby(db['ward5wknn']).transform('sum')\n",
    "     db['area_sum'] = db['area_sqm'].groupby(db['ward5wknn']).transform('sum')\n",
    "     db['area_sum'] = db['area_sum'].apply(lambda x: x/10000) #Hectares\n",
    "     visualize_var= 'ward5wknn_sum_per_meter'\n",
    "     db[visualize_var] = db['ward5wknn_sum'].div(db['area_sum']).round(2)\n",
    "     # select the columns that you with to use for the dissolve and that will be retained\n",
    "     cluster_boundaries = db[[raster_val , visualize_var,\"ward5wknn\", \"geometry\", \"area_sqm\", \"area_sum\", \"ward5wknn_sum\"]]\n",
    "     # dissolve the state boundary by region \n",
    "     cluster_boundaries = cluster_boundaries.dissolve(by=\"ward5wknn\" , aggfunc={raster_val: \"sum\", visualize_var: \"sum\" , \"area_sqm\":\"sum\" , \"area_sum\": \"sum\", \"ward5wknn_sum\": \"sum\",} , )\n",
    "     cluster_boundaries=cluster_boundaries.explode().reset_index() #convert multipolygons to single part polygons. Multis cause errors\n",
    "     cluster_boundaries[\"area_for_dupes\"]=cluster_boundaries.geometry.area\n",
    "     # cluster_boundaries=cluster_boundaries.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"ward5wknn\"])\n",
    "     cluster_boundaries=multi2single(cluster_boundaries)\n",
    "     cluster_boundaries = cluster_boundaries.to_crs(crs)\n",
    "     cluster_boundaries[\"centroid\"] = cluster_boundaries[\"geometry\"].centroid\n",
    "     cluster_boundaries['lon'] = cluster_boundaries['centroid'].x\n",
    "     cluster_boundaries['lat'] = cluster_boundaries['centroid'].y\n",
    "     cluster_boundaries['cluster_area_deg'] = cluster_boundaries.geometry.area\n",
    "     return db , cluster_boundaries, visualize_var\n",
    "\n",
    "\n",
    "def dissolve_categorical(db  , raster_val , WGS84_meters):\n",
    "     db =db.to_crs(WGS84_meters)\n",
    "     db['area_sqm'] = db.geometry.area\n",
    "     cluster_variables = [raster_val]\n",
    "     # db_scaled = robust_scale(db[cluster_variables])\n",
    "     # w = KNN.from_dataframe(db, k=4) #k for four nearest neighbors\n",
    "     # # Specify cluster model with spatial constraint\n",
    "     # model = AgglomerativeClustering(\n",
    "     #      linkage=\"ward\", connectivity=w.sparse, n_clusters=15\n",
    "     # )\n",
    "     # # Fit algorithm to the data\n",
    "     # model.fit(db_scaled)\n",
    "     # db[\"ward5wknn\"] = model.labels_\n",
    "     db[\"ward5wknn\"] = db[raster_val]\n",
    "     db['ward5wknn_sum'] = db[cluster_variables].groupby(db['ward5wknn']).transform('sum')\n",
    "     db['area_sum'] = db['area_sqm'].groupby(db['ward5wknn']).transform('sum')\n",
    "     db['area_sum'] = db['area_sum'].apply(lambda x: x/10000) #Hectares\n",
    "     visualize_var= 'ward5wknn_sum_per_meter'\n",
    "     db[visualize_var] = db['ward5wknn_sum'].div(db['area_sum']).round(2)\n",
    "     # select the columns that you with to use for the dissolve and that will be retained\n",
    "     cluster_boundaries = db[[raster_val , visualize_var,\"ward5wknn\", \"geometry\", \"area_sqm\", \"area_sum\", \"ward5wknn_sum\"]]\n",
    "     # dissolve the state boundary by region \n",
    "     cluster_boundaries = cluster_boundaries.dissolve(by=\"ward5wknn\" , aggfunc={raster_val: \"sum\", visualize_var: \"sum\" , \"area_sqm\":\"sum\" , \"area_sum\": \"sum\", \"ward5wknn_sum\": \"sum\",} , )\n",
    "     cluster_boundaries=cluster_boundaries.explode().reset_index() #convert multipolygons to single part polygons. Multis cause errors\n",
    "     cluster_boundaries[\"area_for_dupes\"]=cluster_boundaries.geometry.area\n",
    "     # cluster_boundaries=cluster_boundaries.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"ward5wknn\"])\n",
    "     cluster_boundaries=multi2single(cluster_boundaries)\n",
    "     cluster_boundaries = cluster_boundaries.to_crs(crs)\n",
    "     cluster_boundaries[\"centroid\"] = cluster_boundaries[\"geometry\"].centroid\n",
    "     cluster_boundaries['lon'] = cluster_boundaries['centroid'].x\n",
    "     cluster_boundaries['lat'] = cluster_boundaries['centroid'].y\n",
    "     cluster_boundaries['cluster_area_deg'] = cluster_boundaries.geometry.area\n",
    "     return db , cluster_boundaries, visualize_var\n",
    "\n",
    "\n",
    "\n",
    "def describe_cluster_trees_pop(cluster_boundaries, category, maps):\n",
    "     df = pd.DataFrame(cluster_boundaries.drop(columns='geometry'))\n",
    "     df = df.sort_values('cluster_area_deg', ascending = False).groupby('ward5wknn').head(2)\n",
    "     df = df.reset_index() \n",
    "     df['Rank'] = df['ward5wknn_sum_per_meter'].rank(method='dense', ascending=False)\n",
    "     df= df.sort_values(['Rank'], ascending=[True])\n",
    "     cent= cluster_boundaries.dissolve().centroid \n",
    "     # White house 38.8977° N, 77.0365° W. centroid of the city.\n",
    "     lat1 = float(cent.centroid.y)\n",
    "     long1= float(cent.centroid.x)\n",
    "     list_statements= []\n",
    "     stats_dict= defaultdict(list)\n",
    "     for index, row in df.iterrows():\n",
    "          lat2 = row['lat'] # cent.centroid.y  #38.8893\n",
    "          long2 =row['lon'] # cent.centroid.x  # -77.0506\n",
    "          # print(f\"lat1 {lat1} , long1 {long1}, lat2 {lat2}, long2 {long2}\")\n",
    "          points = calcNSEW(lat1, long1, lat2, long2)\n",
    "          rank= int(row['Rank'])\n",
    "          rank= ordinal(rank)\n",
    "          value= row['ward5wknn_sum_per_meter']\n",
    "          statement= f\"A cluster that is located in the {points} of the city, is ranked {rank} and has the value of {value}\"\n",
    "          list_statements.append(statement)\n",
    "          try:\n",
    "               polygon = cluster_boundaries.iloc[index]['geometry']\n",
    "               cluster_area_deg= cluster_boundaries.iloc[index]['cluster_area_deg']\n",
    "               ward5wknn= cluster_boundaries.iloc[index]['ward5wknn']\n",
    "               park , education, industry, building=get_services_infrastructure_inside_poly(polygon)\n",
    "               stats_dict[\"cluster\"].append(index)\n",
    "               stats_dict[\"ward5wknn\"].append(ward5wknn)\n",
    "               stats_dict[\"points\"].append(points)\n",
    "               stats_dict[\"cluster_area_deg\"].append(cluster_area_deg)\n",
    "               stats_dict[\"rank\"].append(rank)\n",
    "               stats_dict[\"value\"].append(value)\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               stats_dict[\"park\"].append(len(park))\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               stats_dict[\"education\"].append(len(education))\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               stats_dict[\"industry\"].append(len(industry))\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               stats_dict[\"building\"].append(len(building))\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               \n",
    "               park[\"ward5wknn\"]=ward5wknn\n",
    "               park[\"index\"]=index\n",
    "               park[\"points\"]=points\n",
    "               park[\"rank\"]=rank\n",
    "               park[\"value\"]=value\n",
    "               park.to_csv(f\"{tables}/{ward5wknn}_{category}_park_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               education[\"ward5wknn\"]=ward5wknn\n",
    "               education[\"index\"]=index\n",
    "               education[\"points\"]=points\n",
    "               education[\"rank\"]=rank\n",
    "               education[\"value\"]=value\n",
    "               education.to_csv(f\"{tables}/{ward5wknn}_{category}_education_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               industry[\"ward5wknn\"]=ward5wknn\n",
    "               industry[\"index\"]=index\n",
    "               industry[\"points\"]=points\n",
    "               industry[\"rank\"]=rank\n",
    "               industry[\"value\"]=value\n",
    "               industry.to_csv(f\"{tables}/{ward5wknn}_{category}_industry_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               building[\"ward5wknn\"]=ward5wknn\n",
    "               building[\"index\"]=index\n",
    "               building[\"points\"]=points\n",
    "               building[\"rank\"]=rank\n",
    "               building[\"value\"]=value\n",
    "               building.to_csv(f\"{tables}/{ward5wknn}_{category}_building_df.csv\")               \n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               stats=get_network_stats(polygon)\n",
    "               stats = stats.to_frame(0).T\n",
    "               stats.to_csv(f\"{tables}/{ward5wknn}_{category}_infra_network_stats.csv\")\n",
    "          except:\n",
    "               print(f\"no network for {ward5wknn}\")\n",
    "               pass\n",
    "     \n",
    "     statement_df=pd.DataFrame(dict([ (k,pd.Series(v)) for k,v in stats_dict.items() ]))\n",
    "     statement_df.to_csv(f\"{tables}/{category}_amenities_stats_all.csv\")\n",
    "     return list_statements\n",
    "\n",
    "\n",
    "def describe_db( db_var, category, tables):\n",
    "     try:\n",
    "          cluster_variables= category\n",
    "          k5desc = db_var.groupby('ward5wknn')[cluster_variables].describe()\n",
    "          df1=k5desc.reset_index(0).reset_index(drop=True)\n",
    "          df1.to_csv(f\"{tables}/{category}_summary_stats_df.csv\")\n",
    "     except:\n",
    "          pass\n",
    "\n",
    "# place = 'Nis, Serbia'\n",
    "# get_osm_network_and_travel_time=get_osm_network()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_landsurface_temperature(city, cities_reprojected, max_value , years, month0 ,  month1):\n",
    "\n",
    "    \"\"\"\n",
    "    download_landsurface_temperature: Exports Land Surface temperature rasters for a given country in the cities shapefile\n",
    "\n",
    "    :param country: Name of the country\n",
    "    :param cities_reprojected: geopandas read cities shapefile\n",
    "    :param max_value: cap the number of cities\n",
    "    :param years: years\n",
    "    :param month0 : first hottest month\n",
    "    :param month1 : last hottest month\n",
    "\n",
    "    :return: confirms the export of rasters\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Create a jsondictionary from the geometry of the  shapefile\n",
    "    geom = cities_reprojected['geometry']\n",
    "    jsonDict = eval(geom.to_json())\n",
    "\n",
    "    for index, row in cities_reprojected.iterrows():\n",
    "        #   print(1)\n",
    "        #   city= row['NAME_1'].lower() #instruct the column for city name in the city shapefile\n",
    "          if index  <= max_value:\n",
    "              city_number= index\n",
    "              print(1)\n",
    "              for x in jsonDict['features'][city_number]['geometry']['coordinates']:\n",
    "                  try:\n",
    "                      AOI = ee.Geometry.Polygon(x) #cast polygon as ee geometry\n",
    "                      landsat = ee.ImageCollection(\"LANDSAT/LC08/C02/T1_L2\")\n",
    "                      #   print(3)\n",
    "                      # Filter for hottest months in the past X years\n",
    "                      def filter_hot_month(i):\n",
    "                          return ee.Filter.date(years[i] + '-' + month0 + '-01', years[i] + '-' + month1 + '-01')\n",
    "\n",
    "                      range_list= list(map(filter_hot_month, list(range(0, 10)))) #combination of months and years\n",
    "                      rangefilter = ee.Filter.Or(range_list)\n",
    "                    #   print(4)\n",
    "                      # Define a function to scale the data and mask unwanted pixels\n",
    "                      def maskL457sr(image):\n",
    "                          # Bit 0 - Fill\n",
    "                          # Bit 1 - Dilated Cloud\n",
    "                          # Bit 2 - Cirrus (high confidence)\n",
    "                          # Bit 3 - Cloud\n",
    "                          # Bit 4 - Cloud Shadow\n",
    "                          qaMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0)\n",
    "                          saturationMask = image.select('QA_RADSAT').eq(0)\n",
    "                          # Apply the scaling factors to the appropriate bands.\n",
    "                          thermalBand = image.select('ST_B10').multiply(0.00341802).add(149.0)\n",
    "                          # Replace the original bands with the scaled ones and apply the masks.\n",
    "                          return image.addBands(thermalBand, None, True).updateMask(qaMask).updateMask(saturationMask)\n",
    "                      #   print(5)\n",
    "                      # Apply filter and mask\n",
    "                      collectionSummer = landsat.filter(rangefilter).filterBounds(AOI).map(maskL457sr).select('ST_B10').mean().add(-273.15).clip(AOI)\n",
    "                      # collectionSummer=collectionSummer.toFloat()\n",
    "\n",
    "                      lst_tif = os.path.join(rasters, 'lst_projected.tif')\n",
    "                      geemap.ee_export_image( \n",
    "                                            # collectionSummer.toUint8(), \n",
    "                                            collectionSummer.toFloat(), \n",
    "                                            filename=lst_tif,  \n",
    "                                            crs=EPSG_str,\n",
    "                                            # crs_transform=crs_transform,\n",
    "                                            scale=30, \n",
    "                                            region=AOI\n",
    "                                            # file_per_band=False\n",
    "                                            )\n",
    "                        \n",
    "                    #   print(6)\n",
    "                  except Exception as e:\n",
    "                      print(f\"Error : {e} \")\n",
    "                  return  lst_tif\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "if __name__ ==  '__main__': \n",
    "    if run_tree_and_pop_analysis:\n",
    "        # # City level map\n",
    "        df_filtered= gdf_city\n",
    "        title=\"City Boundary\"\n",
    "        legend_title= 'Boundary'\n",
    "        visualize_column=\"shapeName\"\n",
    "        map_output= os.path.join(maps, f\"{city}_{title}_map.png\")\n",
    "        export = map_city(maps=maps, vector_file=gdf_city, crs=crs, legend_title=legend_title,  \\\n",
    "                            visualize_column=visualize_column, title=title, map_output=map_output)\n",
    "\n",
    "        crs = 4326\n",
    "        fc_filtered  , region= create_fc_collection(gdf_city= gdf_city) \n",
    "\n",
    "        \n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        # #Args\n",
    "        cities_reprojected=gdf_city.reset_index()\n",
    "        max_value=len(cities_reprojected)\n",
    "        years =list(str(i) for i in range(2013, 2024)) #years from 2013-2023\n",
    "        month0 = '06'  # first hottest month  # update for each country\n",
    "        month1 = '09'  # end of hottest month (note that this is exclusive)  # update for each country\n",
    "        country= 'SRB' # replace with the relevent country name\n",
    "        lst_t= download_landsurface_temperature( city, cities_reprojected, max_value, years, month0, month1)\n",
    "        raster_filename= \"land_surface_temperature\"\n",
    "        unprojected_raster =lst_t  \n",
    "        projected_raster = os.path.join(rasters, f'{raster_filename}.tif')\n",
    "        reproject_rast= reproject_rasters(crs=crs, output_dir=rasters, unprojected_raster=unprojected_raster, \\\n",
    "                                projected_raster=projected_raster)\n",
    "\n",
    "        output_shapefile = f'{shapefiles}/polygonized_{raster_filename}.shp'\n",
    "        mask_value=None\n",
    "        variable= raster_filename\n",
    "        tif = projected_raster #os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "        raster_val= variable[0:9]\n",
    "        name = polygonize(raster_file=projected_raster ,output_shapefile=output_shapefile ,crs=EPSG_str,  \\\n",
    "                            mask_value=mask_value, raster_val=raster_val)  \n",
    "\n",
    "        cmap = \"OrRd\"\n",
    "        # cmap = \"Greens\"\n",
    "        # cmap = \"YlGn\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        legend_title= \"LST\"\n",
    "        map_title= f\"Land Surface Temperature across {city}\"\n",
    "        scheme=\"quantiles\"\n",
    "        categorical=False\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        polygonized_shapefile = gpd.read_file(output_shapefile).to_crs(crs)\n",
    "        newdf = overlay(polygonized_shapefile, gdf_city, how=\"intersection\")\n",
    "\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                hexagons=newdf,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                    crs=crs, cmap=cmap , visualize_column=raster_val,  \\\n",
    "                        title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        db , cluster_boundaries, visualize_var =create_ahc_knn_clusters(db = newdf , raster_val=raster_val, \\\n",
    "                                                                        WGS84_meters=WGS84_meters)\n",
    "        map_title= f\"Land Surface Temperature across {city} clusters\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        categorical=True\n",
    "        map_output= f'{maps}/{city}_map_clusters_{raster_filename}.png'\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "            hexagons=db,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                crs=crs, cmap=cmap , visualize_column=visualize_var,  \\\n",
    "                    title=map_title, map_output=map_output,  scheme=scheme,  categorical=categorical)\n",
    "        category=raster_filename\n",
    "        describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "        describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "        delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "\n",
    "\n",
    "        # legend_title= \"Surface Temperature\"\n",
    "        # map_title= f\"Land Surface Temperature across {city} cluster\"\n",
    "        # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "        #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "        #                     visualize_column=raster_val, title=map_title)\n",
    "                        \n",
    "        cluster_boundaries_lst= cluster_boundaries\n",
    "        db_lst = db[[f\"{raster_val}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        tree_cover_tif_esa= download_tree_cover_esa_vito(ee , EPSG_str, region)\n",
    "        # raster_filename= \"tree_cover_esa\"\n",
    "        raster_filename= \"esa_tree_cover\"\n",
    "        unprojected_raster =tree_cover_tif_esa  \n",
    "        projected_raster = os.path.join(rasters, f'{raster_filename}.tif')\n",
    "        reproject_rast= reproject_rasters(crs=crs, output_dir=rasters, unprojected_raster=unprojected_raster, \\\n",
    "                                            projected_raster=projected_raster)\n",
    "\n",
    "        output_shapefile = f'{shapefiles}/polygonized_{raster_filename}.shp'\n",
    "        mask_value=None\n",
    "        variable= raster_filename\n",
    "        tif = projected_raster #os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "        raster_val= variable[0:9]\n",
    "        name = polygonize(raster_file=projected_raster ,output_shapefile=output_shapefile ,crs=EPSG_str,  \\\n",
    "            mask_value=mask_value, raster_val=raster_val)  \n",
    "\n",
    "        # cmap = \"OrRd\"\n",
    "        cmap = \"Greens\"\n",
    "        # cmap = \"YlGn\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        legend_title= \"Tree Cover\"\n",
    "        map_title= f\"Tree Cover (ESA 10 meters) across {city}\"\n",
    "        scheme=\"quantiles\"\n",
    "        categorical=False\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        polygonized_shapefile = gpd.read_file(output_shapefile).to_crs(crs)\n",
    "        newdf = overlay(polygonized_shapefile, gdf_city, how=\"intersection\")\n",
    "\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                hexagons=newdf,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                    crs=crs, cmap=cmap , visualize_column=raster_val,  \\\n",
    "                        title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        db , cluster_boundaries, visualize_var =create_ahc_knn_clusters(db = newdf, raster_val=raster_val, \\\n",
    "                                                                        WGS84_meters=WGS84_meters)\n",
    "        map_title= f\"Tree Cover (ESA 10 meters) across {city} clusters\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        categorical=True\n",
    "        map_output= f'{maps}/{city}_map_clusters_{raster_filename}.png'\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "            hexagons=db,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                crs=crs, cmap=cmap , visualize_column=visualize_var,  \\\n",
    "                    title=map_title, map_output=map_output,  scheme=scheme,  categorical=categorical)\n",
    "        category=raster_filename\n",
    "        describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "        describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "        delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "\n",
    "\n",
    "        # legend_title= \"Tree Cover\"\n",
    "        # map_title= f\"Tree Cover across {city} cluster\"\n",
    "        # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "        #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "        #                     visualize_column=raster_val, title=map_title)\n",
    "\n",
    "\n",
    "\n",
    "        db_esa = db[[f\"{raster_val}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        tree_cover_tif= download_tree_cover(ee , EPSG_str, region)\n",
    "        raster_filename= \"land_tree\"\n",
    "        unprojected_raster =tree_cover_tif  \n",
    "        projected_raster = os.path.join(rasters, f'{raster_filename}.tif')\n",
    "        reproject_rast= reproject_rasters(crs=crs, output_dir=rasters, unprojected_raster=unprojected_raster, \\\n",
    "                                          projected_raster=projected_raster)\n",
    "\n",
    "        output_shapefile = f'{shapefiles}/polygonized_{raster_filename}.shp'\n",
    "        mask_value=None\n",
    "        variable= raster_filename\n",
    "        tif = projected_raster#os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "        raster_val= variable[0:9]\n",
    "        name = polygonize(raster_file=projected_raster ,output_shapefile=output_shapefile ,crs=EPSG_str,  \\\n",
    "            mask_value=mask_value, raster_val=raster_val)  \n",
    "\n",
    "        # cmap = \"OrRd\"\n",
    "        cmap = \"Greens\"\n",
    "        # cmap = \"YlGn\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        legend_title= \"Tree Cover\"\n",
    "        map_title= f\"Tree Cover (Landsat 30 meters) across {city}\"\n",
    "        scheme=\"quantiles\"\n",
    "        categorical=False\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        polygonized_shapefile = gpd.read_file(output_shapefile).to_crs(crs)\n",
    "        newdf = overlay(polygonized_shapefile, gdf_city, how=\"intersection\")\n",
    "\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                hexagons=newdf,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                    crs=crs, cmap=cmap , visualize_column=raster_val,  \\\n",
    "                        title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        db , cluster_boundaries, visualize_var =create_ahc_knn_clusters(db = newdf , raster_val=raster_val, \\\n",
    "                                                                        WGS84_meters=WGS84_meters)\n",
    "        map_title= f\"Tree Cover (Landsat 30 meters) across {city} clusters\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        categorical=True\n",
    "        map_output= f'{maps}/{city}_map_clusters_{raster_filename}.png'\n",
    "\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "            hexagons=db,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                crs=crs, cmap=cmap , visualize_column=visualize_var,  \\\n",
    "                    title=map_title, map_output=map_output,  scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        category=raster_filename\n",
    "        describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "        describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "        delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "\n",
    "\n",
    "        # legend_title= \"Tree Cover (Landsat 30 meters)\"\n",
    "        # map_title= f\"Tree Cover (Landsat 30 meters) across {city} cluster\"\n",
    "        # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "        #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "        #                     visualize_column=raster_val, title=map_title)\n",
    "\n",
    "        db_landsat = db[[f\"{raster_val}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        pop_density_tif= dowload_pop_density(ee , EPSG_str, region)\n",
    "        raster_filename= \"pop_density\"\n",
    "        unprojected_raster =pop_density_tif  \n",
    "        projected_raster = os.path.join(rasters, f'{raster_filename}.tif')\n",
    "        reproject_rast= reproject_rasters(crs=crs, output_dir=rasters, unprojected_raster=unprojected_raster, \\\n",
    "            projected_raster=projected_raster)\n",
    "\n",
    "        output_shapefile = f'{shapefiles}/polygonized_{raster_filename}.shp'\n",
    "        mask_value=None\n",
    "        variable= raster_filename\n",
    "        tif = projected_raster #os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "        raster_val= variable[0:9]\n",
    "        name = polygonize(raster_file=projected_raster ,output_shapefile=output_shapefile ,crs=EPSG_str,  \\\n",
    "            mask_value=mask_value, raster_val=raster_val)  \n",
    "\n",
    "        cmap = \"OrRd\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        legend_title= \"Population\"\n",
    "        map_title= f\"Population (Worldpop 100 meters) across {city}\"\n",
    "        scheme=\"quantiles\"\n",
    "        categorical=False\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        polygonized_shapefile = gpd.read_file(output_shapefile).to_crs(crs)\n",
    "        newdf = overlay(polygonized_shapefile, gdf_city, how=\"intersection\")\n",
    "\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                hexagons=newdf,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                    crs=crs, cmap=cmap , visualize_column=raster_val,  \\\n",
    "                        title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        db , cluster_boundaries, visualize_var =create_ahc_knn_clusters(db = newdf , raster_val=raster_val, \\\n",
    "                                                                        WGS84_meters=WGS84_meters)\n",
    "\n",
    "        map_title= f\"Population (Worldpop 100 meters) \\n across {city} clusters\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        categorical=True\n",
    "        map_output= f'{maps}/{city}_map_clusters_{raster_filename}.png'\n",
    "        export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "            hexagons=db,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                crs=crs, cmap=cmap , visualize_column=visualize_var,  \\\n",
    "                    title=map_title, map_output=map_output,  scheme=scheme,  categorical=categorical)\n",
    "\n",
    "        category=raster_filename\n",
    "        describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "        describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "        delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "\n",
    "        # legend_title= \"Population\"\n",
    "        # map_title= f\"Population (Worldpop 100 meters) across {city} cluster\"\n",
    "        # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "        #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "        #                     visualize_column=raster_val, title=map_title)\n",
    "\n",
    "\n",
    "        db_pop = db[[f\"{raster_val}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        # google_dynamic_worldcover= download_google_dynamic_worldcover(ee , EPSG_str, region)\n",
    "        esri_dynamic_worldcover= download_esri_dynamic_worldcover(ee , EPSG_str, region)\n",
    "\n",
    "        # Esri cover\n",
    "        raster_filename= \"esri_dynamic_worldcover\"\n",
    "        unprojected_raster =esri_dynamic_worldcover  \n",
    "        projected_raster = os.path.join(rasters, f'{raster_filename}.tif')\n",
    "        reproject_rast= reproject_rasters(crs=crs, output_dir=rasters, unprojected_raster=unprojected_raster, \\\n",
    "            projected_raster=projected_raster)\n",
    "\n",
    "        output_shapefile = f'{shapefiles}/polygonized_{raster_filename}.shp'\n",
    "        mask_value=None\n",
    "        variable= raster_filename\n",
    "        tif = projected_raster #os.path.join(rasters, 'tree_cover_projected.tif')\n",
    "        raster_val= variable[0:9]\n",
    "        name = polygonize(raster_file=projected_raster ,output_shapefile=output_shapefile ,crs=EPSG_str,  \\\n",
    "                            mask_value=mask_value, raster_val=raster_val)  \n",
    "\n",
    "        cmap = \"OrRd\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        legend_title= \"Land Cover\"\n",
    "        map_title= f\"Land Cover (ESRI 10 meters) across {city}\"\n",
    "        scheme=\"quantiles\"\n",
    "        categorical=False\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        polygonized_shapefile = gpd.read_file(output_shapefile).to_crs(crs)\n",
    "        newdf = overlay(polygonized_shapefile, gdf_city, how=\"intersection\")\n",
    "        map_output= f'{maps}/{city}_map_{raster_filename}.png'\n",
    "        \n",
    "        export_the_map= map_land_cover(vector_file=gdf_city , \\\n",
    "                                        hexagons=newdf,title=map_title, legend_title=legend_title, \\\n",
    "                                            crs=crs,  visualize_column=raster_val, map_output=map_output)\n",
    "\n",
    "\n",
    "        db , cluster_boundaries, visualize_var =dissolve_categorical(db = newdf , raster_val=raster_val, \\\n",
    "                                                                        WGS84_meters=WGS84_meters)\n",
    "        map_title= f\"Land Cover (ESRI 10 meters) across {city} clusters\"\n",
    "        legends_format= '{:,.0f}'\n",
    "        categorical=True\n",
    "        map_output= f'{maps}/{city}_map_clusters_{raster_filename}.png'\n",
    "        export_the_map= map_land_cover(vector_file=gdf_city , \\\n",
    "                                    hexagons=db,title=map_title, legend_title=legend_title, \\\n",
    "                                        crs=crs,  visualize_column=raster_val, map_output=map_output)\n",
    "\n",
    "        category=raster_filename\n",
    "        describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "        describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "        delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "\n",
    "        # legend_title= \"Land Cover\"\n",
    "        # map_title= f\"Land Cover (ESRI 10 meters) across {city} cluster\"\n",
    "        # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "        #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "        #                     visualize_column=raster_val, title=map_title)\n",
    "\n",
    "\n",
    "        db_lc = db[[f\"{raster_val}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "        # Put everything togther\n",
    "        #-----------------------------------------------------------------------------------#\n",
    "\n",
    "        # cluster_boundaries_lst=cluster_boundaries_lst.to_crs(crs)\n",
    "        # db_lc=db_lc.to_crs(crs)\n",
    "        # db_lc['geometry'] = db_lc['geometry'].centroid\n",
    "        # df = gpd.sjoin(cluster_boundaries_lst, db_lc, how='left', op='contains')\n",
    "        # db_pop= db_pop.to_crs(crs)\n",
    "        # df = gpd.sjoin(df, db_pop, how='left', op='contains')\n",
    "        # db_landsat= db_landsat.to_crs(crs)\n",
    "        # df = gpd.sjoin(df, db_landsat, how='left', op='contains')\n",
    "        # db_esa= db_esa.to_crs(crs)\n",
    "        # df = gpd.sjoin(df, db_esa, how='left', op='contains')\n",
    "\n",
    "        # # # # Place all other dbs i  heat and then do some regs and descriptive stats\n",
    "        # # cluster_boundaries_lst=cluster_boundaries_lst.to_crs(crs)\n",
    "        # # raster_filename= \"land_surface_temperature\"\n",
    "        # # db_var=db_var.to_crs(crs)\n",
    "        # # db_var['geometry'] = db_var['geometry'].centroid\n",
    "        # # df_all= overlay(cluster_boundaries_lst, db_var, how=\"intersection\" , keep_geom_type=False)\n",
    "        # # df_all=df_all.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"geometry\"])\n",
    "        # # df_all\n",
    "        # #-----------------------------------------------------------#\n",
    "        # #-----------------------------------------------------------#\n",
    "        # # Group table by cluster label, keep the variables used \n",
    "        # # for clustering, and obtain their descriptive summary\n",
    "        # cluster_variables=[]\n",
    "        # k5desc = db.groupby('ward5wknn')[cluster_variables].describe()\n",
    "        # # Loop over each cluster and print a table with descriptives\n",
    "        # for cluster in k5desc.T:\n",
    "        #     print('\\n\\t---------\\n\\tCluster %i'%cluster)\n",
    "        #     print(k5desc.T[cluster])\n",
    "            \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section3'></a>\n",
    "<h5><center> <font color='cyan'> Section 3: Visualize Vito heat rasters at city level</font>   </center></h5>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def create_hexbins_and_geom_features(gdf_city, hexbin_res, WGS84):\n",
    "    poly = gdf_city.geometry.unary_union\n",
    "    gdf_boundary = gpd.GeoDataFrame(geometry=[poly],crs=gdf_city.crs)\n",
    "    # gdf_h3 = gdf_boundary.h3.polyfill(9,explode=True)\n",
    "    gdf_h3 = gdf_boundary.h3.polyfill(hexbin_res,explode=True)\n",
    "    gdf_h3 = gdf_h3.set_index('h3_polyfill').h3.h3_to_geo_boundary()\n",
    "    output_shapefile= f'{shapefiles}/h3_grid.shp'\n",
    "    gdf_h3.to_file(output_shapefile)\n",
    "    hexbins_projected = gpd.read_file(output_shapefile).to_crs(WGS84)\n",
    "\n",
    "    geom = gdf_boundary['geometry']\n",
    "    jsonDict = eval(geom.to_json())\n",
    "    for index, row in gdf_boundary.iterrows(): \n",
    "        polygon_list= []\n",
    "        for x in jsonDict['features'][index]['geometry']['coordinates']:\n",
    "            polygon_list.append(x)\n",
    "            region = ee.Geometry.Polygon(polygon_list)\n",
    "            fc_filtered = ee.FeatureCollection(region)  \n",
    "    # gdf_h3.plot()\n",
    "    return fc_filtered  , region , hexbins_projected\n",
    "\n",
    "def clip_and_export_raster(input_gpdf, raster , out_raster, zonal_stat_var, hexbins_projected):\n",
    "    Vector=input_gpdf\n",
    "    with rasterio.open(raster) as src:\n",
    "        Vector=Vector.to_crs(src.crs)\n",
    "        # print(f'Gdf srs: {Vector.crs}')\n",
    "        out_image, out_transform=mask(src,Vector.geometry,crop=True)\n",
    "        out_meta=src.meta.copy() # copy the metadata of the source DEM\n",
    "        \n",
    "    out_meta.update({\n",
    "        \"driver\":\"Gtiff\",\n",
    "        \"height\":out_image.shape[1], # height starts with shape[1]\n",
    "        \"width\":out_image.shape[2], # width starts with shape[2]\n",
    "        \"transform\":out_transform\n",
    "    })\n",
    "                \n",
    "    with rasterio.open(out_raster,'w',**out_meta) as dst:\n",
    "        dst.write(out_image)\n",
    "\n",
    "    # add zonal stats to the hexbin gpd dataframe and replace nans with 0 \n",
    "    hexbins_projected[zonal_stat_var] = zonal_stats(hexbins_projected, out_raster ,stats='mean')\n",
    "    hexbins_projected[zonal_stat_var] = [item['mean'] for item in hexbins_projected[zonal_stat_var]]\n",
    "    hexbins_projected[zonal_stat_var].fillna(0,inplace=True)\n",
    "    max=int(hexbins_projected.loc[hexbins_projected[zonal_stat_var].idxmax()][zonal_stat_var])\n",
    "    \n",
    "    return hexbins_projected , max\n",
    "\n",
    "def get_uh_raster_dirs(input_raster_dir):\n",
    "    raster_dirs=  [x[1] for x in os.walk(input_raster_dir)] # ['Present', 'SSP119' , 'SSP370']\n",
    "    raster_dirs = [x for x in raster_dirs if len(x)>0]\n",
    "    raster_dirs = list(np.concatenate(raster_dirs))\n",
    "    return raster_dirs\n",
    "\n",
    "\n",
    "def process_vito( input_raster_dir, gdf_city , crs):\n",
    "    fc_filtered  , region , hexbins_projected= create_hexbins_and_geom_features(gdf_city, hexbin_res=8, WGS84=WGS84)\n",
    "    extension = '.tif'\n",
    "    for root, dirs_list, files_list in os.walk( input_raster_dir):\n",
    "        print(root)\n",
    "        for file_name in files_list:\n",
    "            print(file_name)\n",
    "            if os.path.splitext(file_name)[-1] == extension:\n",
    "                file_name_path = os.path.join(root, file_name)\n",
    "                file_name = os.path.splitext(file_name)[0]\n",
    "                # file_name = ' '.join(file_name.split(\"_\")[:-2])\n",
    "                file_name = ' '.join(file_name.split(\"_\")[:-2])\n",
    "                out_rst = os.path.join(rasters, f'{file_name}.tif')\n",
    "                hexbins_projected , max= clip_and_export_raster(input_gpdf=gdf_city, raster=file_name_path , out_raster=out_rst, \\\n",
    "                                                                    zonal_stat_var=file_name , hexbins_projected=hexbins_projected)\n",
    "                cmap = \"coolwarm\"\n",
    "                # cmap = \"YlGn\"\n",
    "                legends_format= '{:,.02f}'\n",
    "                legend_title= f\"{file_name}\"\n",
    "                map_title= f\"{file_name} \\n across {city}\"\n",
    "                scheme=\"quantiles\"\n",
    "                categorical=False\n",
    "                map_output= f'{maps}/map_{file_name}.png'\n",
    "                \n",
    "                export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                        hexagons=hexbins_projected,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                            crs=crs, cmap=cmap , visualize_column=file_name,  \\\n",
    "                                title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "\n",
    "                # export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                #         hexagons=hexbins_projected,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                #             crs=crs, cmap=cmap , visualize_column=file_name,  \\\n",
    "                #                 title=map_title, map_output=map_output, scheme=scheme,  categorical=categorical)\n",
    "\n",
    "                db , cluster_boundaries, visualize_var =create_ahc_knn_clusters_vito(db = hexbins_projected , \\\n",
    "                                                                                raster_val=file_name, WGS84_meters=WGS84_meters)\n",
    "                map_title= f\"{file_name} \\n across {city} clusters\"\n",
    "                legends_format= '{:,.02f}'\n",
    "                categorical=True\n",
    "                map_output= f'{maps}/map_clusters_{file_name}.png'\n",
    "                export_the_map= map_func(maps= maps, vector_file=gdf_city , \\\n",
    "                    hexagons=db,legend_title=legend_title, legends_format=legends_format,  \\\n",
    "                        crs=crs, cmap=cmap , visualize_column=visualize_var,  \\\n",
    "                            title=map_title, map_output=map_output,  scheme=scheme,  categorical=categorical)\n",
    "                \n",
    "                category=file_name\n",
    "                describe_cluster_trees= describe_cluster_trees_pop(cluster_boundaries, category, maps)\n",
    "                describe=describe_db( db_var=db, category=category, tables=tables)\n",
    "                # delete_shapefil= delete_shapefile(output_shapefile=output_shapefile)\n",
    "                # legend_title= \"Population\"\n",
    "                # map_title= f\"Population (Worldpop 100 meters) across {city} cluster\"\n",
    "                # map=map_clusters(maps=maps, hexagons=cluster_boundaries, legend_title=legend_title, \\\n",
    "                #                 legends_format=legends_format, crs=crs ,cmap=cmap,  \\\n",
    "                #                     visualize_column=raster_val, title=map_title)\n",
    "\n",
    "                db_var= f\"db_{file_name}\"\n",
    "                db_var = db[[f\"{file_name}\",'geometry' , 'ward5wknn']]\n",
    "\n",
    "                return db_var , cluster_boundaries, visualize_var\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "if __name__ ==  '__main__': \n",
    "    if run_vito_analysis:\n",
    "        # input_raster_dir = f'{base_dir}/data/Phnom Penh/Geotiffs'\n",
    "        input_raster_dir = data\n",
    "        # raster_dirs=get_uh_raster_dirs(input_raster_dir)\n",
    "        db_var , cluster_boundaries, visualize_var=process_vito( input_raster_dir=data, gdf_city=gdf_city , crs=crs)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section4'></a>\n",
    "<h5><center> <font color='cyan'> Section 4: Creat the Heat Scan PPT</font>   </center></h5>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_city_map_slide(img_path ,tables, prs, width, height, var):\n",
    "    prs.slide_height = Inches(height)\n",
    "    prs.slide_width  = Inches(width+4)\n",
    "    blank_slide_layout = prs.slide_layouts[6]\n",
    "    slide = prs.slides.add_slide(blank_slide_layout)\n",
    "    left = Inches(4)\n",
    "    top =Inches(0)\n",
    "    width = Inches(width)\n",
    "    height = Inches(height)\n",
    "    pic = slide.shapes.add_picture(img_path, left, top, width = width, height = height)\n",
    "\n",
    "    # Add key message box.\n",
    "    # df_= points\n",
    "\n",
    "    text_width = Inches(4)\n",
    "    text_height = Inches(4)\n",
    "    left =Inches(0)\n",
    "    top =Inches(1)\n",
    "    txBox = slide.shapes.add_textbox(left, top, text_width, text_height)\n",
    "    tf = txBox.text_frame\n",
    "    tf.word_wrap = True\n",
    "    p = tf.add_paragraph()\n",
    "    run = p.add_run()\n",
    "    run.text = f'Key Message: {var}'\n",
    "    font = run.font\n",
    "    font.name = 'Yu Gothic'\n",
    "    font.size = Pt(12)\n",
    "    font.bold = True\n",
    "    run.font.color.rgb = RGBColor(0, 0, 0)\n",
    "    prs.save(f'{tables}/Heat_Scan.pptx')\n",
    "\n",
    "def add_cluster_map_slide(img_path ,tables, prs, width, height, var):\n",
    "    prs.slide_height = Inches(height)\n",
    "    prs.slide_width  = Inches(width+4)\n",
    "    blank_slide_layout = prs.slide_layouts[6]\n",
    "    slide = prs.slides.add_slide(blank_slide_layout)\n",
    "    left = Inches(4)\n",
    "    top =Inches(0)\n",
    "    width = Inches(width)\n",
    "    height = Inches(height)\n",
    "    pic = slide.shapes.add_picture(img_path, left, top, width = width, height = height)\n",
    "\n",
    "    # Add key message box.\n",
    "    text_width = Inches(4)\n",
    "    text_height = Inches(2)\n",
    "    left =Inches(0)\n",
    "    top =Inches(1)\n",
    "    txBox = slide.shapes.add_textbox(left, top, text_width, text_height)\n",
    "    tf = txBox.text_frame\n",
    "    tf.word_wrap = True\n",
    "    p = tf.add_paragraph()\n",
    "    run = p.add_run()\n",
    "    run.text = f'Key Message: For {var}, there are five clusters in the city.'\n",
    "    font = run.font\n",
    "    font.name = 'Yu Gothic'\n",
    "    font.size = Pt(13)\n",
    "    font.bold = True\n",
    "    run.font.color.rgb = RGBColor(0, 0, 0)\n",
    "    summary_df=pd.read_csv(f\"{tables}/{var}_summary_stats_df.csv\")\n",
    "    # summary_df=pd.read_csv(f\"{tables}/CSDI ssp370 2021 2040_summary_stats_df.csv\")\n",
    "\n",
    "    # for cluster in summary_df[\"ward5wknn\"]:\n",
    "    # for cluster in range(0,3):\n",
    "    for index, poi in summary_df.iterrows():\n",
    "        cluster=poi[\"ward5wknn\"]\n",
    "        try:\n",
    "            # print(\"read cluster\")\n",
    "            # var1= \"CSDI ssp370 2021 2040\"\n",
    "            clusters_level=pd.read_csv(f\"{tables}/{var}_amenities_stats_all.csv\")\n",
    "            clusters_level= clusters_level.sort_values(['rank'], ascending=[True]) \n",
    "            orientation= clusters_level[clusters_level['ward5wknn']==cluster]['points'].item()\n",
    "            rank= clusters_level[clusters_level['ward5wknn']==cluster]['rank'].item()\n",
    "        except:\n",
    "            orientation=\"Not Applicable\"\n",
    "            rank=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            parks= pd.read_csv(f\"{tables}/{cluster}_{var}_park_df.csv\") \n",
    "            parks= len(parks)\n",
    "        except:\n",
    "            parks=0\n",
    "            pass\n",
    "        try:\n",
    "            education= pd.read_csv(f\"{tables}/{cluster}_{var}_education_df.csv\") \n",
    "            education= len(education)\n",
    "        except:\n",
    "            education=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            industry= pd.read_csv(f\"{tables}/{cluster}_{var}_industry_df.csv\") \n",
    "            industry= len(industry)\n",
    "        except:\n",
    "            industry=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            buildings= pd.read_csv(f\"{tables}/{cluster}_{var}_building_df.csv\") \n",
    "            buildings= len(buildings)\n",
    "        except:\n",
    "            buildings=0\n",
    "            pass\n",
    "        try:\n",
    "            road_netwrok= pd.read_csv(f\"{tables}/{cluster}_{var}_infra_network_stats.csv\") \n",
    "            road_netwrok= len(road_netwrok)\n",
    "        except:\n",
    "            road_netwrok=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            mean= poi[\"mean\"].round(2)\n",
    "        except:\n",
    "            mean=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            std= poi[\"std\"].round(2)\n",
    "        except:\n",
    "            std=0\n",
    "            pass\n",
    "        try:\n",
    "            min= poi[\"min\"].round(2)\n",
    "        except:\n",
    "            min=0\n",
    "            pass\n",
    "\n",
    "        try:\n",
    "            max= poi[\"max\"].round(2)\n",
    "        except:\n",
    "            max=0\n",
    "            pass\n",
    "\n",
    "\n",
    "        # Add key message box.\n",
    "        text_width = Inches(4)\n",
    "        text_height = Inches(4)\n",
    "        left =Inches(0)\n",
    "        top =Inches(3+cluster)\n",
    "        txBox = slide.shapes.add_textbox(left, top, text_width, text_height)\n",
    "        tf = txBox.text_frame\n",
    "        tf.word_wrap = True\n",
    "        p = tf.add_paragraph()\n",
    "        run = p.add_run()\n",
    "        run.text = f'Cluster {int(cluster)} located in {orientation} of the city is ranked {rank}. It has a mean of {mean}, Standard Dev of {std}, min {min}, max {max}. It has {parks} parks, {education} education, {industry} industry, {buildings} buildings, {road_netwrok} road_network.'\n",
    "        font = run.font\n",
    "        font.name = 'Yu Gothic'\n",
    "        font.size = Pt(10)\n",
    "        font.bold = True\n",
    "        run.font.color.rgb = RGBColor(0, 0, 0)\n",
    "    prs.save(f'{tables}/Heat_Scan.pptx')\n",
    "\n",
    "\n",
    "def add_title_plus_big_picture_slide(city,tables, prs, width, height):\n",
    "    prs.slide_height = Inches(height)\n",
    "    prs.slide_width  = Inches(width)\n",
    "    blank_slide_layout = prs.slide_layouts[1]\n",
    "    slide = prs.slides.add_slide(blank_slide_layout)\n",
    "\n",
    "    shapes = slide.shapes\n",
    "    title_shape = shapes.title\n",
    "    body_shape = shapes.placeholders[1]\n",
    "    title_shape.text = f'{city} Heat Scan'\n",
    "    title_shape.text_frame.paragraphs[0].alignment = PP_ALIGN.CENTER\n",
    "    # title_shape.left = Inches(4)\n",
    "\n",
    "    tf = body_shape.text_frame\n",
    "    tf.text = 'Find the bullet slide layout'\n",
    "\n",
    "    p = tf.add_paragraph()\n",
    "    p.text = 'Use _TextFrame.text for first bullet'\n",
    "    p.level = 1\n",
    "\n",
    "    p = tf.add_paragraph()\n",
    "    p.text = 'Use _TextFrame.add_paragraph() for subsequent bullets'\n",
    "    p.level = 2\n",
    "    prs.save(f'{tables}/Heat_Scan.pptx')\n",
    "\n",
    "\n",
    "# Add descriptive stats to both city and cluster maps. Title pages too. Add CRP and Bank Logo\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "if __name__ ==  '__main__': \n",
    "    if create_the_ppt:\n",
    "        prs = Presentation() #Create a new presentation. Add the map \n",
    "        add_title_plus_big_picture_slide(city,tables, prs, width, height)\n",
    "        var_list = [f\"{city}_City Boundary_map\",f\"{city}_map_land_surface_temperature\" , f\"{city}_map_pop_density\", f\"{city}_map_tree_cover_esa\", f\"{city}_map_tree_cover_landsat\", f\"{city}_map_pop_density\", f\"{city}_map_esri_dynamic_worldcover\" ]\n",
    "        for var_name in var_list:\n",
    "                if \"clusters\" not in var_name:\n",
    "                    var_name= os.path.basename(var_name)\n",
    "                    # print(var_name)\n",
    "                    var_name=var_name.replace(f'{city}_map_','')\n",
    "                    var_name=var_name.replace('.png','')\n",
    "                    # print(var_name)\n",
    "                    file_list = glob.glob(os.path.join( maps, f\"{city}_map_{var_name}*.png\"))\n",
    "                    for file_path in file_list:\n",
    "                        file_name= os.path.basename(file_path)\n",
    "                        add_city_map_slide( img_path=file_path ,tables=tables, prs=prs, width=width, height=height, var=var_name)\n",
    "                        file_path=os.path.join( maps, f\"{city}_map_clusters_{var_name}.png\")\n",
    "                        var_name= os.path.basename(file_path)\n",
    "                        # print(var_name)\n",
    "                        var_name=var_name.replace(f'{city}_map_clusters_','')\n",
    "                        var_name=var_name.replace('.png','')\n",
    "                        # print(var_name)\n",
    "                        file_name= os.path.basename(file_path)\n",
    "                        add_cluster_map_slide( img_path=file_path ,tables=tables, prs=prs, width=width, height=height, var=var_name)\n",
    "\n",
    "        \n",
    "        var_list = glob.glob(os.path.join( maps, f\"map_*.png\"))\n",
    "        for var_name in var_list:\n",
    "            if \"clusters\" not in var_name:\n",
    "                var_name= os.path.basename(var_name)\n",
    "                # print(var_name)\n",
    "                var_name=var_name.replace('map_','')\n",
    "                var_name=var_name.replace('.png','')\n",
    "                print(var_name)\n",
    "                var=  var_name #\"CSDI ssp370 2021 2040\"    \n",
    "                file_list = glob.glob(os.path.join( maps, f\"map_{var_name}*.png\"))\n",
    "                for file_path in file_list:\n",
    "                    file_name= os.path.basename(file_path)\n",
    "                    add_city_map_slide( img_path=file_path ,tables=tables, prs=prs, width=width, height=height, var=var_name)\n",
    "                    file_path=os.path.join( maps, f\"map_clusters_{var_name}.png\")\n",
    "                    var_name= os.path.basename(file_path)\n",
    "                    # print(var_name)\n",
    "                    var_name=var_name.replace('map_clusters_','')\n",
    "                    var_name=var_name.replace('.png','')\n",
    "                    print(var_name)\n",
    "                    file_name= os.path.basename(file_path)\n",
    "                    add_cluster_map_slide( img_path=file_path ,tables=tables, prs=prs, width=width, height=height, var=var_name)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='section5'></a>\n",
    "<h5><center> <font color='cyan'> Section 5: Population density and tree cover Sub-city level</font>   </center></h5>\n",
    "\n",
    "**Playground**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "uh",
   "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.8.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}