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. 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.
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(
# Allow one metre of error. In fact, the error
# is at most on the order of 1e-8 metres.
lon, lat, PLATE_CARREE)) <= 1.
land_projected = [
# buffer(0) fixes most problems with geometric
# objects, such as self-crossing.
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 the 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))
[(r * math.sin(t), r * math.cos(t))
for t in numpy.linspace(
0, 2 * math.pi, NUM_VERTICES, endpoint=False)])
# The geometric objects in |land_projected|
# only overlap on their boundaries.
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)
# The exception is raised e.g. for (lat, lon) =
# (6, -135) and land from 'ne_110m_land'.
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)
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
# 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:
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:
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:
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():
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:
lat, lon, land, r0=90, r1=0, rstep=0.1):
def CircleIntersectsLand(lat, lon, r):
circle = Circle(lat, lon, r, proj)
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):
if not CircleIntersectsLand(lat, lon, r):
while rstep > 1e-10:
rstep /= 2.
rm = r + rstep
circle = Circle(lat, lon, rm, proj)
if not CircleIntersectsLand(lat, lon, rm):
r = rm
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:
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.
What hemisphere of Earth has the most ocean? Hints: i. use an equal-area projection, ii. Shapely objects have an area property.
What hemisphere has the most land?
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.
How would the results change on an Earth ellipsoid?