How to measure shoreline change between two dates

Problem statement

Two shorelines, two dates, and the question is how far the coast moved. The answer is not the difference in their lengths, because the length of a coastline is not a number. It is a distance measured along a set of transects, perpendicular to the coast, at positions that stay fixed between surveys.

The measurement is straightforward. What makes shoreline change hard is that most of the apparent movement is usually not change: different tidal datums, different source scales, different extraction thresholds and different tide states at the moment of acquisition each produce a shift that looks exactly like erosion.

This guide builds the transect measurement and the checks that separate signal from method.

Quick answer

Cast transects from a fixed baseline and measure where each shoreline crosses:

import numpy as np, geopandas as gpd
from shapely.geometry import LineString, Point

def transects(baseline, spacing=50, length=500):
    out = []
    for d in np.arange(0, baseline.length, spacing):
        p = baseline.interpolate(d)
        q = baseline.interpolate(min(d + 1.0, baseline.length))
        nx, ny = (q.y - p.y), -(q.x - p.x)
        n = np.hypot(nx, ny) or 1.0
        out.append(LineString([(p.x - length * nx / n, p.y - length * ny / n),
                               (p.x + length * nx / n, p.y + length * ny / n)]))
    return gpd.GeoDataFrame({"id": range(len(out))}, geometry=out)

def intersect_distance(transect, shoreline, origin):
    hit = transect.intersection(shoreline)
    if hit.is_empty:
        return np.nan
    pts = [hit] if hit.geom_type == "Point" else list(hit.geoms)
    return min(origin.distance(p) for p in pts if p.geom_type == "Point")

Positive change is accretion, negative is erosion, and the sign convention has to be fixed by the transect's direction rather than inferred.

 year                    time  level
 2017 2017-07-04 15:35:56.461  0.350
 2021 2021-08-02 15:41:26.912 -0.065
 2023 2023-08-12 15:41:31.027  0.433
 2025 2025-08-11 15:41:27.553  0.188
{2017: '11.7 km', 2021: '10.4 km', 2023: '10.5 km', 2025: '10.4 km'}

Same datum, same CRS, a fixed baseline, transects

1. Put both shorelines in the same datum

This is the step that decides whether the answer means anything. A shoreline at MHHW and one at mean sea level differ by whatever horizontal distance the vertical offset covers — on a 1% beach with a 3.1 m range, 310 m. Convert both to a common datum before measuring.

2. Put both in the same projected CRS

Distances are in metres. Use a national grid or a local UTM zone, not degrees.

3. Build a baseline that does not move

A smoothed line offset landward of both shorelines, or an existing fixed feature. It is a coordinate system, not a shoreline, and it must be identical between every pair of dates you compare.

4. Cast transects at a fixed spacing

Fifty metres is a common choice for a beach study and 500 m for a regional one. Every transect keeps its identifier across all dates, which is what makes a time series possible.

5. Handle multiple intersections explicitly

A transect crossing a spit, a barrier island or a river mouth meets the shoreline several times. Taking the nearest crossing to the baseline is one defensible rule; taking the furthest is another. Pick one, and count how often it mattered.

6. Compute the rate, not just the difference

Change per year is comparable between pairs with different intervals. With three or more dates, fit a linear rate per transect and report its standard error.

7. Establish the uncertainty before interpreting anything

The positional uncertainty of each shoreline — georeferencing, digitising, pixel size, tidal correction — combines into a total. Change smaller than that is not detectable, and the fraction of transects above it is the honest headline.

Outer Cape Cod's shoreline, with an uncertainty budget

The data used on this page

The Atlantic shore of outer Cape Cod between Wellfleet and Truro, a coast that is known to be retreating. The shorelines are waterlines from Sentinel-2 L2A scenes — the least-cloudy summer scenes of 2017, 2021, 2023 and 2025, each chosen as the one whose water level at NOAA's Chatham gauge was closest to the 2017 scene's. Chatham is inside the harbour, not on the open beach, so the levels are the nearest measured stand-in. The beach slope comes from the USGS 3DEP DEM.

