Is the shoreline near me advancing or retreating?

The shoreline is just the line where land meets water in a satellite photo. Trace that line in pictures taken years apart and you can see whether the beach has grown out toward the sea (accretion) or been eaten back toward the land (erosion). The catch is the tide. It slides the waterline back and forth all day, so you have to compare like with like.

Which data to use, and how

Use this dataset: HLS Sentinel-2 (short name HLSS30). “HLS” stands for Harmonized Landsat Sentinel, a tidied-up satellite image archive at 30 m per pixel covering 2015 → today. It’s the easiest to start with: one free login, cloud-screened, on a clean grid. Want a longer history? Landsat Collection 2 reaches back to 1982 (also 30 m, but a separate USGS login), and Sentinel-2 L2A gives a sharper 10 m waterline for narrow beaches. Add those only once HLS is working.

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

1. Pick your coast and your years. A small latitude/longitude box along the shore, plus a few date windows (say, one per year you care about).
2. Find the waterline in each picture. Compute a water index (a simple formula that makes water bright and land dark; here MNDWI = (green − shortwave-infrared) / (green + shortwave-infrared)), then cut it at a threshold to get a clean land/water boundary line.
3. Draw lines out from the coast. Lay down evenly spaced transects (straight lines running from land out to sea) and measure where each year’s waterline crosses each one.
4. Measure the movement. How far the line shifted per year along each transect. Seaward means the beach is growing, landward means it’s eroding, and a map of these rates shows where change is fastest.

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 watch each step work.

What you can find out

• Where the waterline sat in each year for any stretch of coast, from the water-index boundary.
• Erosion vs growth, in metres per year along the cross-shore transect lines. The standard approach is called DSAS, the USGS Digital Shoreline Analysis System.
• Where the coast is changing fastest: a hotspot map of the worst retreat and the strongest growth.
• The decades-long trend by extending back to Landsat Collection 2 (1982 onward).
• Season vs long-term change by comparing scenes before and after the monsoon.

What it can’t tell you

• Real erosion vs the tide just being higher that day. The waterline you trace is wherever the sea happened to be at the moment of the photo. Without correcting for tide (or only comparing photos taken at a similar tide), a higher tide can fake erosion that isn’t there.
• Change smaller than a pixel. A 30 m pixel (HLS/Landsat) or 10 m pixel (Sentinel-2) sets the floor. Change narrower than that gets lost in blurry “is it land or water?” pixels. Use 10 m Sentinel-2 for thin beaches.
• How much sand moved (volume). This measures where the line is in 2-D, not the 3-D shape of the beach. For volume you need height data (lidar or a survey).
• Why the coast is changing. Imagery shows the change, not the cause. It can’t tell a sea wall from a port from rising seas from sand mining. That needs extra information.

Gotchas to watch for

• The tide and waves move the line too. The optical waterline is just where the sea edge was at photo-time, so a high-tide photo can look like erosion sitting next to a low-tide one. Compare photos taken at a similar tide stage, and treat very small per-year changes as uncertain.
• The “where’s the water?” cutoff isn’t one fixed number. MNDWI > 0 is a common starting point, but muddy plumes (turbid water), sea foam, wet sand, and dark rooftops all shift it, and a wrong cutoff puts the line in the wrong place. Rather than hard-coding the threshold, check it on one clear scene per coast.
• You can’t see change finer than a pixel. Pixel size sets the floor on what’s detectable (~30 m HLS/Landsat, ~10 m Sentinel-2), so narrow beaches may show nothing real. Reach for 10 m Sentinel-2 L2A on thin beaches.
• Clouds and the monsoon hide the coast. The east coast is cloudy through the northeast-monsoon months, which leaves gaps and partial scenes. Prefer dry-season images and stack several together into a median composite to fill the holes.
• Check before you trust the numbers. Satellite shoreline rates carry several metres of uncertainty, so small rates may just be noise. Cross-check against surveyed shorelines or DSAS reference data before you report anything.

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 HLS scenes, open the green and SWIR bands, build a dry-season MNDWI composite per epoch, threshold to a waterline, measure how far it moved along cross-shore transects, print the rate, and map it. It needs a free Earthdata Login.

import earthaccess, numpy as np, rasterio
import matplotlib.pyplot as plt
from rasterio.warp import transform as warp_xy

earthaccess.login(strategy="netrc")

