Indiana Jones Plane Flying Gif

Great Circle Arcs and Path¶

Overview¶

Imagine a plane flying from Cario to Hong Kong. To a passenger, the plane appears to travel a straight path, but the plane actually curves around the surface of the Earth held down by the gravity of the planet.

Great circles are circles that circumnavigate the globe.

Distance between Points on a Great Circle Arc
Spherical Distance to Degrees
Determine the Bearing of a Great Circle Arc
Generate a Great Circle Arc with Intermediate Points
Determine the Midpoint of a Great Circle Arc
Generate a Great Circle Path
Determine an Antipodal Point
Compare Great Circle Arc to Rhumb Line (TODO)

Prerequisites¶

Concepts	Importance	Notes
Numpy	Necessary	Used to work with large arrays
Pandas	Necessary	Used to read in and organize data (in particular dataframes)
Intro to Cartopy	Helpful	Will be used for adding maps to plotting
Matplotlib	Helpful	Will be used for plotting

Time to learn: 40 minutes

Imports¶

Import Packages
Setup location dataframe with coordinates

import pandas as pd       # reading in data for location information from text file
import numpy as np        # working with arrays, vectors, cross/dot products, and radians

from pyproj import Geod   # working with the Earth as an ellipsod (WGS-84)
import geopy.distance     # working with the Earth as an ellipsod

# Get all Coordinates for Locations
location_df = pd.read_csv("../location_full_coords.txt")
location_df = location_df.rename(columns=lambda x: x.strip()) # strip excess white space from column names and values
location_df.head()

location_df.index = location_df["name"]

Distance Between Points on a Great Circle Arc¶

Determine Distance Between Points Mathematically via Unit Sphere¶

Distance between point A (latA, lonA) and point B (latB, lonB):

d=acos(sin(latA)*sin(latB)+cos(latA)*cos(latB)*cos(lonA-lonB))

(1)

For shorter distances (with less rounding errors):

d=2*asin(\sqrt{sin(\frac{latA-latB}{2})^2 + cos(latA)*cos(latB)*sin(\frac{lonA-lonB}{2})^2})

(2)

Ed Williams: Distance Between Points

def distance_between_points_default(start_point=None, end_point=None):
    earth_radius = 6378137  # meters
    latA = np.deg2rad(location_df.loc[start_point, "latitude"])
    lonA = np.deg2rad(location_df.loc[start_point, "longitude"])
    latB = np.deg2rad(location_df.loc[end_point, "latitude"])
    lonB = np.deg2rad(location_df.loc[end_point, "longitude"])

    distance_default = np.arccos(np.sin(latA)*np.sin(latB)+np.cos(latA)*np.cos(latB)*np.cos(lonA-lonB))
    return distance_default * earth_radius

def distance_between_points_small(start_point=None, end_point=None):
    earth_radius = 6378137  # meters
    latA = np.deg2rad(location_df.loc[start_point, "latitude"])
    lonA = np.deg2rad(location_df.loc[start_point, "longitude"])
    latB = np.deg2rad(location_df.loc[end_point, "latitude"])
    lonB = np.deg2rad(location_df.loc[end_point, "longitude"])

    distance_small = 2 * np.arcsin(np.sqrt((np.sin((latA-latB)/2))**2 + np.cos(latA)*np.cos(latB)*(np.sin((lonA-lonB)/2))**2))
    return distance_small * earth_radius

Additional types of distance measuerments:

Haversine
Vincenty Sphere Great Circle Distance
Vincenty Ellipsoid Great Circle Distance
Meeus Great Circle Distance

Determine Distance Points via Python Package `pyproj`¶

pyproj accounts for different ellipsoids like WGS-84.

First, setup a ellipsoid to represent the Earth (“WGS-84”):

# Distance between Boulder and Boston
geodesic = Geod(ellps="WGS84")
_, _, distance_meter =  geodesic.inv(location_df.loc["boulder", "longitude"],
                                     location_df.loc["boulder", "latitude"],
                                     location_df.loc["boston", "longitude"],
                                     location_df.loc["boston", "latitude"])

print(f"Distance between coordinates (ellipsoid)   = {distance_meter/1000} km")
distance_unit_sphere_default = distance_between_points_default("boulder", "boston")
print(f"Distance between coordinates (unit sphere) = {distance_unit_sphere_default/1000} km")
distance_unit_sphere_small = distance_between_points_small("boulder", "boston")
print(f"Distance between coordinates (unit sphere) = {distance_unit_sphere_small/1000} km")

