How to Measure the Cooling Effect of Parks in Python

On this page

Problem statement

Parks are cooler than the streets around them, and the cooling reaches beyond their edges — the "park cool island". How far it reaches and how strong it is decide whether a new park helps a neighbourhood or only its own lawns, and they are measured the same way in most studies: surface temperature as a function of distance from the park boundary, inside and out.

The profile is simple to compute and easy to over-read. Landsat's thermal pixels are 100 m, so the first 100 m outside a park are partly the park itself; what surrounds a park (a river, a motorway, dense blocks) shapes its profile; and big parks behave differently from small ones.

Quick answer

Rasterise the parks, compute signed distance to their edge, and average LST in distance bands:

import numpy as np, geopandas as gpd, osmnx as ox, rasterio
from rasterio import features
from scipy import ndimage

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"
parks = ox.features_from_bbox((-3.78, 40.37, -3.62, 40.49), {"leisure": "park"})
parks = parks[parks.geom_type.isin(["Polygon", "MultiPolygon"])].to_crs(crs)
parks = parks[parks.area > 20_000]                               # at least 2 ha
print(f"{len(parks)} parks of 2 ha or more, e.g. {', '.join(parks.name.dropna().head(4))}")

inside = features.rasterize(((g, 1) for g in parks.geometry), out_shape=shape, transform=tf) == 1
d_out = ndimage.distance_transform_edt(~inside) * 30
d_in = ndimage.distance_transform_edt(inside) * 30
signed = np.where(inside, -d_in, d_out)

bands = [(-300, -150), (-150, -60), (-60, 0), (0, 60), (60, 150), (150, 300), (300, 600), (600, 1000)]
ok = np.isfinite(lst)
for lo, hi in bands:
    m = ok & (signed > lo) & (signed <= hi)
    print(f"{lo:>5} to {hi:>5} m: {m.sum():>7,} px, median LST {np.nanmedian(lst[m]):5.1f} °C")

185 parks of 2 ha or more, e.g. Parque del Oeste, Parque de Don David Pérez y Don Miguel Ángel Hernández, Jardines del Campo del Moro, Parque de Atenas
 -300 to  -150 m:   3,926 px, median LST  42.5 °C
 -150 to   -60 m:   9,326 px, median LST  42.1 °C
  -60 to     0 m:  10,783 px, median LST  42.1 °C
    0 to    60 m:  17,906 px, median LST  42.5 °C
   60 to   150 m:  30,340 px, median LST  42.8 °C
  150 to   300 m:  44,618 px, median LST  43.0 °C
  300 to   600 m:  66,654 px, median LST  43.1 °C
  600 to  1000 m:  59,021 px, median LST  43.5 °C

Negative distances are inside parks. The profile rises from 42.1 °C just inside the edges to 43.5 °C 600–1,000 m away, and keeps rising slowly across the whole range, so the "reach" of Madrid's parks is not a sharp edge. Two things stand out: the whole contrast is only 1.4 °C at late morning, and the deepest interiors are slightly warmer than the band just inside the edge. Both have the same cause, which the examples below find.

Choose parks, signed distance, bands, intensity and reach

1. Choose the parks

OSM leisure=park includes everything from pocket squares to large wooded parks. Set a minimum area — small parks are sub-pixel for Landsat thermal data.

2. Compute signed distance to the edge

Distance outside, negative distance inside. One layer then describes the whole transition.

3. Summarise by distance band, not pixel

Median LST per band, with pixel counts. Bands narrower than the thermal resolution (about 100 m) are not independent.

4. Estimate the cooling intensity and distance

Intensity: the difference between the park interior and the level the profile reaches far outside. Distance: where the profile flattens — often defined as the first band within some tolerance of the far value.

5. Separate large and small parks

Cooling depends on park size and on what the park contains — trees, watered lawn or dry grass. Profiles averaged over all parks mix very different behaviours.

Inside the park, the edge, and the distance at which the city takes over.

Cooling intensity and reach by park size

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. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 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 10 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-park-cooling-effect-python-2026-10-10.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 — cooling intensity and reach

far = np.nanmedian(lst[ok & (signed > 600) & (signed <= 1000)])
core = np.nanmedian(lst[ok & (signed < -150)])
print(f"park interior (>150 m inside): {core:.1f} °C; city 600–1,000 m away: {far:.1f} °C "
      f"→ cooling intensity {far - core:.1f} °C")
for step in range(0, 600, 60):
    m = ok & (signed > step) & (signed <= step + 60)
    if abs(np.nanmedian(lst[m]) - far) < 0.3:
        print(f"profile within 0.3 °C of the far value from ~{step} m outside the edge")
        break

park interior (>150 m inside): 42.9 °C; city 600–1,000 m away: 43.5 °C → cooling intensity 0.6 °C
profile within 0.3 °C of the far value from ~540 m outside the edge

