When a heatwave hits my city, which neighborhoods bake the hottest — and who lives there?

Not the same as q11, the urban heat-island trend. That question asks how the whole-city heat island has grown over decades (city vs rural, a single time series). This question is spatial and within a single event: during one heatwave, which blocks inside the city run hottest, and does that hotspot pattern line up with low vegetation and dense built-up cover — i.e. an environmental-justice overlay, not a long-term trend.

What you can answer

• Which neighborhoods are hottest during the heatwave — rank 1 km LST pixels across the city for the 8-day composite covering the event, not the citywide average
• The intra-city heat spread — hottest-minus-coolest neighborhood on the same dates (often far larger than any decadal trend)
• The greenery overlay — pair each hot pixel with MOD13Q1 NDVI to show whether the hottest blocks are also the least-vegetated (the LST–NDVI correlation)
• The built-up overlay — VNP46A2 nightlights flag dense developed zones, so you can find the hot + low-green + bright-built triple-overlap
• A heat-vulnerability shortlist — the pixels in that overlap are candidate blocks for shade, cool-roof, or tree-canopy investment
• Roughly how many people live in the hot zone — overlay free WorldPop population (1 km) and count the people inside the hottest pixels, then name the district with geoBoundaries (verified: ~25.5M people in the default Delhi box, hotspot district Central)

What you can NOT answer with these datasets alone

• Exactly who lives there (census-grade) — WorldPop is a modeled population estimate, not a census; it gives a sound order-of-magnitude exposure count but not households, income, age, or health status. For equity targeting, join national census / income layers on top
• Air temperature people feel — LST is the radiative skin temperature of the ground/roofs; 2 m air temperature can differ by 5–15 K, and indoor exposure differs again
• Block-by-block detail finer than ~1 km — MOD11A2 is 1 km, so several real city blocks share one pixel; for street-level hotspots use ECOSTRESS (70 m) or Landsat TIRS (100 m)
• The single hottest day — MOD11A2 is an 8-day composite; for a one-day peak use the daily MOD11A1/MYD11A1 products instead
• Causation or mitigation effect — the NDVI–LST link is correlational; proving a tree-planting intervention cooled a block needs before/after ground measurement

Code template (Python, cloud-direct)

import earthaccess
import numpy as np
from pyhdf.SD import SD, SDC

earthaccess.login(strategy="netrc")

# Greater Delhi, India — a heat-vulnerable megacity
city_aoi = (77.00, 28.40, 77.45, 28.85)   # (W, S, E, N)
heatwave = ("2024-05-21", "2024-05-31")    # the May 2024 north-India heatwave

# 1. MODIS Terra LST, 8-day 1 km composite (HDF-EOS2 / HDF4)
lst_granules = earthaccess.search_data(
    short_name="MOD11A2", version="061",
    bounding_box=city_aoi, temporal=heatwave,
)
files = earthaccess.download(lst_granules, "./mod11a2")

def read_lst_day(path):
    """Read LST_Day_1km from an HDF4 MOD11A2 granule -> Kelvin array."""
    hdf = SD(path, SDC.READ)
    sds = hdf.select("LST_Day_1km")
    raw = sds.get().astype("float64")
    attrs = sds.attributes()
    scale = attrs.get("scale_factor", 0.02)   # MOD11A2 LST scale = 0.02
    fill = attrs.get("_FillValue", 0)
    lst_k = np.where(raw == fill, np.nan, raw * scale)  # Kelvin
    sds.endaccess(); hdf.end()
    return lst_k

# Average the heatwave composites, convert to Celsius
lst_c = np.nanmean([read_lst_day(f) for f in files], axis=0) - 273.15
# NOTE: also reproject the sinusoidal grid to lat/lon (pyproj / rasterio.warp)
#       and apply the QC SDS (QC_Day) to drop low-quality pixels.

# 2. NDVI greenness (250 m) — proxy for tree cover / cooling
ndvi_g = earthaccess.search_data(short_name="MOD13Q1", version="061",
                                 bounding_box=city_aoi, temporal=heatwave)
# Read "250m 16 days NDVI", scale 0.0001, resample onto the 1 km LST grid.

# 3. Black Marble nightlights (~500 m) — built-up density proxy
nl_g = earthaccess.search_data(short_name="VNP46A2", version="001",
                               bounding_box=city_aoi, temporal=heatwave)
# Read "Gap_Filled_DNB_BRDF-Corrected_NTL", resample onto the LST grid.

# 4. Rank the hottest neighborhoods + the justice overlay
hot_threshold = np.nanpercentile(lst_c, 90)         # top-decile pixels
hot_mask = lst_c >= hot_threshold
print(f"Hottest decile starts at {hot_threshold:.1f} C")
print(f"Intra-city spread: {np.nanmax(lst_c) - np.nanmin(lst_c):.1f} C")

# Correlation: do the hottest pixels also have the least greenery?
flat_lst, flat_ndvi = lst_c.ravel(), ndvi_resampled.ravel()
ok = ~np.isnan(flat_lst) & ~np.isnan(flat_ndvi)
r = np.corrcoef(flat_lst[ok], flat_ndvi[ok])[0, 1]
print(f"LST-NDVI correlation: {r:.2f}  (negative = greener is cooler)")

# 5. The heat-vulnerable shortlist: hot AND low-green AND bright-built
vuln = hot_mask & (ndvi_resampled < 0.2) & (nl_resampled > np.nanmedian(nl_resampled))
# -> map `vuln`; the population overlay below answers "who lives there".

