After the fire, is the land actually recovering — and how fast?

A wildfire leaves a dark, charred scar on the land that you can see from space. The more useful question is what happens next: does the green come back, and how fast? Surface reflectance is just how bright the ground looks to a satellite in different colors of light, including colors your eye can’t see. Healthy plants and bare burned dirt reflect those colors very differently, so by watching the reflectance of a single burn scar over months and years you can measure whether it’s greening back.

Which data to use, and how

Use this dataset: HLSL30 (short name HLSL30), NASA’s Harmonized Landsat product. It gives you a fresh picture of the same patch of ground every few days at 30 m resolution (each pixel is about the size of a baseball infield), going back to 2013. Each scene arrives as separate files per color band, so you only download the bands you need.

The trick is a number called the Normalized Burn Ratio (NBR). It’s a one-line formula that combines two reflectance bands (near-infrared and shortwave-infrared). It drops sharply right after a fire and climbs back as plants regrow. For a much longer history, MOD13 NDVI is a coarse, decades-long greenness record you can pair in later.

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

1. Pick your burn scar (a small latitude/longitude box around it) and two dates: one soon after the fire, one a year or more later.
2. Grab two bands. For each date, read just the near-infrared (B05) and shortwave-infrared (B07) bands inside your box.
3. Throw out the clouds. Use the scene’s quality band (Fmask) to blank out cloudy pixels, since a satellite can’t see ground through cloud.
4. Compute NBR and compare. Apply NBR = (B05 − B07) / (B05 + B07), average it over the scar for each date, and subtract. Up means the scar is greening back; flat means recovery has stalled.

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.

Verified locally. Over the 2021 Dixie Fire scar in California (40.0–40.4 °N), HLSL30 NBR, computed as (B05 − B07) / (B05 + B07) and averaged over the area after masking clouds, rose from 0.056 on 6 Jul 2022 (about 11 months after the fire) to 0.132 on 3 Jul 2024 (about 3 years after). That’s an increase of +0.076 NBR, more than a doubling. Real, measurable greening of the burn scar.

What you can find out

• How burned the land was right after the fire. A single post-fire NBR scene; lower means a more severe burn.
• Whether the scar is greening back. Compare NBR (or the greenness index NDVI) from a post-fire scene to one a year or more later over the same area. A rising index means returning plants.
• How fast recovery is going. Build a time series of NBR across many low-cloud scenes and watch the slope. Some scars rebound in a couple of seasons; others stall for years.
• Where recovery is uneven. Because HLS pixels are 30 m, you can map which parts of the scar are bouncing back and which are still bare.
• Long-term context. Pair with the coarse, decades-long MOD13 NDVI record to see how this fire’s recovery compares to the area’s pre-fire baseline.

What it can’t tell you

• What kind of plants came back. NBR and NDVI rise the same way for grass, shrubs, or trees, so a rising number is a proxy for recovery, not proof the original forest is back. Only field surveys or richer (hyperspectral) data can tell you whether the actual forest returned, not just ground cover.
• Whether the ecosystem has truly recovered. Green pixels are not the same as restored soil, wildlife, or the layered structure of a real canopy.
• Why recovery is fast or slow. Rainfall, replanting, how deeply the soil burned, and erosion all matter, and none of them are in the reflectance alone.
• What happened on cloudy or smoky days. Optical satellites can’t see through cloud or thick smoke, so you have to drop those scenes, and the gaps can hide short-term changes.
• Anything smaller than 30 m. A single scorched tree or a narrow strip of regrowth is below the pixel size, so it disappears into its neighbors.

Gotchas to watch for

• Clouds and smoke block the view. On a cloudy date the satellite sees cloud tops, not the ground. The fix: every scene ships a quality band (Fmask) that flags cloud and cloud-shadow pixels; blank those out before you average, and pick low-cloud dates to begin with.
• Compare like with like, same season. Plants are naturally greener in spring/summer and browner in winter, so always compare similar times of year across years (e.g., early July vs early July), never July vs January, or seasonal swing will masquerade as recovery.
• Raw band values need scaling. The files store reflectance as big integers to save space; you multiply by a tiny scale factor (0.0001) to get the real 0–1 reflectance before computing NBR. Skip it and the math still runs but the numbers are meaningless.
• Bands can sit on slightly different grids. The near-infrared and shortwave-infrared bands must line up pixel-for-pixel before you subtract them, so the code snaps one onto the other’s grid (reproject_match). Otherwise you’d be comparing offset locations.
• NBR is a proxy, not proof. It rises whether the regrowth is mature trees, young shrubs, or grass. Keep that honesty when you report a result.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search HLSL30, open the NIR/SWIR COG bands, mask clouds with Fmask, compute NBR over the scar, compare two dates, and plot it. It needs a free Earthdata Login.