import numpy as np
import pandas as pd
import rasterio
import requests
import geopandas as gpd
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.features import shapes
from rasterio.warp import transform_bounds
from scipy import ndimage
from shapely.geometry import LineString, box, shape

BBOX = [-70.03, 41.93, -69.94, 41.995]
catalog = Client.open("https://earth-search.aws.element84.com/v1")


def level_at(when):
    """Water level at Chatham, m NAVD88, at the given UTC minute."""
    reply = requests.get("https://api.tidesandcurrents.noaa.gov/api/prod/datagetter", timeout=60, params={
        "station": "8447435", "product": "water_level", "datum": "NAVD", "units": "metric",
        "time_zone": "gmt", "format": "json", "begin_date": when.strftime("%Y%m%d %H:%M"), "range": 1}).json()
    return float(reply["data"][0]["v"]) if "data" in reply else np.nan


candidates = []
for year in (2017, 2021, 2023, 2025):
    for item in catalog.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime=f"{year}-05-01/{year}-09-30",
                               query={"eo:cloud_cover": {"lt": 5}}).items():
        when = pd.Timestamp(item.properties["datetime"]).tz_convert(None)
        candidates.append({"year": year, "item": item, "time": when, "level": level_at(when)})
candidates = pd.DataFrame(candidates).dropna(subset=["level"])
TARGET = candidates.loc[candidates["year"] == 2017, "level"].min()
chosen = (candidates.assign(off=(candidates["level"] - TARGET).abs()).sort_values("off")
          .groupby("year").head(1).sort_values("year").reset_index(drop=True))


def read_bands(item, bands):
    out = {}
    for band in bands:
        with rasterio.open(item.assets[band].href) as src:
            win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", src.crs, *BBOX),
                                               transform=src.transform).round_offsets().round_lengths()
            if not out:
                grid = (int(win.height), int(win.width))
                out["transform"], out["crs"] = src.window_transform(win), src.crs
            out[band] = src.read(1, window=win, out_shape=grid, resampling=Resampling.bilinear).astype("float32")
    return out


def waterline(item, threshold=0.0):
    b = read_bands(item, ("green", "swir16"))
    total = b["green"] + b["swir16"]
    mndwi = np.where(total > 0, (b["green"] - b["swir16"]) / np.where(total > 0, total, 1), np.nan)
    water = ndimage.binary_closing(ndimage.binary_opening(mndwi > threshold, np.ones((3, 3))), np.ones((3, 3)))
    labels, _ = ndimage.label(water)
    east = np.unique(labels[:, -5:])                         # closing erodes the image border
    ocean = ndimage.binary_fill_holes(np.isin(labels, east[east > 0]))   # whitecaps leave holes
    sea = gpd.GeoSeries([shape(g) for g, v in shapes(ocean.astype("uint8"), transform=b["transform"]) if v],
                        crs=b["crs"]).union_all()
    frame = box(*rasterio.transform.array_bounds(*ocean.shape, b["transform"])).buffer(-30)
    return gpd.GeoDataFrame(geometry=[sea.boundary.intersection(frame)], crs=b["crs"])


lines = {row.year: waterline(row.item) for row in chosen.itertuples()}
print(chosen[["year", "time", "level"]].round(3).to_string(index=False))
print({year: f"{line.length.sum() / 1000:.1f} km" for year, line in lines.items()})

# a baseline that does not move: a straight fit to the 2025 line, 400 m landward
pts = np.vstack([np.asarray(g.coords) for g in lines[2025].geometry.iloc[0].geoms])
a, b0 = np.polyfit(pts[:, 1], pts[:, 0], 1)
y0, y1 = pts[:, 1].min(), pts[:, 1].max()
baseline_line = LineString([(a * y0 + b0 - 400, y0), (a * y1 + b0 - 400, y1)])
s2017, s2025 = lines[2017], lines[2025]

 year                    time  level
 2017 2017-07-04 15:35:56.461  0.350
 2021 2021-08-02 15:41:26.912 -0.065
 2023 2023-08-12 15:41:31.027  0.433
 2025 2025-08-11 15:41:27.553  0.188
{2017: '11.7 km', 2021: '10.4 km', 2023: '10.5 km', 2025: '10.4 km'}

