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

Land-surface temperature (LST) is how hot the ground itself is (roofs, roads, bare soil), measured from orbit by reading the heat the surface gives off as infrared light. It is not the air temperature on a thermometer. It’s the skin temperature of the city. During a heatwave some blocks run far hotter than others, and this question maps which blocks bake the most, and roughly how many people live in them.

Not the same as q11, the urban heat-island trend. That question asks how the whole-city heat island has grown over decades (one number for the city, tracked over years). This question is about one place during one event: during a single heatwave, which blocks inside the city run hottest, and whether those hotspots line up with low greenery and dense buildings. It maps who’s most exposed, not a long-term trend.

Which data to use, and how

Use this dataset: MOD11A2 (short name MOD11A2). It’s NASA’s MODIS satellite land-surface temperature, an 8-day average at 1 km resolution, going back to 2000. “8-day average” means each value is the surface temperature blended over 8 days, which fills in the days clouds blocked the view. 1 km means one number covers a roughly 1 km square, a small cluster of city blocks. Start with this layer; the greenery and population layers get added on top later.

Then do four steps (the runnable code further down does all of this; this is just the idea):

1. Pick your city and the heatwave dates. A small latitude/longitude box around the city, and the 8-day window covering the hot spell.
2. Get a temperature for every 1 km square. Read the LST value out of each MODIS file, throw out the “no data” markers, and convert to real degrees. The raw numbers are stored scaled, so you multiply by 0.02 to get Kelvin, then subtract 273.15 for Celsius.
3. Rank the squares and measure the spread. Find the hottest 10% of squares, and the gap between the hottest and coolest. That gap inside one city is often far bigger than any decade-long citywide trend.
4. Overlay greenery and people. Pair each square with vegetation (greener = cooler) and with a free population layer to estimate how many people sit in the hottest squares, then name the district. (Verified for the default Delhi box: ~25.5 million people, hottest district Central.)

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can see each step work.

What you can find out

• Which neighborhoods are hottest during the heatwave. Rank the 1 km squares for the 8-day window covering the event, instead of just the citywide average.
• The intra-city heat spread, hottest-minus-coolest neighborhood on the same dates. For Delhi’s May 2024 heatwave this was +6.8 °C, often far larger than any decade-long trend.
• The greenery overlay. Pair each hot square with MOD13Q1 NDVI (a greenness index, a stand-in for tree cover) to see whether the hottest blocks are also the least green. For Delhi the link was a strong −0.62 (more green, cooler ground).
• The built-up overlay. VNP46A2 nightlights (how brightly a place glows at night, a stand-in for how built-up it is) flag dense developed zones, so you can find the hot + low-green + bright-built triple overlap.
• A heat-vulnerability shortlist. The squares in that overlap are candidate blocks for shade, cool roofs, or tree-canopy investment.
• Roughly how many people live in the hot zone. Overlay free WorldPop population (1 km) to count people inside the hottest squares, then name the district with geoBoundaries.

What it can’t tell you

• Exactly who lives there (census-grade). WorldPop is a modeled estimate, not a head count. It gives a sound order-of-magnitude exposure number but not households, income, age, or health. For equity targeting you’d join a real census or income layer on top.
• The air temperature people actually feel. LST is the heat of the ground and roofs. The 2 m air temperature on a thermometer can differ by 5–15 °C, and indoors differs again. Treat the map as a relative hotspot ranking, not absolute exposure.
• Block-by-block detail finer than ~1 km. One MOD11A2 square covers several real city blocks, so it blurs them together. For street-level hotspots use sharper sensors (ECOSTRESS at 70 m or Landsat TIRS at 100 m).
• The single hottest day. MOD11A2 is an 8-day average. For a one-day peak use the daily products instead (MOD11A1 / MYD11A1).
• What caused it, or whether a fix worked. The green-vs-hot link is a correlation, not proof. Showing that planting trees actually cooled a block needs before/after measurements on the ground.

Gotchas to watch for

• Ground heat is not air heat. Surface (skin) temperature can run 5–15 °C above the air people breathe, so a reading isn’t telling you “it was 50 °C outside.” Lean on the map to rank hot blocks rather than to pin down absolute exposure.
• It’s an 8-day average, not a day. MOD11A2 blends 8 days, which captures the heatwave’s sustained pattern but smooths out a single peak afternoon. When you need one specific date, switch to the daily product (MOD11A1 / MYD11A1).
• 1 km blurs neighborhoods together. Several blocks share one square, so the fine differences between them vanish. For street resolution, add a sharper sensor like ECOSTRESS or Landsat TIRS.
• Filter out cloud-spoiled pixels first. Cloud-contaminated temperatures can look perfectly reasonable while quietly skewing your ranking. Each file ships a quality layer (the QC_Day quality flag), so drop the low-quality squares before averaging.
• Greenery and lights are stand-ins, not the real thing. NDVI isn’t tree canopy and night brightness isn’t building density. They’re proxies, and squeezing the 250 m / 500 m layers onto the 1 km grid mixes neighbors together, so treat them as supporting evidence rather than measurements.
• There are no people in these satellite files. Every “who lives there” claim comes from a separate population or census layer joined in afterward, each with its own age and accuracy.

