โš™๏ธ How-To Coastal & marine Performance & scale Coastal & Marine ยท step 11 of 20 Updated

How to compute distance to the coast for many points

Problem statement

Distance to the coast is a covariate in half of coastal analysis โ€” exposure, habitat, property value, salinity, climate โ€” and it is computed wrongly more often than almost anything else, because the obvious method gives plausible numbers that are systematically too large.

The obvious method is to compute distances in degrees and multiply by 111.32 km. Tested on 20,000 random points around Great Britain against a proper spherical calculation, that overestimates by a median of 31.1%, and by 68.7% at the 90th percentile. Nothing about the result looks wrong: the pattern is right, the units are kilometres, and every number is too big.

This guide computes it correctly, fast, at the scale where a spatial join stops being viable.

Quick answer

Project the coastline vertices to 3D Cartesian coordinates and use a k-d tree:

import numpy as np, geopandas as gpd
from scipy.spatial import cKDTree

def to_ecef(lon, lat, R=6_371_008.8):
    la, lo = np.radians(lat), np.radians(lon)
    return np.c_[R * np.cos(la) * np.cos(lo),
                 R * np.cos(la) * np.sin(lo),
                 R * np.sin(la)]

verts = np.vstack([np.asarray(g.coords)
                   for g in coast.geometry.explode(index_parts=False)
                   if g.geom_type == "LineString"])
tree = cKDTree(to_ecef(verts[:, 0], verts[:, 1]))

chord, _ = tree.query(to_ecef(points.geometry.x.values, points.geometry.y.values), k=1)
great_circle = 2 * 6_371_008.8 * np.arcsin(chord / (2 * 6_371_008.8))

The k-d tree works in Euclidean space, so the query returns the straight-line chord through the Earth. Converting the chord to a great-circle distance is one line and matters above about a hundred kilometres โ€” the chord is 1.03 m short at 100 km and 1,026 m short at 1,000 km.

Grid comparing degree distance, projected distance, ECEF k-d tree and exact geodesic across accuracy and speed.
Four methods; only the first is both fast and wrong.

Step-by-step solution

1. Decide what "the coast" is

A coastline layer at 1:10 m generalises heavily, so a point in a fjord or an estuary may be several kilometres from the generalised line and metres from the actual water. Choose the source scale to match the distances you care about.

2. Decide whether you want distance to the line or to the water

They differ for a point on an island, in a bay, or up an estuary. Distance to the nearest coastline vertex is not the same as distance to the nearest sea polygon, and inland water bodies are in neither unless you add them.

3. Pick the method from the extent and the point count

  • Projected CRS with sjoin_nearest โ€” correct and simple for a single national extent with a suitable projection.
  • ECEF k-d tree โ€” correct anywhere, fast for millions of points, and the right default.
  • Exact geodesic per pair โ€” the reference, and far too slow for a real point set.
  • Degrees ร— 111.32 โ€” never.

4. Densify the coastline before building the tree

A k-d tree on vertices finds the nearest vertex, not the nearest point on the line. Where vertices are far apart โ€” a generalised coastline, a straight stretch โ€” the distance is overestimated. Densifying the line to a spacing well below your tolerance fixes it.

5. Query, then convert the chord

cKDTree returns Euclidean distance in the coordinate space it was built in. In ECEF that is the chord; converting to a great-circle arc is 2Rยทasin(chord/2R).

6. Validate against an exact geodesic on a sample

Take a few hundred points, compute the exact distance with pyproj.Geod, and compare. Agreement to a few metres means the method and the densification are right.

7. Sign it if you need to distinguish land from sea

Distance to the coast is usually wanted as positive on land and negative at sea, or vice versa. A point-in-polygon test against the land layer supplies the sign.

Bars showing the relative overestimate from degree-based distance at the median and the ninetieth percentile.
Measured on 20,000 points around Great Britain: the naive method is a third too large at the median.

Code examples

Example 1 โ€” the full pipeline with densification 

import numpy as np, geopandas as gpd
from scipy.spatial import cKDTree
from shapely.geometry import LineString

R = 6_371_008.8

def densify(line, max_seg_m):
    """Insert vertices so no segment exceeds max_seg_m (approximate, in degrees)."""
    step = max_seg_m / 111_320.0
    coords = []
    c = np.asarray(line.coords)
    for a, b in zip(c[:-1], c[1:]):
        d = np.hypot(*(b - a))
        n = max(int(d / step), 1)
        for i in range(n):
            coords.append(a + (b - a) * i / n)
    coords.append(c[-1])
    return LineString(coords)

def coast_distance(points_gdf, coast_gdf, max_seg_m=200):
    lines = [g for g in coast_gdf.geometry.explode(index_parts=False)
             if g.geom_type == "LineString"]
    dense = [densify(g, max_seg_m) for g in lines]
    verts = np.vstack([np.asarray(g.coords) for g in dense])

    tree = cKDTree(to_ecef(verts[:, 0], verts[:, 1]))
    p = to_ecef(points_gdf.geometry.x.values, points_gdf.geometry.y.values)
    chord, idx = tree.query(p, k=1, workers=-1)
    arc = 2 * R * np.arcsin(np.clip(chord / (2 * R), 0, 1))
    return arc, verts[idx]

