The pole of inaccessibility

Which point of Great Britain lies farthest from the sea? It is the island’s pole of inaccessibility, the point of a region that lies farthest from its boundary, and the center of the largest circle that fits inside it. Such a circle touches the boundary in at least two places, and in general in three that no half of the circle holds alone, since otherwise it could move away from them and grow. The Ordnance Survey puts Great Britain’s pole at Church Flatts Farm near Coton in the Elms in Derbyshire, 113 km (70 miles) from the mean low water line and equally far from the Wash, the estuary of the Dee and the Severn.

Here, the distance to the coast is measured on a grid over the island, and the largest circle grows from the grid’s best points.

import math

import matplotlib.pyplot as plt
import numpy as np
import planex as px
from scipy.ndimage import maximum_filter
from scipy.optimize import minimize

The coast

The coast is the same as on The coastline paradox, Great Britain’s outline from the 1:10m coastline of Natural Earth, stored as britain-coast.csv next to this page. It is projected to kilometers with an azimuthal equidistant projection around 2° W, 54.5° N, which keeps distances across the island within about 0.1 %. The data is Natural Earth’s and in the public domain under Natural Earth’s terms of use. The licenses of this documentation do not apply to it. Made with Natural Earth.

coast = np.loadtxt("britain-coast.csv", delimiter=",", comments="#")
island = px.Face([px.FromPolyline(coast, closed=True)])
shore = island.outer
print(len(shore.edges), "lines, coast", round(shore.mass), "km")

3705 lines, coast 7320 km

The distance to the coast

shore.Distance(p) is the exact distance from a point to the coast’s lines. On a grid of 10 km over the island, island.Contains keeps the points inside, and their distances map how far each part of the island lies from the sea:

low, high = coast.min(axis=0), coast.max(axis=0)
h, v = np.meshgrid(np.arange(low[0], high[0], 10.0), np.arange(low[1], high[1], 10.0))
inside = island.Contains(np.column_stack([h.ravel(), v.ravel()])).reshape(h.shape)
distance = np.zeros(h.shape)
for i, j in np.argwhere(inside):
    distance[i, j] = shore.Distance((h[i, j], v[i, j]))
print(inside.sum(), "points inside, the farthest", round(distance.max(), 1), "km from the coast")

2180 points inside, the farthest 106.5 km from the coast

Mapped, the distances rise toward the middle of England, where two areas lie more than 100 km from the sea (Figure 1):

Figure 1: The distance to the coast on a grid of 10 km

The largest circle

Where the distance has a maximum, a circle around that point touches the coast in three places at once. The grid comes only within a few kilometers of such a point. Every grid point of more than 80 km that no neighbor beats is therefore a start, and SciPy’s Nelder-Mead method moves the center from there until the distance grows no more. Two starts can lead to the same circle, so the circles are kept by their centers:

def farthest(start):
    found = minimize(lambda p: -shore.Distance(tuple(p)), start, method="Nelder-Mead",
                     options={"xatol": 1e-6, "fatol": 1e-9})
    return tuple(found.x.round(3).tolist()), -found.fun


peaks = np.argwhere((distance == maximum_filter(distance, size=3)) & (distance > 80))
circles = dict(farthest((h[i, j], v[i, j])) for i, j in peaks)
(pole, r), (second, r2) = sorted(circles.items(), key=lambda c: -c[1])[:2]
print("pole", pole, "radius", round(r, 2), "km")
print("second", second, "radius", round(r2, 2), "km,", round(math.dist(pole, second)), "km away")

pole (29.453, -208.194) radius 109.8 km
second (77.85, -270.616) radius 109.19 km, 79 km away

The larger circle is the pole’s. Three points of the coast lie on it, at the heads of the Wash, the Severn and the Mersey near Runcorn:

touching = coast[np.hypot(*(coast - pole).T) < r + 1e-3]
print(touching)

[[ 134.88 -177.5 ]
 [ -27.39 -302.14]
 [ -49.71 -132.1 ]]

Checked against the Ordnance Survey

The farm lies at 52°43.6′ N, 1°37.2′ W, and the same projection as the coast’s puts it into the map:

def project(lon, lat):
    lon, lat, lon0, lat0 = map(math.radians, (lon, lat, -2.0, 54.5))
    c = math.acos(math.sin(lat0) * math.sin(lat) + math.cos(lat0) * math.cos(lat) * math.cos(lon - lon0))
    k = 6371.0088 * c / math.sin(c)
    return (k * math.cos(lat) * math.sin(lon - lon0),
            k * (math.cos(lat0) * math.sin(lat) - math.sin(lat0) * math.cos(lat) * math.cos(lon - lon0)))


farm = project(-(1 + 37.2 / 60), 52 + 43.6 / 60)
print("pole to the farm", round(math.dist(pole, farm), 1), "km")
print("radius", round(r, 1), "km, the Ordnance Survey's 113 km")

pole to the farm 11.7 km
radius 109.8 km, the Ordnance Survey's 113 km

The pole lies 11.7 km from the farm, and its circle is about 3 km smaller than the Ordnance Survey’s distance, because the two measure to different coasts. The Ordnance Survey measured to the mean low water line on its own maps, where the third point lies on the Dee between Neston and Flint. A map at 1:10 million draws the estuaries far more coarsely, and on it the Mersey near Runcorn comes closer to the pole than the Dee.

The second circle shows how much the answer depends on the coast. Centered near Milton Keynes, it touches the Wash, the Thames at Gravesend and the Severn. It is only 0.6 km smaller, 79 km from the pole, so a coast that reached a little farther up the Mersey would make it the pole. Figure 2 shows the pole’s circle with its three points on the coast, the second circle dashed, and the farm as a blue dot.

Figure 2: The largest circle and the second one