The real-data code

The Run it cell above runs the method on synthetic data (no login). Below is the whole recipe against the real archive: search MOD11A2, read and reproject the LST grid, resample NDVI onto it, rank the hottest squares, correlate with greenery, overlay free population, and plot. It needs a free Earthdata Login.

import earthaccess, numpy as np, requests, rasterio
from pyhdf.SD import SD, SDC
from rasterio.warp import reproject, Resampling, calculate_default_transform
from rasterio.windows import from_bounds
from rasterio.transform import from_origin
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

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

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

def read_sds(path, name, scale):
    """Read one SDS from an HDF4 granule -> scaled float array (fill -> NaN)."""
    hdf = SD(path, SDC.READ); sds = hdf.select(name)
    raw, fill = sds.get().astype("float64"), sds.attributes().get("_FillValue", 0)
    out = np.where(raw == fill, np.nan, raw * scale)
    sds.endaccess(); hdf.end()
    return out

# 2. LST (scale 0.02 -> Kelvin) + QC: drop pixels whose QC_Day flags bad quality.
def lst_celsius(path):
    lst = read_sds(path, "LST_Day_1km", 0.02)
    qc = read_sds(path, "QC_Day", 1.0)            # bits 0-1: 00 = good quality
    return np.where((np.nan_to_num(qc).astype("uint8") & 0b11) == 0, lst - 273.15, np.nan)

lst_c = np.nanmean([lst_celsius(f) for f in lst_files], axis=0)   # one C value per 1 km square

# 3. Reproject the MODIS sinusoidal grid to a regular lat/lon grid over the AOI.
SINU = ("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m")
ny, nx = lst_c.shape
# MODIS tile h24v06 sinusoidal extent for these rows/cols (gridded 1 km):
sin_tf = from_origin(7783653.638, 3335851.559, 926.625, 926.625)
dst_tf, dw, dh = calculate_default_transform(SINU, "EPSG:4326", nx, ny,
                                             *rasterio.transform.array_bounds(ny, nx, sin_tf))
lst_ll = np.full((dh, dw), np.nan)
reproject(lst_c, lst_ll, src_transform=sin_tf, src_crs=SINU,
          dst_transform=dst_tf, dst_crs="EPSG:4326", resampling=Resampling.bilinear)

# 4. NDVI greenness (250 m, scale 0.0001) resampled ONTO the LST lat/lon grid.
ndvi_grans = earthaccess.search_data(short_name="MOD13Q1", version="061",
                                     bounding_box=city_aoi, temporal=heatwave)
ndvi_files = earthaccess.download(ndvi_grans, "./mod13q1")
ndvi_sin = np.nanmean([read_sds(f, "250m 16 days NDVI", 0.0001) for f in ndvi_files], axis=0)
ny2, nx2 = ndvi_sin.shape
ndvi_tf = from_origin(7783653.638, 3335851.559, 231.656, 231.656)   # 250 m grid, same tile origin
ndvi = np.full_like(lst_ll, np.nan)
reproject(ndvi_sin, ndvi, src_transform=ndvi_tf, src_crs=SINU,
          dst_transform=dst_tf, dst_crs="EPSG:4326", resampling=Resampling.average)

# 5. Rank the hottest squares + the greenery correlation (the justice fingerprint).
hot_thr = np.nanpercentile(lst_ll, 90)
hot_mask = lst_ll >= hot_thr
spread = np.nanmax(lst_ll) - np.nanmin(lst_ll)
ok = ~np.isnan(lst_ll) & ~np.isnan(ndvi)
r = np.corrcoef(lst_ll[ok], ndvi[ok])[0, 1]

# 6. Who lives there — free WorldPop 1 km population (no NASA login), windowed to the AOI.
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
people = np.nansum(pop)

print(f"Intra-city heat spread: +{spread:.1f} C  |  hottest decile starts at {hot_thr:.1f} C")
print(f"LST-NDVI correlation: {r:+.2f} (negative = greener is cooler)  |  {people:,.0f} people in AOI")

# Plot: the LST map with the hottest-decile neighborhoods outlined.
plt.imshow(lst_ll, cmap="inferno"); plt.colorbar(label="LST (C)")
plt.contour(hot_mask, levels=[0.5], colors="cyan", linewidths=0.8)
plt.title(f"Delhi LST, May 2024 heatwave  (spread +{spread:.1f} C, r={r:+.2f})")
plt.tight_layout(); plt.show()

# Cross-check: LST is the skin temperature of the ground, not the 2 m air people feel
# (can differ 5-15 C) — confirm the ranking against a ground station or ECOSTRESS 70 m.

Where the data comes from

All the NASA layers come through one free login (Earthdata Login). The MODIS land-surface temperature and greenery (MOD11A2, MOD13Q1) live in NASA’s LP DAAC archive, and the VIIRS Black Marble nightlights (VNP46A2) in LAADS DAAC, and earthaccess reaches both with the same login. The “who lives there” layers (WorldPop population, geoBoundaries district names) come from outside NASA, are free, and get joined in on your own machine.

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/