HLSL30 is delivered as per-band Cloud-Optimized GeoTIFFs (COGs) on LP DAAC, so you can read just the slice you need straight from the cloud without downloading whole files. Open the NIR (B05) and SWIR-2 (B07) bands with rioxarray, clip to your burn scar before reading, mask clouds with the Fmask band, then compute NBR = (B05 − B07) / (B05 + B07). Compare an early post-fire scene to a later one.

import warnings, earthaccess, rioxarray, numpy as np
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")

earthaccess.login(strategy="netrc")          # free Earthdata Login

W, S, E, N = -121.4, 40.0, -121.0, 40.4      # Dixie Fire scar, California (2021)
fs = earthaccess.get_fsspec_https_session()

# 1. Locate one early post-fire scene and one ~2 years later, same season (early July).
early = earthaccess.search_data(short_name="HLSL30",
                                temporal=("2022-07-06", "2022-07-07"), bounding_box=(W, S, E, N))
late  = earthaccess.search_data(short_name="HLSL30",
                                temporal=("2024-07-03", "2024-07-04"), bounding_box=(W, S, E, N))

def band_url(g, b):
    return next(u for u in g.data_links() if u.endswith(f".{b}.tif"))

def read_clip(url):                           # read one COG band straight from the cloud, clip to scar
    da = rioxarray.open_rasterio(fs.open(url), masked=True).squeeze()
    return da.rio.reproject("EPSG:4326").rio.clip_box(W, S, E, N)

# 2. Open NIR (B05) + SWIR-2 (B07), snap onto a common grid, scale ints -> reflectance,
#    and 3. blank out cloud / cloud-shadow pixels flagged in Fmask before averaging.
def nbr_map(g):
    b05 = read_clip(band_url(g, "B05")).astype("float32")
    b07 = read_clip(band_url(g, "B07")).astype("float32").rio.reproject_match(b05)
    fm  = read_clip(band_url(g, "Fmask")).rio.reproject_match(b05)
    nir, swir2 = b05.values * 0.0001, b07.values * 0.0001
    nbr = (nir - swir2) / (nir + swir2)
    cloud = ((fm.values.astype("int32") & 0b10) > 0) | ((fm.values.astype("int32") & 0b1000) > 0)
    nbr[cloud | ~np.isfinite(nbr)] = np.nan
    return nbr

# 4. NBR = (B05 - B07)/(B05 + B07), averaged over the scar, then the before/after difference.
m_early, m_late = nbr_map(early[0]), nbr_map(late[0])
e, l = np.nanmean(m_early), np.nanmean(m_late)
verdict = "greening back" if l > e else "stalled / re-burning"
print(f"NBR ~11 mo after fire: {e:.3f}   ~3 yr after: {l:.3f}   recovery: {l-e:+.3f}  -> {verdict}")

# 5. Plot the two NBR maps side by side so you can see *where* the scar recovered.
fig, ax = plt.subplots(1, 2, figsize=(9, 4))
for a, m, t in zip(ax, (m_early, m_late), ("~11 mo after fire", "~3 yr after fire")):
    im = a.imshow(m, vmin=-0.2, vmax=0.5, cmap="RdYlGn"); a.set_title(t); a.axis("off")
fig.colorbar(im, ax=ax, label="NBR (higher = greener)", shrink=0.7)
plt.tight_layout(); plt.show()

# Cross-check: NBR rises the same for grass, shrubs, or trees — it's a proxy, not proof the
# forest returned. Confirm against a field survey or pair with MOD13 NDVI for the longer baseline.

Where the data comes from

Everything here is NASA through one free login (Earthdata Login). HLSL30 lives at the LP DAAC archive and blends NASA’s Landsat with Europe’s Sentinel-2 into one harmonized, every-few-days record. The optional long-history greenness record (MOD13 NDVI) is also NASA. The place name for your scar comes from geoBoundaries, a free open dataset of administrative boundaries — handy for labeling, separate from the satellite data.

Sources

• HLSL30 user guide: https://lpdaac.usgs.gov/products/hlsl30v002/
• MOD13 NDVI: https://lpdaac.usgs.gov/products/mod13a1v061/
• Normalized Burn Ratio (NBR) background: https://www.usgs.gov/landsat-missions/landsat-normalized-burn-ratio
• geoBoundaries (place names): https://www.geoboundaries.org/