---
jupyter:
  description: 'The pole of inaccessibility of Great Britain with planex2d: the exact distance to Natural Earth''s coast
    on a grid, the largest circle inside the coast grown with SciPy, and the Ordnance Survey''s point near Coton in the
    Elms to check against.'
---

# 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.

![Great Britain's coast with the largest circle that fits inside it, touching the Wash, the Severn and the Mersey](../_static/gallery/pole-of-inaccessibility.png){.result}

```python
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
```

```python nbsphinx="hidden"
from docs_figures import FACE, INK, MAP, PLANEX, RED, SOFT, Row
```

## The coast

The coast is the same as on [The coastline paradox](coastline-paradox.md), Great Britain's outline from the 1:10m
coastline of [Natural Earth](https://www.naturalearthdata.com), stored as [britain-coast.csv](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](https://www.naturalearthdata.com/about/terms-of-use/). The licenses of this
documentation do not apply to it. Made with Natural Earth.

```python
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")
```

## 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:

```python
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")
```

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

```python nbsphinx="hide-input"
def mapped(face, plot=True, **_):
    return {"points": [(h[i, j], v[i, j]) for i, j in np.argwhere(inside)], "values": distance[inside],
            "label": "distance to the coast in km"}


Row((island, "The distance to the coast on a grid of 10 km", {"field": mapped, "cmap": MAP.reversed(), "clim": (0, 110),
                                                               "bar": "right"}))
```

## 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:

```python
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")
```

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:

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

## 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:

```python
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")
```

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.

```python nbsphinx="hide-input"
Row((px.Compound([island, px.Circle(pole, r)]), "The largest circle and the second one",
     {"dashed": [px.Circle(second, r2)], "marks": [tuple(p) for p in touching], "dots": [farm]}))
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
fig, ax = plt.subplots(figsize=(5, 7))
island.Plot(ax=ax, face_col=FACE, line_col=INK, legend=False)
ax.add_patch(plt.Circle(pole, r, fill=False, color=INK, lw=2))
ax.add_patch(plt.Circle(second, r2, fill=False, color=SOFT, ls="--"))
ax.plot(*touching.T, "x", color=RED, markersize=9, mew=2)
ax.plot(*farm, "o", color=PLANEX)
ax.axis("off")
save_figure(fig, "pole-of-inaccessibility", bare=True)
```