Distance between coordinates (ellipsoid)   = 2862.5974799145215 km
Distance between coordinates (unit sphere) = 2858.532213639344 km
Distance between coordinates (unit sphere) = 2858.5322136393447 km

Compared to the distance from the associated airports in Denver and Boston (DIA to Logan) which has a distance of 2823 km, using Denver instead of Boulder.

Spherical Distance to Degrees¶

Convert a distance from meters to degrees, measured along the great circle sphere with a constant spherical radius of ~6371 km (mean radius of Earth).

See also: ObsPy kilometer2degrees()

# assumes a spherical Earth
earth_radius = 6378.137 # km

def km_to_degree_distance(distance_km=None):
    return distance_km / (2 * earth_radius * np.pi / 360)

def degree_to_km_distance(distance_degree=None):
    return distance_degree * (2 * earth_radius * np.pi / 360)

print(f"300 km to degrees = {km_to_degree_distance(300)} degrees")
print(f"6.381307 degree to km = {degree_to_km_distance(6.381307)} km")

300 km to degrees = 2.6949458523585643 degrees
6.381307 degree to km = 710.3638458355522 km

Determine the Bearing of a Great Circle Arc¶

Determine the Bearing Mathematically via Unit Sphere¶

Bearing between point A (latA, lonA) and point B (latB, lonB):

x = cos(latA) * sin(latB) - sin(latA) * cos(latB) * cos(lonB - lonA)

(3)

y = sin(lonB - lonA) * cos(latB)

(4)

θ = atan2(y, x)

(5)

Movable Type: Bearing

def bearing_between_points_unit_sphere(start_point=None, end_point=None):
    latA = np.deg2rad(location_df.loc[start_point, "latitude"])
    lonA = np.deg2rad(location_df.loc[start_point, "longitude"])
    latB = np.deg2rad(location_df.loc[end_point, "latitude"])
    lonB = np.deg2rad(location_df.loc[end_point, "longitude"])

    x = np.cos(latA) * np.sin(latB) - np.sin(latA) * np.cos(latB) * np.cos(lonB - lonA)
    y = np.sin(lonB - lonA) * np.cos(latB)
    bearing = np.arctan2(y, x)
    return np.rad2deg(bearing) % 360

Determine the Bearing via Python Package `pyproj`¶

pyproj accounts for different ellipsoids like WGS-84:

def bearing_between_points_ellps(start_point=None, end_point=None):
    geodesic = Geod(ellps="WGS84")
    fwd_bearing, _, _ =  geodesic.inv(location_df.loc[start_point, "longitude"],
                                        location_df.loc[start_point, "latitude"],
                                        location_df.loc[end_point, "longitude"],
                                        location_df.loc[end_point, "latitude"])
    return fwd_bearing

### Compare Unit Sphere and Ellipsoid
beaing_ellps = bearing_between_points_ellps("boulder", "boston")
print(f"forward bearing between coordinates (ellipsoid)   = {beaing_ellps} Degrees")
bearing_us = bearing_between_points_unit_sphere("boulder", "boston")
print(f"forward bearing between coordinates (unit sphere) = {bearing_us} Degrees")

forward bearing between coordinates (ellipsoid)   = 73.51048829569022 Degrees
forward bearing between coordinates (unit sphere) = 73.49180375272644 Degrees

Generating a Great Circle Arc with Intermediates Points¶

Determine Intermediate Points Mathemetically via Unit Sphere¶

Determine the points (lat, lon) a given fraction of a distance (d) between a starting points A (latA, lonA) and the final point B (latB, lonB) where f is a fraction along the great circle arc. f=0 is point A and f=1 is point B.

Note: The points cannot be antipodal because the path is undefined

Where, antipodal is defined by:

latA + latB = 0

(6)

abs(lonA - lonB) = pi

(7)

Where the distance between two points is the angular distance:

