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.
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
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 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
from scipy import optimize import earth def LandIntersection(latlon, land): return earth.GreatCircleLandIntersection( latlon, latlon, 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:
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
- 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?