dist_m, nearest = coast_distance(points, coast)
print(f"{len(points):,} points: median {np.median(dist_m)/1000:.1f} km, "
      f"max {dist_m.max()/1000:.1f} km, within 10 km {np.mean(dist_m < 10_000):.1%}")
20,000 points: median 29.1 km, max 263.3 km, within 10 km 22.5%

workers=-1 uses every core for the query, which is where the time goes on large point sets.

Example 2 โ€” what the naive method costs 

import numpy as np
from scipy.spatial import cKDTree

deg_tree = cKDTree(verts)                       # degrees, treated as a plane
deg_d, _ = deg_tree.query(np.c_[points.geometry.x, points.geometry.y], k=1)
naive_km = deg_d * 111.32

rel = (naive_km - dist_m / 1000) / (dist_m / 1000)
print(f"median relative error {np.median(rel):+.1%}, "
      f"90th percentile {np.percentile(rel, 90):+.1%}")
median relative error +31.1%, 90th percentile +68.7%

The error is always positive and grows with latitude and with how much of the offset is eastโ€“west, so it is not a constant bias that can be calibrated away.

Example 3 โ€” validate against an exact geodesic 

import numpy as np
from pyproj import Geod

geod = Geod(ellps="WGS84")
sample = np.random.default_rng(0).choice(len(points), size=500, replace=False)

_, _, exact = geod.inv(points.geometry.x.values[sample],
                       points.geometry.y.values[sample],
                       nearest[sample, 0], nearest[sample, 1])
err = dist_m[sample] - exact
print(f"vs exact geodesic: mean {err.mean():+.2f} m, "
      f"max |error| {np.abs(err).max():.2f} m")

Residuals of a few metres are the spherical-versus-ellipsoidal difference and the densification spacing. Residuals of hundreds of metres mean the coastline was not densified enough.

Explanation

Why degrees are so badly wrong

A degree of latitude is about 111 km everywhere; a degree of longitude is 111 km at the equator and 111ยทcos(ฯ†) elsewhere โ€” 65 km at 54ยฐN. Treating degree space as a plane therefore stretches every eastโ€“west distance by 1/cos(ฯ†), which at UK latitudes is a factor of 1.70. The measured median error of +31% is that factor diluted by the northโ€“south component of each offset.

Why ECEF plus a k-d tree is the right default

Converting to Earth-centred Cartesian coordinates makes Euclidean distance meaningful โ€” it is the chord โ€” and a k-d tree in three dimensions is fast, exact for nearest-neighbour queries, and has no projection to choose or zone to worry about. It works at the poles and across the antimeridian, which a projected approach does not without care.

Why densification matters more than the sphere-versus-ellipsoid question

The ellipsoidal correction is a few metres in a few hundred kilometres. A generalised coastline with vertices 5 km apart can put the nearest vertex 2.5 km further away than the nearest point on the line. Densifying to 200 m bounds that error at 100 m, which is far below any other uncertainty in the problem.

Why the nearest vertex is not always the nearest coast

A point just inside a narrow estuary is metres from water and tens of kilometres from the open coast along the line. Whether that counts as "close to the coast" is a question about the analysis, not the geometry โ€” and it is why the choice between a coastline layer, a sea polygon and a water-bodies layer has to be made deliberately.

Two scenes contrasting a k-d tree on sparse coastline vertices with the same line densified, showing the nearest-vertex error.
A generalised line with 5 km vertex spacing can overstate the distance by 2.5 km.

Edge cases or notes

  • explode first. MultiLineStrings have no .coords.
  • Use workers=-1 on the query for large point sets.
  • Islands count. A coastline layer includes them; a mainland-only line does not.
  • Estuaries and lakes are in some layers and not others.
  • The antimeridian is harmless in ECEF and a problem in every projected method.
  • Sign the distance with a point-in-polygon test if land and sea must differ.
  • Chord versus arc is 1 m at 100 km and 1 km at 1,000 km.
  • Record the coastline source and scale with the column; the number depends on it.

FAQ

How do I compute distance to the coast in Python?

Convert the coastline vertices and the points to Earth-centred Cartesian coordinates, build a k-d tree, query the nearest vertex, and convert the chord to a great-circle arc.

Why not just use degrees times 111 km?

Because a degree of longitude shrinks with latitude. Measured on 20,000 points around Great Britain, the naive method overestimated by a median of 31.1%.

Do I need to densify the coastline?

Yes. A k-d tree finds the nearest vertex, not the nearest point on the line, so a generalised coastline overestimates by up to half the vertex spacing.

Is a projected CRS good enough?

For a single national extent with a suitable projection, yes. Across continents or near the poles, use the ECEF approach.

How do I get a negative distance at sea?

Sign the result with a point-in-polygon test against the land layer.

Which coastline layer should I use?

One whose generalisation matches your tolerance. A 1:10 m line can be kilometres from the water inside an estuary.