On wet and dry great and small circles

Can you slice Earth flat through its centre while touching no land? If not, then which cut minimizes the length of sliced land? Which cut maximizes it? My favourite search engine returned only partial and inaccessible answers. Which planar cut of Earth that intersects no land is the largest? I have not even found this question online.

With right tools and some brute force, such problems prove easy. I used Cartopy for map projections and Shapely for geometric objects.

Assuming that Earth is a sphere, its intersection with the cutting plane is a great circle when the plane passes through the centre of Earth or a small circle otherwise.
Let us map the surface of Earth with the stereographic projection with centre at the point (lat, lon).

```# The source projection for proj.transform_point()
# and proj.project_geometry(). This one serves as
# a no-op projection.
PLATE_CARREE = crs.PlateCarree()

# Small circles on a sphere are easy. I do not know
# how to specify their equivalent on an ellipsoid.
SPHERICAL_GLOBE = crs.Globe(ellipse='sphere')

def Project(lat, lon, land, klass=crs.Stereographic):
proj = klass(
central_longitude=lon,
central_latitude=lat,
globe=SPHERICAL_GLOBE)
# Allow one metre of error. In fact, the error
# is at most on the order of 1e-8 metres.
assert math.hypot(*proj.transform_point(
lon, lat, PLATE_CARREE)) <= 1.
land_projected = [
# buffer(0) fixes most problems with geometric
# objects, such as self-crossing.
proj.project_geometry(g, PLATE_CARREE).buffer(0)
for g in land]
return proj, land_projected
```

The stereographic projection with centre at (50°N, 20°E)

It tempted me to use the azimuthal equidistant projection known from the United Nations emblem instead. It preserves the ratios of distances along the diameters, and the diameters correspond to great circles. So I could, in theory, get the answer by measuring the length of the intersection of straight lines with land. Unfortunately, in my tests, land on antipodes often disappeared after this projection, which made the intersection undermeasured.

The stereographic projection proved more robust, maybe because it cannot map the antipodes so their neighbourhood is cut from the map. One of circles concentric to its conventional boundary corresponds to a great circle, and the ratios of distances are preserved along it. I measured the length of its intersection with land. The (lat, lon) pairs that determine the centre of the projection were vertices of a more or less evenly spaced mesh covering the sphere.

The code that realizes the above idea is so short that I show it all:

```NUM_VERTICES = 3600

# (lon, clat) is a vertex of a polygon that
# approximates the circle.
# The correction makes the perimeter of the polygon
# equal to the circumference of the circle.
# For a 3600-gon, the correction factor is approximately
# 1 + 7.6e-7.
# The nonlinearity of the projection is negligible.
shift = radius_deg * math.pi / (
NUM_VERTICES * math.sin(math.pi / NUM_VERTICES))
clat = lat + shift
if clat > 90:
clat = 180 - clat
lon = lon + 180
r = math.hypot(*proj.transform_point(
lon, clat, PLATE_CARREE))
return geometry.LinearRing(
[(r * math.sin(t), r * math.cos(t))
for t in numpy.linspace(
0, 2 * math.pi, NUM_VERTICES, endpoint=False)])

def LengthOfIntersectionInKm(
great_circle, land_projected):
# The geometric objects in |land_projected|
# only overlap on their boundaries.
return sum(
g.intersection(great_circle).length
for g in land_projected) / great_circle.length * (

def GreatCircleLandIntersection(lat, lon, land):
proj, land_projected = Project(lat, lon, land)
great_circle = Circle(lat, lon, 90, proj)
try:
return LengthOfIntersectionInKm(
great_circle, land_projected)
except geos.TopologicalError:
# The exception is raised e.g. for (lat, lon) =
# (6, -135) and land from 'ne_110m_land'.
return -1.

def EvaluateForMeshPoints(
fn, lat0, lat1, lon0, lon1, lat_step, land, log):
for lat in numpy.linspace(
lat0, lat1, abs(lat1 - lat0) / lat_step + 1):
# Make the spacing of longitudes approximately
# the same as the spacing of latitudes.
endpoint = (abs(lon1 - lon0) != 360)
num_lons = (
abs(lon1 - lon0) / lat_step *
if abs(lat) != 90 else 1)
for lon in numpy.linspace(
lon0, lon1, num_lons, endpoint=endpoint):
print >>log, lat, lon, fn(lat, lon, land)

# 'ne_110m_land' contains 127 MultiPolygons.
# 'ne_10m_land' contains 1 MultiPolygon with
# 4063 Polygons.
# Intersecting with multiple Polygons proved faster
# than intersecting with one MultiPolygon.
return shapes if len(shapes) > 1 else shapes[0]

def main():
# For this crude approximation, use land data dubbed
# 'scale 1:110 million', with simplified coastline
# and no smaller islands. Read a local copy.
with open('great_circle_110m.txt', 'w') as log:
EvaluateForMeshPoints(
GreatCircleLandIntersection,
0, 90, -180, +180,
1, land, log)
```

On my machine, it was running for about an hour. Let us check out the circles with the shortest land intersection:

