Are the mangroves and coastal wetlands that protect my shore shrinking?

Mangroves are a quiet seawall in front of millions of people. These coastal forests break waves, trap mud, and hold the shoreline together. When they thin out or get cleared, the land behind them is more exposed to storms and erosion. From space you can measure how green that belt of forest is and watch it fill in or fade over the years. Healthy leaves reflect light in a particular way that a satellite can see.

Which data to use, and how

Use this dataset: HLSL30 (short name HLSL30), short for Harmonized Landsat. It’s a tidied-up, cloud-screened view of the land surface at 30 m detail, every few days, going back to 2013. It’s the only NASA dataset you need for a first answer here.

The trick is a number called NDVI (greenness index). Live leaves soak up red light but bounce back near-infrared (which we can’t see). NDVI compares those two bands: high means lush forest, near zero means bare ground, negative means open water. HLSL30 ships each band as a separate image file, so you do the math yourself.

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

1. Pick your coast. A small latitude/longitude box around your estuary, delta, or shoreline, plus a time window. A dry-season month is cleanest, since dry-season leaves are least confused by flooding or fresh rain.
2. Grab two bands. The red (B04) and near-infrared (B05) images for a low-cloud day, clipped to your box.
3. Compute NDVI for every pixel: (B05 − B04) / (B05 + B04). Take the median across your box for one summary number. Dense mangrove reads ~0.6–0.9; open channels read negative.
4. Compare across years. Repeat for the same dry-season window in an earlier year and subtract the two maps. A patch that drops from green to bare is a candidate for real mangrove loss.

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. For the Sundarbans (21.6–22.0 °N, 88.8–89.3 °E) on 4 Jan 2024, one HLSL30 scene gave a median NDVI of 0.64 across the box, with 57% of pixels above 0.6. That’s the signature of a dense mangrove canopy threaded by tidal channels. The open-water channels pull the lowest readings negative, exactly as you’d expect in a delta.

What you can find out

• How dense and green the mangrove canopy is right now. NDVI from HLSL30, at 30 m, anywhere on the coast.
• Whether the canopy is thinning or filling in over the years. Compute NDVI for the same dry-season window each year and difference them to see gains and losses.
• Where along the shore the change is concentrated. Map the NDVI difference to spot specific eroding fronts, cleared patches, or recovering replanting zones.
• Roughly how wide the protective green belt is. Keep only pixels above a greenness threshold (e.g. NDVI > 0.5) to outline the vegetated buffer between open water and inhabited land.
• Open-water vs. land context. Overlay JRC Global Surface Water to separate permanent channels and ponds from the wetland canopy, and Global Mangrove Watch for a known mangrove footprint.

What it can’t tell you

• Whether a green pixel is truly mangrove. NDVI just sees “leafy.” It can’t tell mangrove from rice, palm, or other coastal plants. To confirm species you need Global Mangrove Watch or local field knowledge.
• Seasonal browning vs. real loss. A single low-NDVI scene could be a dry spell, a cloud artifact, or just a high tide flooding the roots. Only a multi-year, same-season comparison separates a true loss from the natural cycle.
• Why it changed. Greenness shows the what, not the cause. A cyclone, clearing for ponds, rising salt, sinking land, or sea-level rise all look the same from above.
• Below-ground health. HLS sees only the canopy top, not the soil carbon or root mass that make mangroves so valuable for storing carbon and steadying the shore.
• Anything finer than 30 m. A thin fringe of mangrove a few metres wide hides inside one 30 m pixel. For narrow belts you need higher-resolution imagery.
• Marshes or seagrass below the waterline. Submerged or half-flooded wetlands barely show up in an optical land index, so treat water-edge pixels with caution.

Gotchas to watch for

