Slicing Earth, Carefully

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')

EARTH_RADIUS_KM = 6370.997

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

Centre in Krakow

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

def Circle(lat, lon, radius_deg, proj):
  # (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 * (
      2 * math.pi * EARTH_RADIUS_KM)

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 *
        math.cos(math.radians(lat)) + endpoint
        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)

def ReadShapes(filename):
  shapes = list(shapereader.Reader(filename).geometries())
  # '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.
  land = ReadShapes('shapes/ne_110m_land')
  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)

METHOD = 'Nelder-Mead'

def main():
  land = earth.ReadShapes('shapes/ne_10m_land')
  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:

 Pacific

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

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:

Torres-Strait

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?
Advertisements
Slicing Earth, Carefully

2 thoughts on “Slicing Earth, Carefully

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s