Heat-scan-workflow_nis.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: Population density and tree cover Sub-city level. NA ](#section4) \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": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages\n",
    "import geopandas as gpd , os , time , math\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Initialize ee and import shapefile and select city\n",
    "ee.Initialize()\n",
    "run_vito_analysis=True\n",
    "# run_vito_analysis=False\n",
    "run_tree_and_pop_analysis=False\n",
    "# run_tree_and_pop_analysis=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",
    "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",
    "\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. Each municipality has its own assembly (elected every four years in local elections), a municipal president, public service property and a budget. Municipalities usually have more than 10,000 inhabitants.[2]\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": 38,
   "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)"
   ]
  },
  {
   "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": 39,
   "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",
    "\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",
    "        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",
    "\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)"
   ]
  },
  {
   "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": 40,
   "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=2\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=3\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",
    "\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",
    "          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",
    "               park[\"index\"]=index\n",
    "               park[\"points\"]=points\n",
    "               park[\"rank\"]=rank\n",
    "               park[\"value\"]=value\n",
    "               park.to_csv(f\"{tables}/{index}_{category}_park_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               education[\"index\"]=index\n",
    "               education[\"points\"]=points\n",
    "               education[\"rank\"]=rank\n",
    "               education[\"value\"]=value\n",
    "               education.to_csv(f\"{tables}/{index}_{category}_education_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               industry[\"index\"]=index\n",
    "               industry[\"points\"]=points\n",
    "               industry[\"rank\"]=rank\n",
    "               industry[\"value\"]=value\n",
    "               industry.to_csv(f\"{tables}/{index}_{category}_industry_df.csv\")\n",
    "          except:\n",
    "               pass\n",
    "\n",
    "          try:\n",
    "               building[\"index\"]=index\n",
    "               building[\"points\"]=points\n",
    "               building[\"rank\"]=rank\n",
    "               building[\"value\"]=value\n",
    "               building.to_csv(f\"{tables}/{index}_{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}/{index}_{category}_infra_network_stats.csv\")\n",
    "          except:\n",
    "               print(f\"no network for {index}\")\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",
    "# place = 'Nis, Serbia'\n",
    "# get_osm_network_and_travel_time=get_osm_network()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "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": 42,
   "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",
    "        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",
    "        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",
    "        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",
    "        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",
    "        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": 45,
   "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\"Vito variable: {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\"Vito variable: {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",
    "\n",
    "\n",
    "# def get_vito_raster_path(raster_dirs , input_raster_dir, gdf_city):\n",
    "#     fc_filtered  , region , hexbins_projected= create_hexbins_and_geom_features(gdf_city, hexbin_res=8, WGS84=WGS84)\n",
    "#     for subdir in raster_dirs:\n",
    "#         extension = '.tif'\n",
    "#         for root, dirs_list, files_list in os.walk( os.path.join(input_raster_dir, subdir)):\n",
    "#             for file_name in files_list:\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\"Vito variable: {file_name} \\n across {city}\"\n",
    "#                     scheme=\"quantiles\"\n",
    "#                     categorical=False\n",
    "#                     map_output= f'{maps}/map_{file_name}.png'\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\"Vito variable: {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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "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",
    "        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: Population density and tree cover Sub-city level</font>   </center></h5>\n",
    "\n",
    "**Playground**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assemble the pieces**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pptx import Presentation\n",
    "# Create a new presentation\n",
    "presentation = Presentation()\n",
    "# Adding a Slide with Layout 3 \n",
    "# In python-pptx, slide layouts are numbered from 0. Slide layout 3 is the ‘Title and Content’ layout. Let’s add a slide with this layout.\n",
    "slide_layout = presentation.slide_layouts[3]\n",
    "slide = presentation.slides.add_slide(slide_layout)\n",
    "# Adding a Title\n",
    "# The ‘Title and Content’ layout has two placeholders: one for the title and one for the content. Let’s add a title to our slide.\n",
    "title = slide.shapes.title\n",
    "title.text = \"Dataframe in PPT\"\n",
    "\n",
    "# Adding a Dataframe\n",
    "# Now, let’s add a dataframe to our slide. We’ll create a dataframe using pandas and then add it to the slide.\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "# Create a dataframe\n",
    "df = pd.DataFrame({\n",
    "    'Name': ['Alice', 'Bob', 'Charlie'],\n",
    "    'Age': [25, 30, 35],\n",
    "    'Occupation': ['Engineer', 'Doctor', 'Teacher']\n",
    "})\n",
    "\n",
    "# Add the dataframe to the slide\n",
    "# content = slide.placeholders[1]\n",
    "# # table = content.insert_table(df.shape[0] + 1, df.shape[1]).table\n",
    "\n",
    "# Add the column names\n",
    "# for i, column_name in enumerate(df.columns):\n",
    "#     table.cell(0, i).text = column_name\n",
    "\n",
    "# # Add the data\n",
    "# for i in range(df.shape[0]):\n",
    "#     for j in range(df.shape[1]):\n",
    "#         table.cell(i + 1, j).text = str(df.iloc[i, j])\n",
    "# Adding Text\n",
    "# Finally, let’s add some text to our slide. We’ll add a new slide with layout 3 and then add some text to it.\n",
    "\n",
    "# Add a new slide with layout 3\n",
    "slide = presentation.slides.add_slide(slide_layout)\n",
    "\n",
    "# Add a title\n",
    "title = slide.shapes.title\n",
    "title.text = \"Text in PPT\"\n",
    "\n",
    "# Add some text\n",
    "content = slide.placeholders[1]\n",
    "content.text = \"This is some text.\"\n",
    "# Saving the Presentation \n",
    "# Once we’re done adding content to our slides, we can save the presentation.\n",
    "\n",
    "presentation.save(f'{tables}/presentation.pptx')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df.plot()\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",
    "# df=df.drop('index_right', axis=1)\n",
    "# db_pop= db_pop.to_crs(crs)\n",
    "# df1 = gpd.sjoin(df, db_pop, how='left', op='contains')\n",
    "# df1=df1.drop('index_right', axis=1)\n",
    "# db_landsat= db_landsat.to_crs(crs)\n",
    "# df2 = gpd.sjoin(df1, db_landsat, how='left', op='contains')\n",
    "# df2=df2.drop('index_right', axis=1)\n",
    "# db_esa= db_esa.to_crs(crs)\n",
    "# df3 = gpd.sjoin(df2, db_esa, how='left', op='contains')\n",
    "# df3\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",
    "# # raster_filename= \"tree_cover_esa\"\n",
    "# cluster_boundaries_lst=cluster_boundaries_lst.to_crs(crs)\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 = gpd.sjoin(cluster_boundaries_lst, db_var, how='left', op='contains')\n",
    "# df = gpd.sjoin(df, db_var, how='left', op='contains')\n",
    "# df\n",
    "# df_all=df_all.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"geometry\"])\n",
    "# df_all\n",
    "\n",
    "# maindf = cluster_boundaries_lst \n",
    "# # inputframes = [db_esri_dyna, db_pop_densi, db_esa_tree_, db_land_tree]\n",
    "# inputframes = [ db_var]\n",
    "# resultframes = [] #An empty list to hold the results\n",
    "\n",
    "# for frame in inputframes:\n",
    "#     frame =frame.to_crs(crs)\n",
    "#     frame['geometry'] = frame['geometry'].centroid\n",
    "#     result= overlay(maindf, frame, how=\"intersection\" , keep_geom_type=False)\n",
    "#     result=result.sort_values('area_for_dupes', ascending=False).drop_duplicates([\"geometry\"])\n",
    "#     # result = gpd.overlay(df1=maindf, df=frame, how=\"intersection\")\n",
    "#     result['area'] = result.geometry.area\n",
    "#     resultframes.append(result)\n",
    "# resultframes\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "For the zonal h3 have the routing and correlations. Apply h3 to the above just in case.\n",
    "Multivariate hotspots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# var dataset = ee.ImageCollection(\"ESA/WorldCover/v100\").first();\n",
    "\n",
    "# var areaImage = ee.Image.pixelArea().addBands(\n",
    "#       dataset)\n",
    "\n",
    "# var areas = areaImage.reduceRegion({\n",
    "#       reducer: ee.Reducer.sum().group({\n",
    "#       groupField: 1,\n",
    "#       groupName: 'landcover_class',\n",
    "#     }),\n",
    "#     geometry: area,\n",
    "#     scale: 10,\n",
    "#     maxPixels: 1e13\n",
    "#     }); \n",
    "\n",
    "\n",
    "# import numpy as np\n",
    "# import geopandas as gpd\n",
    "\n",
    "# df = gpd.read_file('/home/bera/geodata/Rectangle_with_hole.shp')\n",
    "# g = [i for i in df.geometry]\n",
    "\n",
    "# all_coords = []\n",
    "# for b in g[0].boundary: # for first feature/row\n",
    "#     coords = np.dstack(b.coords.xy).tolist()\n",
    "#     all_coords.append(*coords)                 \n",
    "\n",
    "# all_coords\n",
    "\n",
    "# import math\n",
    "# from shapely.affinity import rotate\n",
    "# from shapely.geometry import LineString, Point\n",
    "\n",
    "# start = Point(0, 0)\n",
    "# length = 1\n",
    "# angle = math.pi / 3\n",
    "\n",
    "# end = Point(start.x + length, start.y)\n",
    "# line = LineString([start, end])\n",
    "# line = rotate(line, angle, origin=start, use_radians=True)\n",
    "# print(line)\n",
    "\n",
    "# import geopandas as gpd\n",
    "# from shapely.geometry import Point, LineString\n",
    "# import matplotlib.pyplot as plt\n",
    "# pA1 = Point(1.5, 1.75)\n",
    "# # cluster_boundaries\n",
    "# # handle all points and create relating lines\n",
    "# pA1 = Point(1.5, 1.75)\n",
    "# pA1 = Point(cluster_boundaries.centroid[0])\n",
    "# pA2 = Point(2, 2)\n",
    "# line_A = LineString([[pA1.x, pA1.y], [pA2.x, pA2.y]])\n",
    "# pB1 = Point(3.0, 2.0)\n",
    "# pB2 = Point(3, 4)\n",
    "# line_B = LineString([[pB1.x, pB1.y], [pB2.x, pB2.y]])\n",
    "# pC1 = Point(2.5, 1.25)\n",
    "# pC2 = Point(1, 1)\n",
    "# line_C = LineString([[pC1.x, pC1.y], [pC2.x, pC2.y]])\n",
    "\n",
    "# # create a geodataframe,\n",
    "# # assigning the column containing `LineString` as its geometry\n",
    "# pts_and_lines = gpd.GeoDataFrame([['A', pA1, pA2, 14, line_A],\n",
    "#             ['B', pB1, pB2, 18, line_B],\n",
    "#             ['C', pC1, pC2, 19, line_C]],\n",
    "#             columns=['id', 'beg_pt', 'end_pt', 'value', 'LineString_obj'], \n",
    "#             geometry='LineString_obj')  # declare LineString (last column) as the `geometry`\n",
    "\n",
    "\n",
    "# # make a plot of the geodataframe obtained\n",
    "# f, ax = plt.subplots(1, figsize = [4, 4])\n",
    "# pts_and_lines.plot(ax=ax, column = 'value');\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysal\n",
    "import numpy as np\n",
    "from pysal.model.spreg import ols, ml_error, ml_lag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pysal package contains many sample data files that can be used to demonstrate the package's abilities. For this example we will be analyzing Columbus home values with relation to income and crime by neighborhood. First we will run an ordinary least squares linear regression model to analyze the relationship between these variables.\n",
    "\n",
    "This first section of code will read in the home values (dependent variable) into an array y and the income and crime values (independent variables) into a two dimmensional array X."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = pysal.lib.io.open(pysal.lib.examples.get_path(\"columbus.dbf\"),'r')\n",
    "y = np.array(f.by_col['HOVAL'])\n",
    "\n",
    "y.shape = (len(y),1)\n",
    "X= []\n",
    "X.append(f.by_col['INC'])\n",
    "X.append(f.by_col['CRIME'])\n",
    "X = np.array(X).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have stored the values we are analyzing we can perform ordinary least squares (OLS) regression. This is done with the pysal.spreg module. Our instance of OLS, named ls, has many useful tools for reviewing the results of our test. In this case, we will use ls.summary to obtain of full summary of the results, but for more specific results you can look at some of the other options on the pysal.spreg page linked above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "REGRESSION\n",
      "----------\n",
      "SUMMARY OF OUTPUT: ORDINARY LEAST SQUARES\n",
      "-----------------------------------------\n",
      "Data set            :    Columbus\n",
      "Weights matrix      :        None\n",
      "Dependent Variable  :    home val                Number of Observations:          49\n",
      "Mean dependent var  :     38.4362                Number of Variables   :           3\n",
      "S.D. dependent var  :     18.4661                Degrees of Freedom    :          46\n",
      "R-squared           :      0.3495\n",
      "Adjusted R-squared  :      0.3212\n",
      "Sum squared residual:   10647.015                F-statistic           :     12.3582\n",
      "Sigma-square        :     231.457                Prob(F-statistic)     :   5.064e-05\n",
      "S.E. of regression  :      15.214                Log likelihood        :    -201.368\n",
      "Sigma-square ML     :     217.286                Akaike info criterion :     408.735\n",
      "S.E of regression ML:     14.7406                Schwarz criterion     :     414.411\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "            Variable     Coefficient       Std.Error     t-Statistic     Probability\n",
      "------------------------------------------------------------------------------------\n",
      "            CONSTANT      46.4281827      13.1917570       3.5194844       0.0009867\n",
      "              Income       0.6289840       0.5359104       1.1736736       0.2465669\n",
      "               Crime      -0.4848885       0.1826729      -2.6544086       0.0108745\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "REGRESSION DIAGNOSTICS\n",
      "MULTICOLLINEARITY CONDITION NUMBER           12.538\n",
      "\n",
      "TEST ON NORMALITY OF ERRORS\n",
      "TEST                             DF        VALUE           PROB\n",
      "Jarque-Bera                       2          39.706           0.0000\n",
      "\n",
      "DIAGNOSTICS FOR HETEROSKEDASTICITY\n",
      "RANDOM COEFFICIENTS\n",
      "TEST                             DF        VALUE           PROB\n",
      "Breusch-Pagan test                2           5.767           0.0559\n",
      "Koenker-Bassett test              2           2.270           0.3214\n",
      "================================ END OF REPORT =====================================\n"
     ]
    }
   ],
   "source": [
    "\n",
    "ls = ols.OLS(y, X, name_y = 'home val', name_x = ['Income', 'Crime'], name_ds = 'Columbus')\n",
    "print(ls.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate spatial autocorrelation in the residuals with Moran's I test. Our first step in that process is to create a spatial weights matrix. PySAL's example data has a GAL file that we can read in directly to create this matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = pysal.lib.io.open(pysal.lib.examples.get_path(\"columbus.gal\")).read()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can pass this weights matrix to the pysal.Moran function, along with our model residuals (ls.u).\n",
    "\n",
    "The observed value for I is much higher than the value we would expect if there was no spatial dependence. The p-value is the probability that we would observe the value of I that we did (or one greater) if there were no spatial dependence.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed I: 0.17131015816921516 \n",
      "Expected I: -0.020833333333333332 \n",
      "   p-value: 0.018930427521201776\n"
     ]
    }
   ],
   "source": [
    "\n",
    "mi = pysal.explore.esda.Moran(ls.u, w, two_tailed=False)\n",
    "print('Observed I:', mi.I, '\\nExpected I:', mi.EI, '\\n   p-value:', mi.p_norm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a spatial regression model to account for spatial non-independence. The spreg module has several different functions for creating a spatial regression model. In this example, we will use a spatial error model, but the implementation of a spatial lag model is similar.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "shapes (49,1) and (49,1) not aligned: 1 (dim 1) != 49 (dim 0)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32mc:\\Users\\Aziz\\Dropbox\\CRP\\UHI\\scripts\\Heat-scan-workflow_nis.ipynb Cell 33\u001b[0m line \u001b[0;36m1\n\u001b[1;32m----> <a href='vscode-notebook-cell:/c%3A/Users/Aziz/Dropbox/CRP/UHI/scripts/Heat-scan-workflow_nis.ipynb#X45sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m spat_err \u001b[39m=\u001b[39m ml_error\u001b[39m.\u001b[39;49mML_Error(y, X, w, \n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/Aziz/Dropbox/CRP/UHI/scripts/Heat-scan-workflow_nis.ipynb#X45sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m                              name_y\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mhome value\u001b[39;49m\u001b[39m'\u001b[39;49m, name_x\u001b[39m=\u001b[39;49m[\u001b[39m'\u001b[39;49m\u001b[39mincome\u001b[39;49m\u001b[39m'\u001b[39;49m,\u001b[39m'\u001b[39;49m\u001b[39mcrime\u001b[39;49m\u001b[39m'\u001b[39;49m], \n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/Aziz/Dropbox/CRP/UHI/scripts/Heat-scan-workflow_nis.ipynb#X45sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m                              name_w\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mcolumbus.gal\u001b[39;49m\u001b[39m'\u001b[39;49m, name_ds\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mcolumbus\u001b[39;49m\u001b[39m'\u001b[39;49m)\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/Aziz/Dropbox/CRP/UHI/scripts/Heat-scan-workflow_nis.ipynb#X45sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m \u001b[39mprint\u001b[39m(spat_err\u001b[39m.\u001b[39msummary)\n",
      "File \u001b[1;32mc:\\Users\\Aziz\\anaconda3\\envs\\uh\\lib\\site-packages\\pysal\\model\\spreg\\ml_error.py:466\u001b[0m, in \u001b[0;36mML_Error.__init__\u001b[1;34m(self, y, x, w, method, epsilon, spat_diag, vm, name_y, name_x, name_w, name_ds)\u001b[0m\n\u001b[0;32m    464\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39maic \u001b[39m=\u001b[39m DIAG\u001b[39m.\u001b[39makaike(reg\u001b[39m=\u001b[39m\u001b[39mself\u001b[39m)\n\u001b[0;32m    465\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mschwarz \u001b[39m=\u001b[39m DIAG\u001b[39m.\u001b[39mschwarz(reg\u001b[39m=\u001b[39m\u001b[39mself\u001b[39m)\n\u001b[1;32m--> 466\u001b[0m SUMMARY\u001b[39m.\u001b[39;49mML_Error(reg\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m, w\u001b[39m=\u001b[39;49mw, vm\u001b[39m=\u001b[39;49mvm, spat_diag\u001b[39m=\u001b[39;49mspat_diag)\n",
      "File \u001b[1;32mc:\\Users\\Aziz\\anaconda3\\envs\\uh\\lib\\site-packages\\pysal\\model\\spreg\\summary_output.py:243\u001b[0m, in \u001b[0;36mML_Error\u001b[1;34m(reg, w, vm, spat_diag, regimes)\u001b[0m\n\u001b[0;32m    241\u001b[0m reg\u001b[39m.\u001b[39m__summary \u001b[39m=\u001b[39m {}\n\u001b[0;32m    242\u001b[0m \u001b[39m# compute diagnostics and organize summary output\u001b[39;00m\n\u001b[1;32m--> 243\u001b[0m beta_diag(reg, robust\u001b[39m=\u001b[39;49m\u001b[39mNone\u001b[39;49;00m)\n\u001b[0;32m    244\u001b[0m reg\u001b[39m.\u001b[39m__summary[\u001b[39m'\u001b[39m\u001b[39msummary_r2\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m%-20s\u001b[39;00m\u001b[39m:\u001b[39m\u001b[39m%12.3f\u001b[39;00m\u001b[39m                \u001b[39m\u001b[39m%-22s\u001b[39;00m\u001b[39m:\u001b[39m\u001b[39m%12.3f\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39m%\u001b[39m (\n\u001b[0;32m    245\u001b[0m     \u001b[39m'\u001b[39m\u001b[39mSigma-square ML\u001b[39m\u001b[39m'\u001b[39m, reg\u001b[39m.\u001b[39msig2, \u001b[39m'\u001b[39m\u001b[39mLog likelihood\u001b[39m\u001b[39m'\u001b[39m, reg\u001b[39m.\u001b[39mlogll)\n\u001b[0;32m    246\u001b[0m reg\u001b[39m.\u001b[39m__summary[\u001b[39m'\u001b[39m\u001b[39msummary_r2\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m%-20s\u001b[39;00m\u001b[39m:\u001b[39m\u001b[39m%12.3f\u001b[39;00m\u001b[39m                \u001b[39m\u001b[39m%-22s\u001b[39;00m\u001b[39m:\u001b[39m\u001b[39m%12.3f\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39m%\u001b[39m (\n\u001b[0;32m    247\u001b[0m     \u001b[39m'\u001b[39m\u001b[39mS.E of regression\u001b[39m\u001b[39m'\u001b[39m, np\u001b[39m.\u001b[39msqrt(reg\u001b[39m.\u001b[39msig2), \u001b[39m'\u001b[39m\u001b[39mAkaike info criterion\u001b[39m\u001b[39m'\u001b[39m, reg\u001b[39m.\u001b[39maic)\n",
      "File \u001b[1;32mc:\\Users\\Aziz\\anaconda3\\envs\\uh\\lib\\site-packages\\pysal\\model\\spreg\\summary_output.py:714\u001b[0m, in \u001b[0;36mbeta_diag\u001b[1;34m(reg, robust)\u001b[0m\n\u001b[0;32m    712\u001b[0m reg\u001b[39m.\u001b[39mstd_err \u001b[39m=\u001b[39m diagnostics\u001b[39m.\u001b[39mse_betas(reg)\n\u001b[0;32m    713\u001b[0m reg\u001b[39m.\u001b[39mz_stat \u001b[39m=\u001b[39m diagnostics\u001b[39m.\u001b[39mt_stat(reg, z_stat\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n\u001b[1;32m--> 714\u001b[0m reg\u001b[39m.\u001b[39mpr2 \u001b[39m=\u001b[39m diagnostics_tsls\u001b[39m.\u001b[39;49mpr2_aspatial(reg)\n\u001b[0;32m    715\u001b[0m \u001b[39m# organize summary output\u001b[39;00m\n\u001b[0;32m    716\u001b[0m reg\u001b[39m.\u001b[39m__summary[\u001b[39m'\u001b[39m\u001b[39msummary_std_err\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m robust\n",
      "File \u001b[1;32mc:\\Users\\Aziz\\anaconda3\\envs\\uh\\lib\\site-packages\\pysal\\model\\spreg\\diagnostics_tsls.py:217\u001b[0m, in \u001b[0;36mpr2_aspatial\u001b[1;34m(tslsreg)\u001b[0m\n\u001b[0;32m    215\u001b[0m y \u001b[39m=\u001b[39m tslsreg\u001b[39m.\u001b[39my\n\u001b[0;32m    216\u001b[0m predy \u001b[39m=\u001b[39m tslsreg\u001b[39m.\u001b[39mpredy\n\u001b[1;32m--> 217\u001b[0m pr \u001b[39m=\u001b[39m pearsonr(y, predy)[\u001b[39m0\u001b[39m]\n\u001b[0;32m    218\u001b[0m pr2_result \u001b[39m=\u001b[39m \u001b[39mfloat\u001b[39m(pr \u001b[39m*\u001b[39m\u001b[39m*\u001b[39m \u001b[39m2\u001b[39m)\n\u001b[0;32m    219\u001b[0m \u001b[39mreturn\u001b[39;00m pr2_result\n",
      "File \u001b[1;32mc:\\Users\\Aziz\\anaconda3\\envs\\uh\\lib\\site-packages\\scipy\\stats\\_stats_py.py:4463\u001b[0m, in \u001b[0;36mpearsonr\u001b[1;34m(x, y, alternative)\u001b[0m\n\u001b[0;32m   4459\u001b[0m     msg \u001b[39m=\u001b[39m (\u001b[39m\"\u001b[39m\u001b[39mAn input array is nearly constant; the computed \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   4460\u001b[0m            \u001b[39m\"\u001b[39m\u001b[39mcorrelation coefficient may be inaccurate.\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m   4461\u001b[0m     warnings\u001b[39m.\u001b[39mwarn(stats\u001b[39m.\u001b[39mNearConstantInputWarning(msg))\n\u001b[1;32m-> 4463\u001b[0m r \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39;49mdot(xm\u001b[39m/\u001b[39;49mnormxm, ym\u001b[39m/\u001b[39;49mnormym)\n\u001b[0;32m   4465\u001b[0m \u001b[39m# Presumably, if abs(r) > 1, then it is only some small artifact of\u001b[39;00m\n\u001b[0;32m   4466\u001b[0m \u001b[39m# floating point arithmetic.\u001b[39;00m\n\u001b[0;32m   4467\u001b[0m r \u001b[39m=\u001b[39m \u001b[39mmax\u001b[39m(\u001b[39mmin\u001b[39m(r, \u001b[39m1.0\u001b[39m), \u001b[39m-\u001b[39m\u001b[39m1.0\u001b[39m)\n",
      "File \u001b[1;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36mdot\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: shapes (49,1) and (49,1) not aligned: 1 (dim 1) != 49 (dim 0)"
     ]
    }
   ],
   "source": [
    "\n",
    "spat_err = ml_error.ML_Error(y, X, w, \n",
    "                             name_y='home value', name_x=['income','crime'], \n",
    "                             name_w='columbus.gal', name_ds='columbus')\n",
    "print(spat_err.summary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "\u001b[1;31mUnable to start Kernel 'uh (Python 3.8.17)' due to a timeout waiting for the ports to get used. \n",
      "\u001b[1;31mView Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "from shapely.geometry import MultiPoint\n",
    "points2 = polys.copy()\n",
    "points2.geometry = points2.geometry.apply(lambda x: MultiPoint(list(x.exterior.coords)))"
   ]
  }
 ],
 "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
}