• Tides change the picture. A mangrove flooded at high tide looks less green than the same forest at low tide, because the water hides the leaves. Compare a high-tide scene to a low-tide one and you can fake a “loss” that isn’t there. Prefer scenes near the same tidal stage, and lean on the median rather than any single pixel.
• Always compare the same season. Mangroves brown a little in the dry season and green up after rain, just like any forest, so comparing January to July can invent a change that’s really just the calendar. Always difference the same month across years (dry-season to dry-season).
• Clouds are sneaky. HLSL30 is already cloud-screened, but the coast is cloudy and a thin haze or a missed cloud edge can drag NDVI down. Pick a low-cloud scene and don’t trust one bad day; back it up with another date in the same window.
• Water reads negative, and that’s correct. Tidal channels and ponds give negative NDVI, which pulls the average down in a delta. A low average can just mean lots of water, not sick forest, so look at the fraction of pixels above 0.6 (canopy) instead of the average alone.

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 red + NIR bands of a low-cloud scene in two dry-season years, compute median NDVI and canopy fraction for each, difference them, and plot it. It needs a free Earthdata Login.

HLSL30 ships as Cloud-Optimized GeoTIFFs (one file per band). Open the red (B04) and near-infrared (B05) bands with rioxarray, clip to your bounding box before computing, and take the median NDVI. Dense mangrove reads ~0.6–0.9; open water reads negative.

import os, re, warnings, earthaccess, rioxarray, numpy as np
import matplotlib.pyplot as plt
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 = 88.8, 21.6, 89.3, 22.0          # your coast (Sundarbans)

# 1. One dry-season NDVI map per year: search HLSL30, open B04 (red) + B05 (NIR)
#    COGs from the first scene, clip to the box, and compute per-pixel NDVI.
def ndvi_map(year):
    res = earthaccess.search_data(short_name="HLSL30",
                                  temporal=(f"{year}-01-01", f"{year}-03-31"),
                                  bounding_box=(W, S, E, N), count=60)
    urls = res[0].data_links()
    b04 = next(u for u in urls if "B04" in u and u.endswith(".tif"))
    b05 = next(u for u in urls if "B05" in u and u.endswith(".tif"))
    red, nir = earthaccess.open([b04, b05])
    red = rioxarray.open_rasterio(red, masked=True).rio.clip_box(W, S, E, N, crs="EPSG:4326")
    nir = rioxarray.open_rasterio(nir, masked=True).rio.clip_box(W, S, E, N, crs="EPSG:4326")
    ndvi = (nir.values - red.values) / (nir.values + red.values)
    return np.where(np.isfinite(ndvi), ndvi, np.nan).squeeze()

old, new = ndvi_map(2018), ndvi_map(2024)    # same dry-season window, six years apart

# 2. Summary numbers: median NDVI and the fraction of dense-canopy pixels (>0.6),
#    which is more honest than the mean in a delta full of negative-NDVI water.
def summarize(a):
    v = a[np.isfinite(a)]
    return float(np.nanmedian(v)), float((v > 0.6).mean())

med_old, canopy_old = summarize(old)
med_new, canopy_new = summarize(new)

# 3. Difference the two same-season maps; a persistent drop is candidate loss.
diff = new - old
canopy_change = canopy_new - canopy_old
verdict = "shrinking" if canopy_change < -0.02 else "growing" if canopy_change > 0.02 else "stable"
print(f"Mangrove canopy is {verdict}: dense-canopy fraction {canopy_old:.2f} -> {canopy_new:.2f} "
      f"({canopy_change:+.2f}); median NDVI {med_old:.2f} -> {med_new:.2f}")

# 4. Map where it changed — red = greened up, blue = lost canopy.
plt.imshow(diff, cmap="RdBu", vmin=-0.4, vmax=0.4)
plt.colorbar(label="ΔNDVI (2024 − 2018)"); plt.title("Same-season NDVI change")
plt.tight_layout(); plt.show()

# Before you trust it: overlay Global Mangrove Watch so a "loss" pixel is real mangrove
# (not rice/palm), and check both scenes are near the same tidal stage (see the gotchas).

Where the data comes from

HLSL30 comes from NASA through one free login (Earthdata Login). It’s served by the LP DAAC, NASA’s land-data archive, which harmonizes the Landsat satellites into one consistent 30 m product. The two optional baselines are free and separate. JRC Global Surface Water (from the European Commission’s Joint Research Centre) gives you a ready-made map of where water permanently sits. Global Mangrove Watch gives you a known footprint of where mangroves actually are, useful for checking that your green pixels are really mangrove and not some other crop.