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.8%, and by 68.6% 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.

145 coastline features, 20,000 points

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. Expect a steady shortfall of a few tenths of a per cent — that is the sphere against the ellipsoid, not a bug. Compute the same check on the sphere as well: it should agree to the centimetre, and if it does not, the method is wrong.

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.

Code examples

The data used on this page

The coastline is Natural Earth's 1:10m coastline, kept for 15° W to 9° E and 45 to 64° N so that every point has its true nearest coast in the set. The points are 20,000 random locations, seed 0, uniform in longitude and latitude over 8° W to 2° E and 50 to 59° N — Great Britain, Ireland's east coast, and the seas around them, so both land and sea.

import os

import geopandas as gpd
import numpy as np
import requests

if not os.path.exists("ne_10m_coastline.zip"):
    reply = requests.get("https://naciscdn.org/naturalearth/10m/physical/ne_10m_coastline.zip", timeout=300)
    reply.raise_for_status()
    with open("ne_10m_coastline.zip", "wb") as out:
        out.write(reply.content)
coast = gpd.read_file("ne_10m_coastline.zip").cx[-15:9, 45:64]

rng = np.random.default_rng(0)
points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(rng.uniform(-8, 2, 20_000),
                                                      rng.uniform(50, 59, 20_000)), crs=4326)
print(f"{len(coast)} coastline features, {len(points):,} points")

145 coastline features, 20,000 points

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.2 km, max 243.4 km, within 10 km 23.1%

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 - great_circle / 1000) / (great_circle / 1000)     # same vertices, same points
print(f"median relative error {np.median(rel):+.1%}, "
      f"90th percentile {np.percentile(rel, 90):+.1%}")

median relative error +31.8%, 90th percentile +68.6%

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")
rel_err = err / exact
print(f"as a share of the distance: median {np.median(rel_err):+.2%}, "
      f"range {rel_err.min():+.2%} to {rel_err.max():+.2%}")

_, _, exact_sphere = Geod(a=R, b=R).inv(points.geometry.x.values[sample], points.geometry.y.values[sample],
                                         nearest[sample, 0], nearest[sample, 1])
print(f"vs the same calculation on the sphere: max |error| {np.abs(dist_m[sample] - exact_sphere).max():.2f} m")

vs exact geodesic: mean -107.27 m, max |error| 733.90 m
as a share of the distance: median -0.25%, range -0.36% to -0.04%
vs the same calculation on the sphere: max |error| 0.00 m

All of the residual here is the sphere against the WGS 84 ellipsoid: the same calculation on a sphere agrees to 0.00 m, and against the ellipsoid the spherical distance runs 0.04–0.36% short, up to 734 m at a couple of hundred kilometres. If that matters, correct it with Geod.inv on the nearest vertex, as here. A residual that does not scale with distance means something else is wrong, usually the densification.

Example 4 — the other figures on this page

print("chord against arc")
for arc_km in (100, 1000):
    arc = arc_km * 1000
    chord = 2 * R * np.sin(arc / (2 * R))
    print(f"  {arc_km:5,} km arc: chord {arc - chord:,.2f} m shorter")

undense = great_circle
print(f"\nundensified vertices against 200 m densified: median {np.median(undense - dist_m):.0f} m, "
      f"95th percentile {np.percentile(undense - dist_m, 95):,.0f} m, max {np.max(undense - dist_m):,.0f} m")

chord against arc
    100 km arc: chord 1.03 m shorter
  1,000 km arc: chord 1,026.22 m shorter

undensified vertices against 200 m densified: median 0 m, 95th percentile 416 m, max 7,706 m

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.

Densification and the ellipsoid are two different errors

The ellipsoidal correction is proportional to distance: a quarter of a per cent at the median here, up to 734 m for the farthest points. The vertex error is not proportional to anything. 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, even for a point a few hundred metres offshore. On this data, skipping densification changed the median distance by nothing but put 5% of points more than 416 m too far away, and one 7.7 km too far. Densifying to 200 m bounds that error at 100 m.

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.

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.

Internal links

• The coastline paradox explained for GIS work — why the source scale changes the answer
• How to measure distance accurately in Python — the general version
• How to find the nearest point in GeoPandas — the simple projected approach
• Nearest neighbour joins explained — when sjoin_nearest is enough
• Maritime boundaries explained: EEZ, territorial sea and baselines — distances from the same line
• Projected versus geographic CRS explained — why degrees are not metres
• How to use a spatial index in GeoPandas — indexing for the projected route
• How to extract a coastline from a raster in Python — producing the line

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.8%.

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.