# A short, straight stretch of the Andhra–Tamil Nadu east coast as (W, S, E, N).
aoi = (80.18, 13.20, 80.34, 13.40)
epochs = ["2015", "2020", "2025"]   # one dry-season window per epoch

# 1. For each epoch: search HLS Sentinel-2, open green (B03) + SWIR1 (B11),
#    median-composite over the window, and threshold MNDWI to a land/water mask.
def water_mask(year):
    grans = earthaccess.search_data(
        short_name="HLSS30", bounding_box=aoi,
        temporal=(f"{year}-01-01", f"{year}-04-30"), cloud_cover=(0, 10))
    g_stack, s_stack, ref = [], [], None
    for g in grans:
        assets = {l["Name"][-3:]: l["URL"]
                  for l in g["umm"]["RelatedUrls"] if l.get("URL", "").endswith(".tif")}
        if "B03" not in assets or "B11" not in assets:
            continue
        with rasterio.open(assets["B03"]) as fg, rasterio.open(assets["B11"]) as fs:
            win = fg.window(*fg.bounds) if ref is None else ref
            g_stack.append(fg.read(1).astype("f4"))
            s_stack.append(fs.read(1).astype("f4"))
            if ref is None:
                ref, transform_, crs = True, fg.transform, fg.crs
    green = np.nanmedian(np.stack(g_stack), axis=0)   # composite suppresses clouds/waves
    swir = np.nanmedian(np.stack(s_stack), axis=0)
    mndwi = (green - swir) / (green + swir + 1e-9)
    return (mndwi > 0), transform_, crs            # MNDWI>0 = water (verify per coast)

masks = {y: water_mask(y) for y in epochs}
mask0, transform0, crs0 = masks[epochs[0]]
px = abs(transform0.a)                              # pixel size in metres (HLS = 30 m)

# 2. Cast shore-normal transects (here image rows), find the land->water crossing
#    column per row per epoch = the shoreline intercept, in metres from the box edge.
def intercepts(mask):
    out = []
    for row in mask:                               # first water pixel scanning seaward
        w = np.where(row)[0]
        out.append(w[0] * px if w.size else np.nan)
    return np.array(out)

pos = {y: intercepts(masks[y][0]) for y in epochs}

# 3. End-point rate (m/yr): seaward (+) = accretion, landward (-) = erosion.
yrs = int(epochs[-1]) - int(epochs[0])
rate = (pos[epochs[-1]] - pos[epochs[0]]) / yrs    # +ve col-shift = waterline moved seaward
rate = rate[np.isfinite(rate)]
mean_rate = np.nanmean(rate)
verdict = "accreting (growing)" if mean_rate > 0 else "eroding (retreating)"
print(f"Shoreline is {verdict}: {mean_rate:+.1f} m/yr mean "
      f"({np.nanpercentile(rate,10):+.1f} to {np.nanpercentile(rate,90):+.1f} m/yr along-shore)")

# 4. Map the per-transect change rate along the coast.
plt.figure(figsize=(7, 3))
plt.plot(rate, lw=1)
plt.axhline(0, c="k", lw=0.6); plt.axhline(px/yrs, ls=":", c="grey", label="±1 pixel/record")
plt.axhline(-px/yrs, ls=":", c="grey")
plt.ylabel("rate (m/yr)"); plt.xlabel("alongshore transect"); plt.legend()
plt.title(f"{epochs[0]}→{epochs[-1]} shoreline change"); plt.tight_layout(); plt.show()

# Before you trust it: rates within ±1 pixel/record are noise, and the waterline is the
# instantaneous tide line — cross-check against surveyed / DSAS shorelines at a matched tide stage.

Where the data comes from

HLS (HLSS30 and its Landsat twin HLSL30) comes from NASA’s LP DAAC archive through one free login (Earthdata Login). Going further back with Landsat Collection 2 L2 pulls from USGS, which uses separate credentials, and the sharper 10 m Sentinel-2 L2A comes from ESA/Copernicus. Each extra source you add is a new account, not a new password on the same one.

Sources

• HLS (HLSS30/HLSL30): https://lpdaac.usgs.gov/products/hlss30v002/
• Landsat Collection 2: https://www.usgs.gov/landsat-missions/landsat-collection-2
• DSAS (Digital Shoreline Analysis System): https://www.usgs.gov/centers/whcmsc/science/digital-shoreline-analysis-system-dsas
• Water indices (NDWI/MNDWI/AWEI) — McFeeters 1996; Xu 2006; Feyisa et al. 2014