The interior, more than 150 m inside, is only 0.6 °C cooler than the city 600–1,000 m away — less than the band just inside the edge, which is 1.4 °C cooler. The reach estimate, about 540 m here, depends on the tolerance and on the 100 m thermal pixels: the first ring outside a park always contains some park signal. Report the reach as "within about N hundred metres", not to the metre.

Example 2 — large against small parks

lab, n = ndimage.label(inside)
size_ha = np.bincount(lab.ravel()) * 900 / 1e4
for lo, hi in ((2, 10), (10, 50), (50, 10_000)):
    ids = np.where((size_ha >= lo) & (size_ha < hi))[0]
    ids = ids[ids > 0]
    m_in = ok & np.isin(lab, ids) & (signed < -30)
    if m_in.sum() > 50:
        print(f"parks {lo}–{hi} ha ({len(ids)}): interior median {np.nanmedian(lst[m_in]):5.1f} °C, "
              f"{far - np.nanmedian(lst[m_in]):.1f} °C below the city")

parks 2–10 ha (95): interior median  41.4 °C, 2.1 °C below the city
parks 10–50 ha (32): interior median  42.2 °C, 1.3 °C below the city
parks 50–10000 ha (12): interior median  42.6 °C, 0.9 °C below the city

It is the other way round from the textbook: the 95 small parks are 2.1 °C cooler inside than the city, the 12 largest only 0.9 °C. Small parks are close to the thermal pixel size and are measured partly as their surroundings, which should make them look warmer, not cooler — so the difference is real. What the parks contain is the likely explanation, in the next example.

Example 3 — what is inside matters

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.78, 40.37, -3.62, 40.49], 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)
for code, name in ((10, "tree cover"), (30, "grassland"), (60, "bare")):
    m = ok & inside & (wc == code)
    if m.sum() > 100:
        print(f"park pixels mapped as {name:<10}: {m.sum():>6,}, median LST {np.nanmedian(lst[m]):5.1f} °C")

park pixels mapped as tree cover: 14,164, median LST  40.8 °C
park pixels mapped as grassland :  6,382, median LST  44.8 °C
park pixels mapped as bare      :    318, median LST  43.0 °C

A park is not uniformly cool. Inside Madrid's parks, tree-covered pixels are at 40.8 °C and grassland at 44.8 °C — a degree hotter than the city 600–1,000 m away. In a Mediterranean summer unirrigated grass heats like a street, so a park's cooling depends on how much of it is trees rather than grass. Splitting the size classes by tree share is the check for whether that explains why the largest parks are warmest here. What cools is trees and water, not the word "park".

Median LST by distance from the nearest park edge, negative inside.

Why parks cool beyond their edges

Cool air over a park drifts into surrounding streets, especially on calm nights; shade from trees at the edge reaches the pavement. For surface temperature at Landsat's resolution, part of the apparent reach is also the thermal pixel footprint straddling the edge.

Why size and composition matter

Evaporative cooling scales with the area of transpiring vegetation, and shade with the area of canopy. A large wooded park has both; a paved square with a few trees has little of either. Averaging all "parks" together hides the difference that design decisions depend on.

Why surface and air cooling differ

Surface temperatures differ by many degrees between tree canopy and asphalt at midday; air temperatures differ by one or two, and most at night. Surface-based park cooling is an indicator of where air cooling is likely, not a measure of it.

Small parks were cooler inside than large ones; grass inside parks was hotter than the city.

Water in parks and overlapping influences

Water in parks (lakes, the Manzanares river) makes parks look cooler than their vegetation alone would.
Overlapping influences: near other parks, profiles never reach a "city" value.
Edge definition: OSM park boundaries include paved areas and car parks in some cases.
Direction: prevailing winds carry cool air further on the downwind side; profiles by direction reveal it.
Resolution: use bands of at least 60–100 m and treat the first outside band as mixed.

Internal links

Urban cooling and green space explained — the mechanisms
How to relate LST to vegetation — the pixel-level relationship
LST resolution misread — why the first 100 m are mixed
How to build a summer LST composite — the input
How to measure distance to the grid — the same distance transform elsewhere
How to calculate LST by neighbourhood — area summaries

FAQ

How do I measure park cooling from satellite data?

Compute distance to park edges, then median land surface temperature in distance bands inside and outside the parks.

How far does park cooling reach?

Typically a few hundred metres for large parks; estimate it from the profile and report it approximately.

Do small parks cool their surroundings?

They can: in Madrid the small parks were cooler inside than the large ones. They are near the thermal pixel size, so satellite estimates of their reach are uncertain.

Is surface cooling the same as air cooling?

No. Surface differences are much larger; air cooling is smaller and strongest at night.

Why are parts of my park hot?

Dry grass, bare ground and paving heat up in summer; tree cover and water stay cool.

Which parks cool most?

Ones with a lot of tree canopy and water, whatever their size; in a dry summer, open grass cools little.