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

A wildfire leaves a dark, charred scar that you can see from space — but the more useful question is what happens next. Does the green come back? How quickly? NASA’s HLSL30 (Harmonized Landsat) gives you 30 m surface reflectance every few days, so you can watch a single burn scar over months and years and measure its recovery with the Normalized Burn Ratio (NBR) — a simple index that drops sharply right after a fire and climbs back as vegetation regrows.

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), an increase of +0.076 NBR — more than a doubling. That is real, measurable greening of the burn scar. NBR is a proxy for recovery, not proof that the original forest is back: it rises whether the regrowth is mature trees, young shrubs, or grass.

What you can answer

• How burned the land was right after the fire — a single post-fire NBR scene, lower = more severe burn
• Whether the scar is greening back — compare NBR (or NDVI) from a post-fire scene to one a year or more later over the same area; a rising index means returning vegetation
• 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 is 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 you can NOT answer with these datasets alone

• What kind of vegetation came back — NBR/NDVI rise for grass, shrubs, or trees alike; only field surveys or hyperspectral data tell you whether the forest (not just ground cover) is returning
• Whether the ecosystem has truly recovered — green pixels are not the same as restored soil, wildlife, or canopy structure
• Why recovery is fast or slow — rainfall, replanting, soil burn severity, and erosion all matter and aren’t in the reflectance alone
• What happened on cloudy or smoky dates — optical sensors can’t see through cloud or thick smoke, so you must drop those scenes (and gaps can hide short-term changes)
• Sub-30 m detail — a single scorched tree or a narrow regrowth strip is below HLS resolution

Code template (Python, cloud-direct)

Verified locally. HLSL30 is delivered as per-band Cloud-Optimized GeoTIFFs (COGs) on LP DAAC. 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 os, re, warnings, earthaccess, rioxarray, numpy as np
warnings.filterwarnings("ignore")

# load Earthdata creds from .env without `source` (passwords can break the shell)
for line in open(".env"):
    m = re.match(r'\s*(?:export\s+)?([A-Z0-9_]+)\s*=\s*(.*)\s*$', line)
    if m: os.environ.setdefault(m.group(1), m.group(2).strip().strip('"').strip("'"))
earthaccess.login(strategy="environment")   # free Earthdata Login

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

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, clip to the scar
    da = rioxarray.open_rasterio(fs.open(url), masked=True).squeeze()
    return da.rio.reproject("EPSG:4326").rio.clip_box(W, S, E, N)

def mean_nbr(g):
    b05 = read_clip(band_url(g, "B05")).astype("float32")           # NIR
    b07 = read_clip(band_url(g, "B07")).astype("float32").rio.reproject_match(b05)  # SWIR2
    fm  = read_clip(band_url(g, "Fmask")).rio.reproject_match(b05)
    nir, swir2 = b05 * 0.0001, b07 * 0.0001
    nbr = ((nir - swir2) / (nir + swir2)).values.astype("float32")
    cloud = ((fm.values.astype("int32") & 0b10) > 0) | ((fm.values.astype("int32") & 0b1000) > 0)
    nbr[cloud] = np.nan
    nbr[~np.isfinite(nbr)] = np.nan
    return float(np.nanmean(nbr))

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))
e, l = mean_nbr(early[0]), mean_nbr(late[0])
print(f"NBR  ~11 mo after fire: {e:.3f}   ~3 yr after: {l:.3f}   recovery: {l-e:+.3f}")
# -> NBR  ~11 mo after fire: 0.056   ~3 yr after: 0.132   recovery: +0.076