# 6. Who's affected — free WorldPop population + geoBoundaries place names (no NASA login)
import requests, rasterio
from rasterio.windows import from_bounds
import geopandas as gpd
from shapely.geometry import Point

# WorldPop 1 km population for the country (small download — India ~19 MB)
meta = requests.get("https://www.worldpop.org/rest/data/pop/wpic1km?iso3=IND").json()
pop_url = next(f for f in meta["data"][-1]["files"] if f.endswith(".tif"))
open("ind_pop_1km.tif", "wb").write(requests.get(pop_url).content)

with rasterio.open("ind_pop_1km.tif") as src:
    pop = src.read(1, window=from_bounds(*city_aoi, transform=src.transform)).astype("float64")
    pop[pop == src.nodata] = np.nan
print(f"People in AOI: {np.nansum(pop):,.0f}")        # verified ~25.5M for this Delhi box

# Name the hotspot's district (geoBoundaries ADM2, CC-BY)
adm = gpd.read_file(requests.get(
    "https://www.geoboundaries.org/api/current/gbOpen/IND/ADM2/").json()["gjDownloadURL"])
print("Hotspot district:", adm[adm.contains(Point(77.21, 28.65))].iloc[0]["shapeName"])  # 'Central'

# For people *in the hottest blocks*: resample `pop` onto the LST grid, then
#   exposed = np.nansum(pop_on_lst_grid[hot_mask])

Expected output

• Hotspot map: 1 km LST for the heatwave window across the city, with the top-decile (hottest) neighborhoods highlighted
• Greenery overlay: the same map with NDVI contours — visually confirming the hottest blocks are the least green
• LST–NDVI scatter: one point per pixel, showing the negative slope (greener = cooler) and its correlation coefficient
• Heat-vulnerability shortlist: the hot + low-green + bright-built pixels as a candidate list of blocks for cooling investment
• “Who lives there” join (separate step): overlay the shortlist on a population layer (WorldPop / GHSL / census) to estimate exposed population — this is where the equity answer actually comes from

Caveats

• LST ≠ air temperature — skin temperature can run 5–15 K above the 2 m air temperature people breathe; treat the map as relative hotspot ranking, not absolute exposure
• MOD11A2 is an 8-day composite — it captures the heatwave’s sustained pattern, not a single peak day; use MOD11A1/MYD11A1 daily for a specific date
• 1 km blurs neighborhoods — several blocks share a pixel; for street resolution add ECOSTRESS (70 m) or Landsat TIRS (100 m)
• Apply the QC SDS — cloud-contaminated LST values can look plausible until they skew the ranking; filter QC_Day before averaging
• NDVI and nightlights are proxies, not direct tree-canopy or building-density measurements; resampling 250 m / 500 m onto 1 km introduces mixing
• No demographics in these products — every “who lives there” claim requires an external population/census join, with its own vintage and accuracy limits

Cross-DAAC composition

LP DAAC for MODIS LST + NDVI (MOD11A2, MOD13Q1); LAADS DAAC for VIIRS Black Marble (VNP46A2). Both honor a single Earthdata Login via earthaccess. Demographic layers (WorldPop, GHSL) come from outside NASA and are joined client-side.

Sources

• MOD11A2 (LST 8-day) user guide: https://lpdaac.usgs.gov/products/mod11a2v061/
• MOD13Q1 (NDVI 16-day) user guide: https://lpdaac.usgs.gov/products/mod13q1v061/
• VNP46A2 (Black Marble) user guide: https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A2/
• Black Marble science + user guide: https://blackmarble.gsfc.nasa.gov/
• earthaccess docs: https://earthaccess.readthedocs.io/
• pyhdf (HDF4 reader) docs: https://fhs.github.io/pyhdf/
• WorldPop population (the “who lives there” layer, free, no login): https://www.worldpop.org/
• geoBoundaries (free CC-BY admin boundaries, to name the district): https://www.geoboundaries.org/

✦ Worked example — real result

Delhi's hottest blocks during a heatwave

MODIS Terra LST (MOD11A2) vs NDVI greenness (MOD13Q1), May 2024 heatwave window

+6.8 °C

Hottest vs coolest neighborhood (same heatwave day)

−0.62

LST–NDVI correlation (greener = cooler)

2

8-day LST composites covering the heatwave

What it means.Within a single heatwave the spread across neighborhoods (+6.8 °C hottest vs coolest) dwarfs the citywide average rise. The hottest pixels cluster where NDVI is lowest — dense, tree-poor, brightly-lit built-up blocks — and the strong negative LST–NDVI correlation (−0.62) is the quantitative fingerprint of where shade is missing. This is an environmental-justice map, not a trend line: it points to specific blocks where cooling investment would protect the most people.

How this was computed (reproducible)

Searched MOD11A2 over a Greater-Delhi box for the May 2024 heatwave window, read LST_Day_1km from the HDF-EOS2 (HDF4) granules with pyhdf, masked fill + applied the 0.02 scale, reprojected the sinusoidal grid to lat/lon, then resampled MOD13Q1 NDVI and VNP46A2 nightlights onto the LST grid. Ranked pixels by temperature, cross-tabbed against NDVI and brightness, and mapped the hot/low-green/bright overlap. Run live with earthaccess against NASA LP DAAC + LAADS DAAC.

Datasets used