Example 1 — the full measurement

import numpy as np, geopandas as gpd, pandas as pd
from shapely.geometry import Point

def shoreline_change(baseline, shorelines: dict, spacing=50, length=500):
    """shorelines: {date: GeoDataFrame of lines}, all in the same projected CRS."""
    t = transects(baseline, spacing, length)
    origins = [baseline.interpolate(d) for d in np.arange(0, baseline.length, spacing)]

    rows = []
    for date, gdf in shorelines.items():
        line = gdf.union_all()
        for tid, (tr, o) in enumerate(zip(t.geometry, origins)):
            rows.append({"transect": tid, "date": pd.Timestamp(date),
                         "distance_m": intersect_distance(tr, line, o)})
    df = pd.DataFrame(rows)
    wide = df.pivot(index="transect", columns="date", values="distance_m")

    first, last = wide.columns.min(), wide.columns.max()
    years = (last - first).days / 365.25
    wide["change_m"] = wide[last] - wide[first]
    wide["rate_m_yr"] = wide["change_m"] / years
    return wide

result = shoreline_change(baseline_line, {"2017-07-04": s2017, "2025-08-11": s2025}, length=800)
print(result[["change_m", "rate_m_yr"]].describe().round(2))
print(f"transects with no intersection: {result['change_m'].isna().sum()}")

date   change_m  rate_m_yr
count    156.00     156.00
mean     -23.13      -2.85
std       22.47       2.77
min      -76.72      -9.47
25%      -38.49      -4.75
50%      -21.92      -2.70
75%      -10.52      -1.30
max       21.92       2.70
transects with no intersection: 3

Example 2 — the uncertainty budget

import numpy as np

def position_uncertainty(georef_m, digitise_m, pixel_m, tidal_m):
    """Root sum of squares of the independent contributions."""
    terms = np.array([georef_m, digitise_m, pixel_m / 2, tidal_m])
    return float(np.sqrt((terms ** 2).sum()))

# georeferencing: how far the two scenes are offset over land, by phase correlation
from skimage.registration import phase_cross_correlation

nir = {year: read_bands(chosen.loc[chosen["year"] == year, "item"].iloc[0], ("nir",))["nir"]
       for year in (2017, 2025)}