d = \frac{\text{total distance of arc}}{\text{earth's radius}}

(8)

The intermediate points (lat, lon) along a given path starting point to end point:

A = sin(\frac{(1-f) * d}{sin(d)}

(9)

B = \frac{sin(f*d)}{sin(d)}

(10)

x = A * cos(latA) * cos(lonA) + B * cos(latB) * cos(lonB)

(11)

y = A * cos(latA) * sin(lonA) + B * cos(latB) * sin(lonB)

(12)

z = A * sin(latA) + B * sin(latB)

(13)

lat = atan2(z, \sqrt{x^2 + y^2})

(14)

lon = atan2(y, x)

(15)

def intermediate_points(start_point=None, end_point=None,
                        fraction=None, distance=None):
    earth_radius = 6378137  # meters
    total_distance = distance / earth_radius
    latA = np.deg2rad(location_df.loc[start_point, "latitude"])
    lonA = np.deg2rad(location_df.loc[start_point, "longitude"])
    latB = np.deg2rad(location_df.loc[end_point, "latitude"])
    lonB = np.deg2rad(location_df.loc[end_point, "longitude"])

    A = np.sin((1-fraction) * total_distance) / np.sin(total_distance)
    B = np.sin(fraction * total_distance) / np.sin(total_distance)
    x = (A * np.cos(latA) * np.cos(lonA)) + (B * np.cos(latB) * np.cos(lonB))
    y = (A * np.cos(latA) * np.sin(lonA)) + (B * np.cos(latB) * np.sin(lonB))
    z = (A * np.sin(latA)) + (B * np.sin(latB))
    lat = np.arctan2(z, np.sqrt(x**2 + y**2))
    lon = np.arctan2(y, x)
    return (np.rad2deg(lat), np.rad2deg(lon))

def calculate_intermediate_pts(start_point=None, end_point=None,
                               fraction=None, total_distance_meter=None):
    fractions = np.arange(0, 1+fraction, fraction)
    intermediate_lat_lon = []
    for fractional in fractions:
        intermediate_pts = intermediate_points(start_point, end_point,
                                                fractional, total_distance_meter)
        intermediate_lat_lon.append(intermediate_pts)
    return intermediate_lat_lon

Determine Intermediate Points via Python Package `pyproj` and `geopy`¶

Interpolate with N total equally spaced number of points
Interpolate every N meters
Interpolate a fractional distance along arc

def interpolate_points_along_gc(lat_start,
                                lon_start,
                                lat_end,
                                lon_end,
                                distance_between_points_meter): 
    lat_lon_points = [(lat_start, lon_start)]
    
    # move to next point when distance between points is less than the equal distance
    move_to_next_point = True
    while(move_to_next_point):
        forward_bearing, _, distance_meters = geodesic.inv(lon_start,
                                                            lat_start, 
                                                            lon_end,
                                                            lat_end)
        if distance_meters < distance_between_points_meter:
            # ends before overshooting
            move_to_next_point = False
        else:
            start_point = geopy.Point(lat_start, lon_start)
            distance_to_move = geopy.distance.distance(
                            kilometers=distance_between_points_meter /
                            1000)  # distance to move towards the next point
            final_position = distance_to_move.destination(
                            start_point, bearing=forward_bearing)
            lat_lon_points.append((final_position.latitude, final_position.longitude))
            # new starting position is newly found end position
            lon_start, lat_start = final_position.longitude, final_position.latitude
    lat_lon_points.append((lat_end, lon_end))
    return lat_lon_points

Plot Arcs as Points on a World Map¶

import matplotlib.pyplot as plt
from cartopy import crs as ccrs, feature as cfeature

def plot_coordinate(lst_of_coords=None, title=None):
    # Set up world map plot on the United States
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    lon_west, lon_east, lat_south, lat_north = -130, -60, 20, 60
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.STATES, edgecolor="black")
    
    # Plot Latitude/Longitude Location
    longitudes = [x[1] for x in lst_of_coords] # longitude
    latitudes = [x[0] for x in lst_of_coords] # latitude
    plt.plot(longitudes, latitudes)
    plt.scatter(longitudes, latitudes)
    
    # Setup Axis Limits and Title/Labels
    plt.title(title)
    plt.show()

Interpolate with N Total Equally Spaced Points¶

n_total_points = 10 # total points (n points)
distance_between_points_meter = distance_meter / (n_total_points + 1)
print(f"Each point will be separated by {distance_between_points_meter} meters ({distance_between_points_meter/1000} km)")

Each point will be separated by 260236.13453768377 meters (260.23613453768377 km)

lat_start, lon_start = location_df.loc[["boulder"]]["latitude"].iloc[0], location_df.loc[["boulder"]]["longitude"].iloc[0]
lat_end, lon_end = location_df.loc[["boston"]]["latitude"].iloc[0], location_df.loc[["boston"]]["longitude"].iloc[0]

intermediate_geodesic = interpolate_points_along_gc(lat_start,
                                          lon_start,
                                          lat_end,
                                          lon_end,
                                          distance_between_points_meter)
print(f"{len(intermediate_geodesic)} Total Points")
intermediate_geodesic

12 Total Points

[(np.float64(40.015), np.float64(-105.2705)),
 (40.64283438472448, -102.32002071588883),
 (41.19386139956729, -99.31719425393653),
 (41.665293789240074, -96.2672998277903),
 (42.05464865958041, -93.17653047007545),
 (42.35980367525435, -90.05192021556942),
 (42.579048241302566, -86.9012334462751),
 (42.71112689737456, -83.73281874084786),
 (42.75527239726804, -80.5554326250441),
 (42.711226442193585, -77.37804142647055),
 (42.579246749547636, -74.20961159223961),
 (np.float64(42.3601), np.float64(-71.0589))]

plot_coordinate(intermediate_geodesic,
                title=f"Interpolate {n_total_points} Points")

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/formatters.py:402, in BaseFormatter.__call__(self, obj)
    400     pass
    401 else:
--> 402     return printer(obj)
    403 # Finally look for special method names
    404 method = get_real_method(obj, self.print_method)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    167     from matplotlib.backend_bases import FigureCanvasBase
    168     FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
    171 data = bytes_io.getvalue()
    172 if fmt == 'svg':

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2154     # we do this instead of `self.figure.draw_without_rendering`
   2155     # so that we can inject the orientation
   2156     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157         self.figure.draw(renderer)
   2158 if bbox_inches:
   2159     if bbox_inches == "tight":

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     92 @wraps(draw)
     93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94     result = draw(artist, renderer, *args, **kwargs)
     95     if renderer._rasterizing:
     96         renderer.stop_rasterizing()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
   3254             # ValueError can occur when resizing a window.
   3256     self.patch.draw(renderer)
-> 3257     mimage._draw_list_compositing_images(
   3258         renderer, self, artists, self.suppressComposite)
   3260     renderer.close_group('figure')
   3261 finally:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes.draw(self, renderer, **kwargs)
    504         self.imshow(img, extent=extent, origin=origin,
    505                     transform=factory.crs, *factory_args[1:],
    506                     **factory_kwargs)
    507 self._done_img_factory = True
--> 509 return super().draw(renderer=renderer, **kwargs)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
   3223 if artists_rasterized:
   3224     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
   3227     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3229 renderer.close_group('axes')
   3230 self.stale = False

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/feature_artist.py:184, in FeatureArtist.draw(self, renderer)
    179     geoms = self._feature.geometries()
    180 else:
    181     # For efficiency on local maps with high resolution features (e.g
    182     # from Natural Earth), only create paths for geometries that are
    183     # in view.
--> 184     geoms = self._feature.intersecting_geometries(extent)
    186 stylised_paths = {}
    187 # Make an empty placeholder style dictionary for when styler is not
    188 # used.  Freeze it so that we can use it as a dict key.  We will need
    189 # to unfreeze all style dicts with dict(frozen) before passing to mpl.

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:312, in NaturalEarthFeature.intersecting_geometries(self, extent)
    305 """
    306 Returns an iterator of shapely geometries that intersect with
    307 the given extent.
    308 The extent is assumed to be in the CRS of the feature.
    309 If extent is None, the method returns all geometries for this dataset.
    310 """
    311 self.scaler.scale_from_extent(extent)
--> 312 return super().intersecting_geometries(extent)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:112, in Feature.intersecting_geometries(self, extent)
    109 if extent is not None and not np.isnan(extent[0]):
    110     extent_geom = sgeom.box(extent[0], extent[2],
    111                             extent[1], extent[3])
--> 112     return (geom for geom in self.geometries() if
    113             geom is not None and extent_geom.intersects(geom))
    114 else:
    115     return self.geometries()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:297, in NaturalEarthFeature.geometries(self)
    295     if reader.crs is not None:
    296         self._crs = reader.crs
--> 297     geometries = tuple(reader.geometries())
    298     _NATURAL_EARTH_GEOM_CACHE[key] = geometries
    299 else:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/shapereader.py:173, in BasicReader.geometries(self)
    161 def geometries(self):
    162     """
    163     Return an iterator of shapely geometries from the shapefile.
    164 
   (...)    171 
    172     """
--> 173     for shape in self._reader.iterShapes(bbox=self._bbox):
    174         # Skip the shape that can not be represented as geometry.
    175         if shape.shapeType != shapefile.NULL:
    176             yield sgeom.shape(shape)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:2806, in Reader.iterShapes(self, bbox)
   2802 if self.numShapes:
   2803     # Iterate exactly the number of shapes from shx header
   2804     for i in range(self.numShapes):
   2805         # MAYBE: check if more left of file or exit early?
-> 2806         shape = self.__shape(oid=i, bbox=bbox)
   2807         if shape:
   2808             yield shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:2684, in Reader.__shape(self, oid, bbox)
   2681 shapeType = unpack("<i", b_io.read(4))[0]
   2683 ShapeClass = SHAPE_CLASS_FROM_SHAPETYPE[shapeType]
-> 2684 shape = ShapeClass.from_byte_stream(
   2685     shapeType, b_io, recLength_bytes, oid=oid, bbox=bbox
   2686 )
   2688 # Seek to the end of this record as defined by the record header because
   2689 # the shapefile spec doesn't require the actual content to meet the header
   2690 # definition.  Probably allowed for lazy feature deletion.
   2691 # f.seek(next_shape)
   2693 return shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1178, in _CanHaveBBox.from_byte_stream(cls, shapeType, b_io, next_shape, oid, bbox)
   1171     return None
   1173 nParts: Optional[int] = (
   1174     _CanHaveParts._read_nparts_from_byte_stream(b_io)
   1175     if shapeType in _CanHaveParts_shapeTypes
   1176     else None
   1177 )
-> 1178 nPoints: int = cls._read_npoints_from_byte_stream(b_io)
   1179 # Previously, we also set __zmin = __zmax = __mmin = __mmax = None
   1181 if nParts:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1126, in _CanHaveBBox._read_npoints_from_byte_stream(b_io)
   1124 @staticmethod
   1125 def _read_npoints_from_byte_stream(b_io: ReadableBinStream) -> int:
-> 1126     return unpack("<i", b_io.read(4))[0]

error: unpack requires a buffer of 4 bytes

<Figure size 1500x1000 with 2 Axes>

Interpolate every N meters¶

distance_between_points_meter = 112000
print(f"Each point will be separated by {distance_between_points_meter} meters ({distance_between_points_meter/1000} km)")

Each point will be separated by 112000 meters (112.0 km)

lat_start, lon_start = location_df.loc["boulder", "latitude"], location_df.loc["boulder", "longitude"]
lat_end, lon_end = location_df.loc["boston", "latitude"], location_df.loc["boston", "longitude"]

intermediate_geodesic = interpolate_points_along_gc(lat_start,
                                          lon_start,
                                          lat_end,
                                          lon_end,
                                          distance_between_points_meter)
print(f"{len(intermediate_geodesic)} Total Points")
intermediate_geodesic

27 Total Points

[(np.float64(40.015), np.float64(-105.2705)),
 (40.29443152420481, -104.00739372929635),
 (40.55994883031889, -102.73410083358552),
 (40.811312209497714, -101.45097760323463),
 (41.04829010633232, -100.15841313140038),
 (41.27065999040772, -98.8568290209118),
 (41.47820922569627, -97.54667886388195),
 (41.67073593071231, -96.22844748540258),
 (41.848049822017195, -94.90264994540797),
 (42.00997303342117, -93.56983029584387),
 (42.156340903095064, -92.23056009359964),
 (42.287002720790326, -90.88543667319323),
 (42.40182242747838, -89.53508118686946),
 (42.50067925996674, -88.1801364234981),
 (42.583468333429785, -86.82126442135142),
 (42.65010115530707, -85.45914389340653),
 (42.700506064664665, -84.09446748716023),
 (42.73462859187762, -82.72793890397021),
 (42.75243173435998, -81.36026990555973),
 (42.75389614502746, -79.99217723746453),
 (42.73902023120852, -78.62437950079696),
 (42.707820162798164, -77.25759400470386),
 (42.66032978955507, -75.89253363227026),
 (42.596600468550456, -74.52990375235937),
 (42.51670080386281, -73.1703992089883),
 (42.42071630165334, -71.81470141834599),
 (np.float64(42.3601), np.float64(-71.0589))]

plot_coordinate(intermediate_geodesic,
                title=f"Interpolate every {distance_between_points_meter/1000} km")

Interpolate a fractional distance along arc¶

fraction = 1/10

distance_between_points_meter = fraction * distance_meter
print(f"Each point will be separated by {distance_between_points_meter} meters ({distance_between_points_meter/1000} km)")

Each point will be separated by 286259.74799145217 meters (286.2597479914522 km)

lat_start, lon_start = location_df.loc["boulder", "latitude"], location_df.loc["boulder", "longitude"]
lat_end, lon_end = location_df.loc["boston", "latitude"], location_df.loc["boston", "longitude"]

intermediate_ellipsoid = interpolate_points_along_gc(lat_start,
                                          lon_start,
                                          lat_end,
                                          lon_end,
                                          distance_between_points_meter)
print(f"{len(intermediate_ellipsoid)} Total Points")
intermediate_ellipsoid

11 Total Points

[(np.float64(40.015), np.float64(-105.2705)),
 (40.70144152851926, -102.0220139666611),
 (41.29459964930597, -98.7107954391427),
 (41.790822409112124, -95.34405626799516),
 (42.18692017366623, -91.93032741511365),
 (42.4802543416051, -88.47932636547252),
 (42.66881568690329, -85.00175846666659),
 (42.75128706952139, -81.50905885116296),
 (42.72708599453554, -78.01308797522954),
 (42.59638380227174, -74.52579917182065),
 (np.float64(42.3601), np.float64(-71.0589))]

plot_coordinate(intermediate_ellipsoid,
                title=f"(Ellipsoid) Interpolate on Fraction {fraction}")

distance_unit_sphere_default = distance_between_points_default("boulder", "boston")
intermediate_unit_sphere = calculate_intermediate_pts("boulder", "boston",
                                               fraction, distance_unit_sphere_default)
print(f"{len(intermediate_unit_sphere)} Total Points")
intermediate_unit_sphere

11 Total Points

[(np.float64(40.015), np.float64(-105.2705)),
 (np.float64(40.69956840796515), np.float64(-102.0224051615289)),
 (np.float64(41.291305596308824), np.float64(-98.71150732723615)),
 (np.float64(41.786544424077), np.float64(-95.3450112546059)),
 (np.float64(42.1820804843035), np.float64(-91.93144054027105)),
 (np.float64(42.475261312703864), np.float64(-88.48050569778584)),
 (np.float64(42.66406520473303), np.float64(-85.00290620334833)),
 (np.float64(42.74716436193266), np.float64(-81.51007310432796)),
 (np.float64(42.723967826794556), np.float64(-78.01386515980239)),
 (np.float64(42.59464096704978), np.float64(-74.52623685097386)),
 (np.float64(42.3601), np.float64(-71.0589))]

plot_coordinate(intermediate_unit_sphere,
                title=f"(Unit Sphere) Interpolate on Fraction {fraction}")

# compare math and geodesic outputs
for i in range(len(intermediate_ellipsoid)):
    _, _, distance_m = geodesic.inv(intermediate_ellipsoid[i][0], intermediate_ellipsoid[i][1],
                                   intermediate_unit_sphere[i][0], intermediate_unit_sphere[i][1])
    if np.isnan(distance_m): distance_m = 0
    print(f"Distance between ellipsoid/unit sphere defined point {i}: {distance_m} meters")

Distance between ellipsoid/unit sphere defined point 0: 0 meters
Distance between ellipsoid/unit sphere defined point 1: 0 meters
Distance between ellipsoid/unit sphere defined point 2: 0 meters
Distance between ellipsoid/unit sphere defined point 3: 0 meters
Distance between ellipsoid/unit sphere defined point 4: 0 meters
Distance between ellipsoid/unit sphere defined point 5: 132.55154860348821 meters
Distance between ellipsoid/unit sphere defined point 6: 136.26443374199349 meters
Distance between ellipsoid/unit sphere defined point 7: 132.0972186842191 meters
Distance between ellipsoid/unit sphere defined point 8: 112.95654352480192 meters
Distance between ellipsoid/unit sphere defined point 9: 71.29186326482345 meters
Distance between ellipsoid/unit sphere defined point 10: 0.0 meters

Determine the Midpoint of a Great Circle Arc¶

The midpoint of an arc can be determiend as a fractional distance along an arc where f = 0.5.

midpoint = distance_meter / 2
lat_start, lon_start = location_df.loc["boulder", "latitude"], location_df.loc["boulder", "longitude"]
lat_end, lon_end = location_df.loc["boston", "latitude"], location_df.loc["boston", "longitude"]

intermediate_geodesic = interpolate_points_along_gc(lat_start,
                                          lon_start,
                                          lat_end,
                                          lon_end,
                                          midpoint)
print(f"{len(intermediate_geodesic)} Total Points")
print(intermediate_geodesic)
print(f"Midpoint = {intermediate_geodesic[1]}")

3 Total Points
[(np.float64(40.015), np.float64(-105.2705)), (42.48025434160511, -88.4793263654725), (np.float64(42.3601), np.float64(-71.0589))]
Midpoint = (42.48025434160511, -88.4793263654725)

distance_unit_sphere_default = distance_between_points_default("boulder", "boston")
intermediate_unit_sphere = calculate_intermediate_pts("boulder", "boston",
                                               1/2, distance_unit_sphere_default)
print(f"{len(intermediate_unit_sphere)} Total Points")
print(intermediate_unit_sphere)
print(f"Midpoint = {intermediate_unit_sphere[1]}")

3 Total Points
[(np.float64(40.015), np.float64(-105.2705)), (np.float64(42.475261312703864), np.float64(-88.48050569778584)), (np.float64(42.3601), np.float64(-71.0589))]
Midpoint = (np.float64(42.475261312703864), np.float64(-88.48050569778584))

# Compare geodesic and unit sphere
_, _, distance_m = geodesic.inv(intermediate_geodesic[1][0], intermediate_geodesic[1][1],
                                   intermediate_unit_sphere[1][0], intermediate_unit_sphere[1][1])
print(f"Distance between geodesic/unint sphere defined midpoint = {distance_m} meters")

Distance between geodesic/unint sphere defined midpoint = 132.55154860505925 meters

Generate a Great Circle Path¶

Generate points along a great circle path bewtween two points.

# Generate Latitude Coordinates based on Longitude Coordinates
def generate_latitude_along_gc(start_point=None, end_point=None, number_of_lon_pts=360):
    lat1 = np.deg2rad(location_df.loc[start_point, "latitude"])
    lon1 = np.deg2rad(location_df.loc[start_point, "longitude"])
    lat2 = np.deg2rad(location_df.loc[end_point, "latitude"])
    lon2 = np.deg2rad(location_df.loc[end_point, "longitude"])

    # Verify not meridian (longitude passes through the poles)
    if np.sin(lon1 - lon2) == 0:
        print("Invalid inputs: start/end points are meridians")
        # plotting meridians at 0 longitude through all latitudes
        meridian_lat = np.arange(-90, 90, 180/len(longitude_lst)) # split in n number
        meridians = []
        for lat in meridian_lat:
            meridians.append((lat, 0))
        return meridians

    # verify not anitpodal (diametrically opposite, points)
    if lat1 + lat2 == 0 and abs(lon1-lon2) == np.pi:
        print("Invalid inputs: start/end points are antipodal")
        return []

    # note: can be expanded to handle input of np arrays by filter out antipodal/merdiain points

    # generate n total number of longitude points along the great circle
    # https://github.com/rspatial/geosphere/blob/master/R/greatCircle.R#L18C3-L18C7
    gc_lon_lst = []
    for lon in range(1, number_of_lon_pts+1):
        new_lon = (lon  * (360/number_of_lon_pts) - 180)
        gc_lon_lst.append(np.deg2rad(new_lon))

    # Intermediate points on a great circle: https://edwilliams.org/avform147.htm"
    gc_lat_lon = []
    for gc_lon in gc_lon_lst:
        num = np.sin(lat1)*np.cos(lat2)*np.sin(gc_lon-lon2)-np.sin(lat2)*np.cos(lat1)*np.sin(gc_lon-lon1)
        den = np.cos(lat1)*np.cos(lat2)*np.sin(lon1-lon2)
        new_lat = np.arctan(num/den)
        gc_lat_lon.append((np.rad2deg(new_lat), np.rad2deg(gc_lon)))
    return gc_lat_lon

def arc_points(start_lat=None,
               start_lon=None,
               end_lat=None,
               end_lon=None,
               n_total_points=10):
    _, _, distance_meter =  geodesic.inv(start_lon,
                                        start_lat,
                                        end_lon,
                                        end_lat)
        
    distance_between_points_meter = distance_meter / (n_total_points + 1)

    
    points_along_arc = interpolate_points_along_gc(start_lat,
                                              start_lon,
                                              end_lat,
                                              end_lon,
                                              distance_between_points_meter)
    return points_along_arc

def plot_coordinate(lat_lon_lst=None,
                    start_point=None, end_point=None,
                    title=None):
    # Set up world map plot
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    lon_west, lon_east, lat_south, lat_north = -180, 180, -90, 90
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.BORDERS, edgecolor='grey')
    ax.add_feature(cfeature.STATES, edgecolor="grey")
        
    # Plot Latitude/Longitude Location of great circle path
    longitudes = [x[1] for x in lat_lon_lst] # longitude
    latitudes = [x[0] for x in lat_lon_lst] # latitude
    plt.plot(longitudes, latitudes, c="cornflowerblue")
    plt.scatter(longitudes, latitudes, c="cornflowerblue")

    # Overly great circle with arc rom start/end point
    start_end_lat_lon = arc_points(location_df.loc[start_point, "latitude"],
                                   location_df.loc[start_point, "longitude"],
                                   location_df.loc[end_point, "latitude"],
                                   location_df.loc[end_point, "longitude"],
                                   n_total_points=20)
    longitudes = [x[1] for x in start_end_lat_lon] # longitude
    latitudes = [x[0] for x in start_end_lat_lon] # latitude
    plt.plot(longitudes, latitudes, c="red")
    plt.scatter(longitudes, latitudes, c="red")
    
    # Setup Axis Limits and Title/Labels
    plt.title(title)
    plt.show()

start_pt = "boulder"
end_pt = "boston"
n_pts = 360
lat_lon_pts = generate_latitude_along_gc(start_pt, end_pt, number_of_lon_pts=n_pts)
plot_coordinate(lat_lon_pts, start_pt, end_pt,
                f"Plot Great Circle, made from the arc {start_pt.title()} to {end_pt.title()}")

start_pt = "arecibo"
end_pt = "greenwich"
n_pts = 360
lat_lon_pts = generate_latitude_along_gc(start_pt, end_pt, number_of_lon_pts=n_pts)
plot_coordinate(lat_lon_pts, start_pt, end_pt,
                f"Plot Great Circle, made from the arc {start_pt.title()} to {end_pt.title()}")

start_pt = "zambezi"
end_pt = "svalbard"
n_pts = 360
lat_lon_pts = generate_latitude_along_gc(start_pt, end_pt, number_of_lon_pts=n_pts)
plot_coordinate(lat_lon_pts, start_pt, end_pt,
                f"Plot Great Circle, made from the arc {start_pt.title()} to {end_pt.title()}")

Determine an Antipodal Point¶

Antipodal is the point on the globe that is on the exact opposite side of the Earth.

Antipodal latitude is defined as:

\text{antipodal latitude} = -1 * \text{latitude}

(16)

Antipodal longitude is defined as:

\text{anitpodal longitude} = (\text{longitude} + 180) \text{ if longitude} \le 0

(17)

\text{anitpodal longitude} = (\text{longitude} - 180) \text{ if longitude} \gt 0

(18)

Antipodes Map

def antipodal(start_point=None):
    anti_lat = -1 * location_df.loc[start_point, "latitude"]
    ref_lon = location_df.loc[start_point, "longitude"]
    if ref_lon > 0:
        anti_lon = ref_lon - 180
    else:
        anti_lon = ref_lon + 180
    #if anti_lon >= 180:
    #    anti_lon = -1 * (anti_lon % 180)
    return (anti_lat, anti_lon)

def is_antipodal(start_point=None, end_point=None):
    lon1 = np.deg2rad(location_df.loc[start_point, "longitude"])
    lat1 = np.deg2rad(location_df.loc[start_point, "latitude"])
    lon2 = np.deg2rad(location_df.loc[end_point, "longitude"])
    lat2 = np.deg2rad(location_df.loc[end_point, "latitude"])
    return lat1 + lat2 == 0 and abs(lon1-lon2) == np.pi

def plot_antipodal(start_point=None):
    # Set up world map plot
    fig = plt.subplots(figsize=(15, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    lon_west, lon_east, lat_south, lat_north = -180, 180, -90, 90
    ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
    ax.coastlines(color="black")
    ax.add_feature(cfeature.BORDERS, edgecolor='grey')
    ax.add_feature(cfeature.STATES, edgecolor="grey")
        
    # Plot Start point
    plt.scatter(location_df.loc[start_point, "longitude"],
                location_df.loc[start_point, "latitude"],
                s=100, c="cornflowerblue", label=start_point.title())

    # Plot Antipodal Point
    antipodal_point = antipodal(start_point)
    plt.scatter(antipodal_point[1], antipodal_point[0], s=100, c="red", label="Antipodal")
    
    # Setup Axis Limits and Title/Labels
    plt.title(f"{start_point.title()} and Antipodal Point {antipodal_point}")
    plt.legend(loc="lower right")
    plt.show()

print(antipodal("boulder"))
plot_antipodal("boulder")

(np.float64(-40.015), np.float64(74.7295))

print(antipodal("svalbard"))
plot_antipodal("svalbard")

(np.float64(-77.875), np.float64(-159.0248))

print(antipodal("cairo"))
plot_antipodal("cairo")

(np.float64(-30.0444), np.float64(-148.7643))

Summary¶

Calculating and mapping the midpoint and intermediate points along the great circle arc and closed circle path.

What’s next?¶

With a great circle arc defined, determine if a third point is along the arc or at what distance it sits from the great circle arc and path.

Resources and references¶

Foundations

Coordinate Types

Tutorials

Great Circles and a Point

Great Circle Arcs and Path