```\$ sort -k3 -n great_circle_110m.txt | head
4.0 -108.802228412 2637.84025153
6.0 -116.648044693 2675.44337382
6.0 -114.636871508 2706.60790556
4.0 -109.805013928 2747.17793236
6.0 -115.642458101 2784.60240991
5.0 -114.636871508 2789.91418816
5.0 -113.631284916 2839.06706021
6.0 -112.625698324 2868.19346786
6.0 -113.631284916 2880.62812684
4.0 -110.807799443 2939.329449
```

and those with the longest land intersection:

```\$ sort -k3 -n -r great_circle_110m.txt | head
4.0 -154.930362117 22662.3962712
3.0 -154.930362117 22560.3219242
3.0 -155.933147632 22438.262132
4.0 -153.927576602 22430.6255353
3.0 -153.927576602 22401.2028276
4.0 -155.933147632 22367.8729212
5.0 -153.854748603 22325.8557016
5.0 -152.849162011 22318.0202642
6.0 -153.854748603 22241.6673472
5.0 -154.860335196 22239.5669657
```

For ‘ne_110m_land’ data coming from Natural Earth, the error in the measured length may be on the order of a few kilometres, at least on the equator where we have duplicate data. For instance, the points (0.0 −160.0) and (0.0 20.0) determine the same great circle seen from two sides:

```\$ grep ^0 great_circle_110m.txt | sort -k3 -n -r | head
0.0 -160.0 22153.8112556
0.0 20.0 22153.1791523
0.0 21.0 21444.598626
0.0 -159.0 21441.5582888
0.0 -161.0 21171.2386781
0.0 19.0 21170.0842492
0.0 -156.0 20870.1373834
0.0 24.0 20868.6096133
0.0 22.0 20795.1111135
0.0 -158.0 20794.7358011
```

When we use the most accurate land shapes, ‘ne_10m_land’, the error decreases to a few dozen metres:

```\$ sort -k3 -n -r great_circle_through_poles_10m.txt | head
0.0 -160.0 21984.2636822
0.0 20.0 21984.1709172
0.0 21.0 21459.8373531
0.0 -159.0 21459.7796093
0.0 19.0 21080.7771226
0.0 -161.0 21080.7567727
0.0 -158.0 20878.5033907
0.0 22.0 20878.4813585
0.0 26.0 20706.5668797
0.0 -154.0 20706.5664931
```

Another program searched more thoroughly for great circles using projections with centres in six regions that had the lowest local minimums of the overlap fraction: (4–8°N 108–118°W), (23–25°N 78–81°E), (51–64°N 45–50°E), (49–53°N 3–6°W), (42–45°N 153–157°E), and (76–79°N 98–103°E), and in the region with the highest local maximum: (0–7°N 151–161°W). The latitude step was 30′ and the ’ne_10m_land’ shapes were used. Then I ran another search with latitude step 6′ in the most promising subregions. Finally, I applied `scipy.optimize.minimize()`:

```from scipy import optimize
import earth

def LandIntersection(latlon, land):
return earth.GreatCircleLandIntersection(
latlon[0], latlon[1], land)

def NegativeLandIntersection(latlon, lands):
return -LandIntersection(latlon, lands)

def main():
for x0 in [
(5.2, -114.45), (24.1, 79.17), (57.5, 46.75),
(50.7, -4.76), (44.1, 154.79), (77.4, 100.5)]:
print optimize.minimize(
LandIntersection, x0=x0,
args=(land,), method=METHOD)
print optimize.minimize(
NegativeLandIntersection, x0=(3.4, -154.14),
args=(land,), method=METHOD)
```

Here is the gallery of the locally optimal great circles, in the azimuthal equidistant projection. I swapped the centres of the first and the last map with their antipodes to keep the Old World near the centre.

Onwards to the question about the largest small circle that does not pass through land. For crude results, it is enough to use ‘ne_110m_land’ and evaluate this function at the vertices of the mesh:

```def LargestCircleThroughOcean(
lat, lon, land, r0=90, r1=0, rstep=0.1):

def CircleIntersectsLand(lat, lon, r):
circle = Circle(lat, lon, r, proj)
return any(
g.intersects(circle) for g in land_projected)

proj, land_projected = Project(lat, lon, land)
for r in numpy.linspace(
r0, r1, abs(r1 - r0) / rstep + 1):
try:
if not CircleIntersectsLand(lat, lon, r):
break
except geos.TopologicalError:
return -1.
else:
return r
while rstep > 1e-10:
rstep /= 2.
rm = r + rstep
circle = Circle(lat, lon, rm, proj)
if not CircleIntersectsLand(lat, lon, rm):
r = rm
return r
```

I fine-tuned the results similarly to the way described above. You may have guessed that there is a large circle in the Pacific Ocean between Queensland in Australia and Baja California Sur in Mexico, shown here in the orthographic projection:

But I bet you would not imagine that its radius is beaten by this circle around Australia and Antarctica:

This is its tour de force south of the Rote Island, between the Thursday Island and the Moa Island in the Torres Strait, and through the Louisiade Archipelago:

Here is the link to the code and some open problems: one hard, one tricky, and two easy ones. Please leave a comment if you solve any of them.

1. What hemisphere of Earth has the most ocean? Hints: i. use an equal-area projection, ii. Shapely objects have an `area` property.
2. What hemisphere has the most land?
3. What is the largest small circle that intersects only land? You can confine the search to Asia and find the radii of candidate circles through plain binary search because the largest circle does not contain the Caspian Sea.
4. How would the results change on an Earth ellipsoid?