inland = np.s_[:, : nir[2017].shape[1] // 2]              # the western half is land
shift, _, _ = phase_cross_correlation(nir[2017][inland], nir[2025][inland], upsample_factor=20)
georef = float(np.hypot(*shift) * 10)

# extraction: how far the line moves if the MNDWI threshold moves from 0.0 to 0.2
alt = shoreline_change(baseline_line, {"2025-08-11": waterline(chosen["item"].iloc[-1], 0.2),
                                       "2025-08-12": s2025}, length=800)
digitise = float(alt["change_m"].abs().median())

# tide: the level difference over the beach slope, both measured
DEM = "https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n42w070/USGS_13_n42w070.tif"
from pyproj import Transformer

tr = transects(baseline_line, 50, 800)
with rasterio.open(DEM) as src:
    win = rasterio.windows.from_bounds(*transform_bounds("EPSG:32619", src.crs, *tr.total_bounds),
                                       transform=src.transform).round_offsets().round_lengths()
    dem = src.read(1, window=win, boundless=True, fill_value=-999999).astype("float32")
    dem[dem < -9000] = np.nan
    dem_transform, to_dem = src.window_transform(win), Transformer.from_crs(32619, src.crs, always_xy=True)
slopes = []
for geom in tr.geometry:
    d = np.arange(0, geom.length, 5.0)
    xy = [geom.interpolate(x) for x in d]
    rows, cols = rasterio.transform.rowcol(dem_transform, *to_dem.transform([p.x for p in xy], [p.y for p in xy]))
    z = dem[np.clip(rows, 0, dem.shape[0] - 1), np.clip(cols, 0, dem.shape[1] - 1)]
    above = np.nonzero(z >= 2.0)[0]                       # the +2 m and 0 m NAVD88 contours
    if len(above) and len(below := np.nonzero(z <= 0.0)[0][np.nonzero(z <= 0.0)[0] > above.max()]):
        slopes.append(2.0 / (d[below.min()] - d[above.max()]))
slope = float(np.median(slopes))
levels = chosen.set_index("year")["level"]
tidal = abs(levels[2025] - levels[2017]) / slope

u_2017 = position_uncertainty(georef_m=georef, digitise_m=digitise, pixel_m=10.0, tidal_m=0.0)
u_2025 = position_uncertainty(georef_m=georef, digitise_m=digitise, pixel_m=10.0, tidal_m=tidal)
u_change = float(np.hypot(u_2017, u_2025))
print(f"co-registration offset {georef:.1f} m, threshold sensitivity {digitise:.1f} m, "
      f"beach slope 1 in {1 / slope:.0f} ({len(slopes)} profiles), tide difference "
      f"{abs(levels[2025] - levels[2017]):.2f} m = {tidal:.1f} m of beach")
print(f"2017 ±{u_2017:.1f} m, 2025 ±{u_2025:.1f} m, change ±{u_change:.1f} m")
unmatched = (candidates.loc[candidates["year"] == 2025, "level"] - levels[2017]).abs().max() / slope
print(f"had the 2025 scene been picked for cloud alone, the tide term could have been {unmatched:.1f} m")

detectable = result["change_m"].abs() > u_change
print(f"transects with detectable change: {detectable.sum()} of {len(result)} "
      f"({detectable.mean():.1%})")

co-registration offset 0.7 m, threshold sensitivity 11.0 m, beach slope 1 in 15 (98 profiles), tide difference 0.16 m = 2.4 m of beach
2017 ±12.1 m, 2025 ±12.3 m, change ±17.2 m
had the 2025 scene been picked for cloud alone, the tide term could have been 15.0 m
transects with detectable change: 94 of 159 (59.1%)

Which term dominates depends on what you controlled. Here the scenes were chosen for matching water levels, so the tide contributed 2.4 m, and the largest term was the extraction itself: moving the MNDWI threshold from 0 to 0.2 moved the line 11 m. Picked for cloud cover alone, the 2025 scene could have added a 15 m tide term. On a gently sloping beach an uncorrected tide state of half a metre is tens of metres of horizontal position, which is why it is the term most often left out and least often small.

Example 3 — a rate from three or more dates

import numpy as np, pandas as pd

def linear_rates(wide_dates: pd.DataFrame):
    """wide_dates: index=transect, columns=Timestamps, values=distance in metres."""
    dates = wide_dates.columns
    x = np.array([(d - dates[0]).days / 365.25 for d in dates])
    out = []
    for tid, row in wide_dates.iterrows():
        y = row.values.astype(float)
        ok = np.isfinite(y)
        if ok.sum() < 3:
            out.append({"transect": tid, "rate_m_yr": np.nan, "se": np.nan, "n": int(ok.sum())})
            continue
        slope, intercept = np.polyfit(x[ok], y[ok], 1)
        resid = y[ok] - (slope * x[ok] + intercept)
        se = float(np.sqrt((resid ** 2).sum() / (ok.sum() - 2)
                           / ((x[ok] - x[ok].mean()) ** 2).sum()))
        out.append({"transect": tid, "rate_m_yr": float(slope), "se": se, "n": int(ok.sum())})
    df = pd.DataFrame(out)
    df["significant"] = df["rate_m_yr"].abs() > 1.96 * df["se"]
    return df

four = shoreline_change(baseline_line, {f"{year}-08-01": line for year, line in lines.items()}, length=800)
rates = linear_rates(four[[c for c in four.columns if isinstance(c, pd.Timestamp)]])
print(rates[["rate_m_yr", "se"]].describe().round(2).loc[["count", "mean", "50%"]])
print(f"significant at 95%: {rates['significant'].sum()} of {rates['rate_m_yr'].notna().sum()} transects")

       rate_m_yr      se
count     156.00  156.00
mean       -3.19    2.33
50%        -3.03    2.26
significant at 95%: 50 of 156 transects

Reporting the standard error per transect turns a map of rates into a map of rates you can defend, and the significant column is usually the one a reviewer asks for.

Transects over lengths and areas

Why transects and not lengths or areas

Length is scale-dependent and therefore not comparable between sources. Area change between two shorelines is well defined but gives one number for a whole cell, hiding erosion in one place and accretion in another. Transects give a signed distance at a fixed location, which is what "the beach retreated 12 m" means and is what a time series needs.

Why the tidal correction can dominate

Horizontal position is the vertical offset divided by the beach slope. On a 2% slope — a typical sandy beach — a 0.5 m difference in water level is 25 m of horizontal shoreline position. Two images at different tide states therefore show an apparent change of tens of metres with nothing happening on the ground, which is larger than most real annual rates.

Why the baseline must be identical between comparisons

The measurement is a distance from the baseline, so a baseline that moves adds its own movement to every transect. That includes regenerating the baseline from one of the shorelines, which is a common shortcut and makes the first date's errors appear in every later measurement.

Why to report detectability rather than change

A map of change values with no uncertainty implies that every arrow is real. On outer Cape Cod the combined uncertainty came to ±17.2 m, and 59% of transects moved by more than that — a more honest summary than the median waterline retreat of 21.9 m on its own. It is a waterline, too: swash and the seasonal beach profile move it as well as erosion, which is another reason to state what was measured.

Concave coasts, multiple crossings and nourishment

• Transects cross each other on a concave coast. Space them from a smoothed baseline.
• Multiple intersections need a rule. Nearest or furthest, applied consistently.
• Spits and inlets migrate along the coast, which transects perpendicular to it measure badly.
• Nourished beaches are not natural signals. Record intervention dates.
• A storm between surveys is an event, not a rate.
• Different sources mean different scales. Generalise both, or compare only like with like.
• Keep transect ids stable across every date, or the time series is meaningless.
• Vegetation lines and wet-dry lines are proxies, not datums; state which you digitised.

Internal links

• How to extract a coastline from a raster in Python — producing the lines to compare
• Tidal datums explained: which shoreline is the shoreline — the correction that dominates the budget
• The coastline paradox explained for GIS work — why not to compare lengths
• How to fetch and align tide gauge data in Python — getting the tide state at acquisition
• How to detect change between satellite images in Python — the raster approach
• How to measure distance accurately in Python — the distance underneath
• Sea level rise data explained — the longer-term signal
• Spatial data quality dimensions — building an uncertainty budget

FAQ

How do I measure shoreline change?

Cast transects perpendicular to a fixed baseline, measure where each shoreline crosses each transect, and difference the distances. Report the rate per year.

Can I just compare the two coastline lengths?

No. Length depends on the source scale, so two lengths from different sources are not comparable and neither is the difference.

What dominates the uncertainty?

Usually the tide state at acquisition. On a 2% beach slope, 0.5 m of water level is 25 m of horizontal position.

How do I handle a transect that crosses the shoreline twice?

Pick a rule — nearest crossing to the baseline, or furthest — apply it everywhere, and count how often it mattered.

How much change is real?

Only change larger than the combined positional uncertainty of both dates. Report the fraction of transects above that threshold.

Do I need more than two dates?

For a rate you can defend, yes. Three or more lets you fit a linear rate per transect and report its standard error.