Farie së ndeerme…

Nicola Chetta z albańskich okolic,
którego z chwałą w Contessie zrodzono,
by jak ptak nagi w gnieździe móc się szkolić,
wstąpił w Palermo w Albańczyków grono.

Ich dom się rzeźwić w cieniu mu pozwolił
jak zeschła gałąź pod wina osłoną,
ubrał i taktem, mądrością okolił;
dziś, ksiądz, w Kościele żyje jak mąż z żoną.

Rozpostarł skrzydła, zbłąkana ptaszyna,
w równej w Palermo i w Contessie mierze;
cześć Albańczyków w pismach przypominał.

Niczym jedwabnik zaprzątnięty szczerze
ten skarb splótł, wyszył, spisał — to danina;
szczęściu Albanii składa ją w ofierze.

Coat of arms of Albania
Coat of arms of Albania. Source: Wikipedia

This is the Polish translation of the first sonnet in Albanian, written in 1777 by Kolë Ketta (Nicola Chetta), an Arbëresh (Italian Albanian) from Contessa Entellina in Sicily. I translated and published it to celebrate November 28, the Independence Day of Albania. Here is the original:

Farie së ndeerme në Kuntisë u bii
Kolë Ketta, vllastaar t’Arbrit dhee:
Shkoi në Palermë praa tek e Arbrit shpii,
Ç’e reshti, si zogu rep në folee.

E veshi e e ngjeshi me zakon, me urtësii,
Për në vapët e përtriijti ndënë hjee,
Si të veshkurin rremp stolis një dhrii!
E nani prift klisha kurorë e vee.

Si zok e sbjerrë praa t’du krahët çoi
Në Palermë e n’Kuntisë, po ktei e atei:
Ndeer n’e Arbreshet të gjith gramët kërkoi.

Si krymp mëndafshi gjith e svis vetëhei
E ktë vistaar tuar, kjëndisi e shkroi,
Se të kjosëj gjithë Arbrin ndjeer përtei.

If you speak Albanian and would like to judge my version, here is its literal re-translation into English:

Nicola Chetta from the Albanian region
who was born with glory in Contessa,
to be able to educate himself like a naked bird in a nest,
entered in Palermo the company of Albanians.

Their house let him refresh in shadow
like a dry branch under grapevine’s cover,
dressed [him] and surrounded with tact, wisdom;
today, a priest, he lives in the Church like a husband with [his] wife.

He spread [his] wings, a stray bird,
in Palermo and in Contessa in equal degree;
he reminded of the honour of Albanians in [his] writings.

Like a truly busy silkworm
he plaited, embroidered, wrote down this treasure — it is a tribute;
he offers it as a sacrifice to the fortune of Albania.

Farie së ndeerme…

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.

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