How to Summarise Land Surface Temperature by Neighbourhood in Python
On this page
Problem statement
Heat maps are pixels; decisions are made by district. A municipality planning cool pavements, shade and tree planting needs to rank its districts, compare them with social data and track them over years. That is zonal statistics — and the choice of statistic matters more than it looks. A district's mean surface temperature is dominated by its large surfaces (a park, an airport, a rail yard); its 90th percentile describes its hot spots; the share of its area above a threshold describes how much of it is hot. The rankings differ between them.
Quick answer
Rasterise the districts onto the LST grid and compute several statistics per district:
import numpy as np, pandas as pd, osmnx as ox, rasterio
from rasterio import features
with rasterio.open("madrid_lst_jja2023.tif") as src: # from the composite guide
lst = src.read(1)
tf, crs, shape = src.transform, src.crs, (src.height, src.width)
ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
d = ox.features_from_place("Madrid, Spain", {"boundary": "administrative", "admin_level": "9"})
d = d[(d.admin_level == "9") & d.geom_type.isin(["Polygon", "MultiPolygon"])][["name", "geometry"]]
city = ox.geocode_to_gdf("Madrid, Spain").geometry.iloc[0]
d = d[d.representative_point().within(city)] # drop neighbouring towns' districts
d = d.to_crs(crs).reset_index(drop=True)
zone = features.rasterize(((g, i + 1) for i, g in enumerate(d.geometry)), out_shape=shape,
transform=tf, fill=0)
HOT = 45.0
rows = []
for i, name in enumerate(d.name):
v = lst[(zone == i + 1) & np.isfinite(lst)].astype(float)
if v.size:
rows.append({"district": name, "median": np.median(v), "p90": np.percentile(v, 90),
f"area > {HOT:.0f} °C": (v > HOT).mean(), "km²": v.size * 900 / 1e6})
t = pd.DataFrame(rows).set_index("district").sort_values("median", ascending=False)
print(f"{len(t)} districts")
print(t.round(2).head(6).to_string())
print("…")
print(t.round(2).tail(3).to_string())21 districts
median p90 area > 45 °C km²
district
Villa de Vallecas 47.78 50.91 0.82 51.47
Barajas 47.02 50.60 0.71 41.68
Fuencarral-El Pardo 46.84 50.05 0.74 237.78
Vicálvaro 45.86 49.52 0.67 34.85
Villaverde 45.63 48.08 0.62 20.16
Latina 44.97 50.67 0.50 25.41
…
median p90 area > 45 °C km²
district
Ciudad Lineal 41.68 43.86 0.02 11.41
Retiro 41.56 43.77 0.01 5.45
Chamartín 41.49 43.72 0.02 9.20
Three statistics per district, and a geography that surprises people who expect the centre to be hottest. At late morning the hottest districts are on the edge — Villa de Vallecas, Barajas, Fuencarral-El Pardo — with runways, dry farmland and new estates; the coolest are central Chamartín, Retiro and Ciudad Lineal, six degrees lower. By median and by hot-area share the order is largely the same; the 90th percentile picks out districts with small, very hot patches, such as Latina, sixth by median but among the highest by 90th percentile.
Districts on the raster's CRS, several statistics
1. Get districts in the raster's CRS
Official boundaries or OSM administrative polygons; reproject them to the LST grid's CRS before rasterising.
2. Rasterise once
An integer zone array lets every statistic be a grouped operation. Pixels whose centre falls in a district belong to it.
3. Compute several statistics
Median for the typical surface, a high percentile for hot spots, and the share of area above a threshold for extent. Add the count of valid pixels.
4. Keep water and parks in mind
Large parks and reservoirs lower a district's median; they are real parts of the district, but they may mask very hot built-up areas beside them. Consider statistics over built-up pixels only as a second view.
5. Rank with caution
Differences of a fraction of a degree between districts are within the noise of a composite. Group districts into bands rather than reading exact ranks.
Rankings that disagree, and built-up surfaces only
The data used on this page
The summer 2023 land surface temperature composite of Madrid from the composite guide: the median of every Landsat 8 and 9 Collection 2 Level-2 scene from June to August 2023 on the Microsoft Planetary Computer with more than a fifth of the area clear, on a 30 m UTM 30N grid, with the count of clear observations as a second band. The setup rebuilds it if madrid_lst_jja2023.tif is missing. OpenStreetMap comes through osmnx, via the maps.mail.ru Overpass mirror while overpass-api.de was refusing connections, and is read as the mirror answered on 9 October 2026.
import io
import os
import requests
import zipfile
import numpy as np
import planetary_computer
import pystac_client
import rasterio
from rasterio.transform import from_origin
from rasterio.warp import Resampling, reproject, transform_bounds
# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
reply = requests.get("https://data.spatialworkflow.io/osm-cache-lst-by-neighbourhood-python-2026-10-09.zip", timeout=300)
reply.raise_for_status()
with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
archive.extractall("cache")
if not os.path.exists("madrid_lst_jja2023.tif"): # the composite guide's code, without the prints
AOI = (-4.05, 40.20, -3.35, 40.65)
DST_CRS = "EPSG:32630"
x0, y0, x1, y1 = transform_bounds("EPSG:4326", DST_CRS, *AOI)
x0, y1 = np.floor(x0 / 30) * 30, np.ceil(y1 / 30) * 30
W, H = int(np.ceil((x1 - x0) / 30)), int(np.ceil((y1 - y0) / 30))
DST_TF = from_origin(x0, y1, 30, 30)
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
items = [i for i in cat.search(collections=["landsat-c2-l2"], bbox=list(AOI),
datetime="2023-06-01/2023-08-31").items()
if i.properties["platform"] in ("landsat-8", "landsat-9") and i.properties["eo:cloud_cover"] < 60]
def on_grid(item, asset):
out = np.zeros((H, W), dtype="uint16")
with rasterio.open(planetary_computer.sign(item.assets[asset].href)) as src:
reproject(rasterio.band(src, 1), out, dst_transform=DST_TF, dst_crs=DST_CRS,
resampling=Resampling.nearest, src_nodata=0 if asset == "lwir11" else 1, dst_nodata=0)
return out
layers = []
for it in sorted(items, key=lambda i: i.datetime):
dn, qa = on_grid(it, "lwir11"), on_grid(it, "qa_pixel")
t = np.where(dn == 0, np.nan, dn * 0.00341802 + 149.0 - 273.15)
t[(((qa >> 1) & 1) | ((qa >> 3) & 1) | ((qa >> 4) & 1)) == 1] = np.nan
if np.isfinite(t).mean() > 0.2:
layers.append(t)
stack = np.stack(layers)
profile = {"driver": "GTiff", "height": H, "width": W, "count": 2, "dtype": "float32",
"crs": DST_CRS, "transform": DST_TF, "nodata": np.nan, "compress": "deflate"}
with rasterio.open("madrid_lst_jja2023.tif", "w", **profile) as dst:
dst.write(np.nanmedian(stack, axis=0).astype("float32"), 1)
dst.write(np.isfinite(stack).sum(axis=0).astype("float32"), 2)
dst.set_band_description(1, "median LST °C, Landsat 8/9, Jun–Aug 2023")
dst.set_band_description(2, "clear observations")
with rasterio.open("madrid_lst_jja2023.tif") as src:
print(src.shape, src.crs, src.descriptions or "", f"median {np.nanmedian(src.read(1)):.1f} °C")
import osmnx as ox
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api" # overpass-api.de refused us(1682, 1994) EPSG:32630 ('median LST °C, Landsat 8/9, Jun–Aug 2023', 'clear observations') median 46.9 °C
Example 1 — rankings that disagree
ranks = t[["median", "p90", f"area > {HOT:.0f} °C"]].rank(ascending=False)
disagree = ranks.max(axis=1) - ranks.min(axis=1)
print(f"districts whose rank differs by 5+ places between statistics: {(disagree >= 5).sum()}")
print(ranks.loc[disagree.sort_values(ascending=False).index[:4]].astype(int).to_string())districts whose rank differs by 5+ places between statistics: 5
median p90 area > 45 °C
district
Centro 12 14 19
Tetuán 13 17 20
Moncloa-Aravaca 14 8 9
Moratalaz 18 13 13
Five districts move five or more places between statistics. A district can be middling by median and near the top by 90th percentile, as Moncloa-Aravaca is — one hot patch, an industrial estate, rail land or a stretch of dry grass, in otherwise ordinary fabric. Which ranking to publish depends on whether the question is typical exposure or worst spots.
Example 2 — built-up surfaces only
import pystac_client, planetary_computer
from rasterio.warp import reproject, Resampling
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
wc = np.zeros(shape, dtype="uint8")
for it in cat.search(collections=["esa-worldcover"], bbox=[-3.90, 40.30, -3.50, 40.60], datetime="2021").items():
with rasterio.open(planetary_computer.sign(it.assets["map"].href)) as src:
tmp = np.zeros(shape, dtype="uint8")
reproject(rasterio.band(src, 1), tmp, dst_transform=tf, dst_crs=crs, resampling=Resampling.mode)
wc = np.where(wc == 0, tmp, wc)
zid = {n: i + 1 for i, n in enumerate(d.name)}
t["built median"] = [float(np.nanmedian(lst[(zone == zid[n]) & (wc == 50)])) for n in t.index]
t["tree share"] = [((zone == zid[n]) & (wc == 10)).sum() / (zone == zid[n]).sum() for n in t.index]
print(t[["median", "built median", "tree share"]].round(2).sort_values("built median",
ascending=False).head(5).to_string())
print(f"correlation of tree share with district median LST: "
f"{np.corrcoef(t['tree share'], t['median'])[0, 1]:.2f}") median built median tree share
district
Barajas 47.02 46.87 0.09
Villa de Vallecas 47.78 45.61 0.03
Villaverde 45.63 45.45 0.10
Vicálvaro 45.86 45.00 0.04
Carabanchel 43.98 44.10 0.15
correlation of tree share with district median LST: -0.38
Restricting to built-up pixels removes the parks' influence and compares the fabric itself. Tree share relates moderately to the all-surface median (r = −0.38) — partly because trees are cool, partly because a park is a lot of tree share. On built-up pixels alone, Barajas moves ahead of Villa de Vallecas, whose median falls by 2.2 °C: much of Vallecas's heat is open land, not buildings.
Example 3 — export for a report
out = d.set_index("name").join(t, how="inner")
out.to_file("madrid_districts_lst.gpkg")
print(f"wrote {len(out)} districts with {len(t.columns)} heat statistics")wrote 21 districts with 6 heat statistics
A district GeoPackage with the statistics is what a planning team can join to their own social and health data.
Several statistics, hidden variation, relative thresholds
Why several statistics
A district is a mix of surfaces. The median describes its typical surface; a high percentile describes its hottest; the share above a threshold describes how much of it is hot. Programmes aimed at hot spots, at typical streets or at extent need different ones.
Why districts hide what matters
District averages mix a park with the blocks beside it. Heat risk is local: the block without trees next to the park is not helped by the park's contribution to the average. District statistics are for allocation and monitoring; targeting needs the pixel map or small areas.
Why thresholds should be relative for comparisons across years
An absolute threshold (45 °C) depends on the summer's weather. For trends, use thresholds relative to each year's citywide median, or relative composites.
Small districts, edge pixels and admin levels
- Small districts have few thermal pixels (100 m native); their statistics are noisy.
- Edge pixels straddle district boundaries; with 30 m resampled thermal data this blurs boundaries by about 100 m.
- OSM admin levels vary by country; Madrid's 21 districts are level 9 in OSM, its neighbourhoods level 10.
- osmnx tag dictionaries combine keys with OR; filter the admin level after downloading.
- Population weighting changes the ranking again — see the exposure guide.
Internal links
- How to calculate zonal statistics in Python — the general method
- How to estimate population exposure to heat — weighting by people
- How to build a summer LST composite — the input
- Heat exposure explained — from hot surfaces to health
- How to relate LST to vegetation — why some districts are cooler
- Modifiable areal unit problem explained — why the zones matter
FAQ
How do I calculate average surface temperature per district?
Rasterise the districts onto the LST grid and compute the median (and other statistics) of valid pixels per district.
Median or mean?
Median, plus a high percentile and a hot-area share. The mean is pulled by large parks and water.
Why do district rankings change with the statistic?
Districts differ in how their hot surfaces are distributed — a few very hot patches or widespread moderate heat.
Should parks be included?
Yes for describing the district; also compute built-up-only statistics to compare the fabric.
Can I compare districts across years?
Use relative statistics or the same composite method; absolute temperatures depend on the summer.
What threshold defines a hot surface?
There is no standard; choose one from the distribution (e.g. the citywide 90th percentile) and report it.