Heat Vulnerability Index in Manhattan

1. Imports, setting up, data sources, data preparation

import os
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point
import osmnx as ox
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

colors = ["#FFA500", "#FF4500"]  # Orange to red
orange_red_cmap = LinearSegmentedColormap.from_list("OrangeRed", colors)

# https://data.cityofnewyork.us/Health/Modified-Zip-Code-Tabulation-Areas-MODZCTA-/pri4-ifjk/about_data
# https://data.cityofnewyork.us/Health/Heat-Vulnerability-Index-Rankings/4mhf-duep
# https://data.cityofnewyork.us/Housing-Development/Building-Footprints-Deprecated-/nqwf-w8eh

df = pd.read_csv('Heat_Vulnerability_Index_Rankings_20250115.csv')
admin = ox.geocode_to_gdf('Manhattan')
admin.plot()
display(df.head(3))


gdf = gpd.read_file('Modified Zip Code Tabulation Areas (MODZCTA)_20250115.geojson')
gdf.modzcta = gdf.modzcta.astype(int)
gdf = gpd.GeoDataFrame(df.merge(gdf, left_on = 'ZIP Code Tabulation Area (ZCTA) 2020', right_on = 'modzcta'))
gdf.plot(cmap = orange_red_cmap)

footprints = gpd.read_file('Building Footprints_20250115.geojson')
footprints.heightroof = footprints.heightroof.astype(float)
print(len(footprints))
footprints = gpd.overlay(admin, footprints)
print(len(footprints))
footprints.head(10).keys()

footprints2 = footprints[
    (footprints['heightroof'] > 49) & 
    (footprints.is_valid) & 
    (footprints.geometry.type == 'Polygon')
]

footprints2.plot()

footprints3 = gpd.sjoin(footprints2, gdf)

# Normalize the 'Heat Vulnerability Index (HVI)' values to the range [0, 1]
hvi_min = footprints3['Heat Vulnerability Index (HVI)'].min()
hvi_max = footprints3['Heat Vulnerability Index (HVI)'].max()
footprints3['HVI_normalized'] = (footprints3['Heat Vulnerability Index (HVI)'] - hvi_min) / (hvi_max - hvi_min)
footprints3['colors'] = footprints3['HVI_normalized'].apply(lambda x: orange_red_cmap(x))
footprints3['colors'] = footprints3['colors'].apply(lambda rgba: tuple(int(255 * v) for v in rgba[:3]))
footprints3 = footprints3[['colors', 'geometry', 'heightroof']]

len(footprints3)

2. Creating the map

import pydeck as pdk


# Calculate min and max for 'heightroof' column
min_height = footprints3['heightroof'].min()
max_height = footprints3['heightroof'].max()

# Update the pydeck layer
layer = pdk.Layer(
    "PolygonLayer",
    data=footprints3,
    get_polygon="geometry.coordinates",
    extruded=True,
    get_elevation="heightroof",
    get_fill_color = "colors",
)

view_state = pdk.ViewState(
    latitude=40.75755,
    longitude=-73.9801,
    zoom=13,
    pitch=45, # controls the angle of the camera relative to the ground
    bearing = -5 # controls the rotation of the camera around the 
)

# Create the 3D map
r = pdk.Deck(layers=[layer], initial_view_state=view_state)
r.to_html('DAY25.html') 
r