The Earth Polychromatic Imaging Camera, or EPIC, captured this image of smoke from a wildfire in North America on Aug.15, 2018. Credit: NASA Goddard/ Katy Mersmann

Great Circles and a Point¶

Overview¶

A plane traveling across the country suddenly discovers it is low on fuel! It can no longer make it to the planned airport, instead it has to find the closest airport to its current position that it can make it with its remaining fuel.

Determine the distance of a point to a great circle arc (cross-track and along-track distance)
Determine if a point lies on a great circle arc and path (with and without tolerances)
Determine the distance of a point to a great circle path (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: 30 Minutes

Imports¶

Import Packages
Setup location dataframe with coordinates

import pandas as pd                                    # read in data text file
import numpy as np                                     # working with degrees and radians

from pyproj import Geod                                # working with the Earth as an ellipsod (WGS-84)
import geopy.distance                                  # moving along a known distance on the Earth's ellipsoid surface

import matplotlib.pyplot as plt                        # plotting a graph
from cartopy import crs as ccrs, feature as cfeature   # plotting a world map

# 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.index = location_df["name"]
location_df.head()

Determine the distance of a point to a great circle arc¶

The cross-track distance, sometimes known as cross track error, determines the distance from a point to a great circle arc and can be determined with vectors (typically simpler too).

geodesic = Geod(ellps="WGS84")
earth_radius = 6378137  # meters

Cross track distance: angular distance from point P to great circle path
Along track distance: angular distance along the great circle path from A to B before hitting a point that is closest to point P

Cross Track Distance¶

Distance of a point to a great circle arc is defined as:

dxt = asin( sin(δ13) ⋅ sin(θ13 − θ12) ) * R

(1)

δ13 (delta_13) is (angular) distance from start point to third point
θ13 (theta_13) is (initial) bearing from start point to third point
θ12 (theta_12) is (initial) bearing from start point to end point
R is the earth’s radius

d_xt = np.arcsin(np.sin(delta_13)*np.sin(theta_13 - theta_12)) * R

(2)

XTD =asin(sin(dist_AD)*sin(crs_AD-crs_AB))

(3)

Positive Cross-Track Distance: Point lies in the hemisphere to the left of the great circle
Negative Cross-Track Distance: Point lies in the hemiphere to the right of the great circle

If the point A is the N. or S. Pole replace crs_AD-crs_AB with lon_D-lon_B or lon_B-lon_D, respectively

def cross_track_distance(start_point=None, end_point=None, new_point=None):
    fwd_bearing_start_end, _, _ = 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"])
    fwd_bearing_start_new, _, distance_m_start_new = geodesic.inv(location_df.loc[start_point, "longitude"],
                                                                  location_df.loc[start_point, "latitude"],
                                                                  location_df.loc[new_point, "longitude"],
                                                                  location_df.loc[new_point, "latitude"]) 

    angular_distance_start_new = distance_m_start_new / earth_radius
    ct_distance = np.arcsin(np.sin(angular_distance_start_new) * np.sin(np.deg2rad(fwd_bearing_start_new - fwd_bearing_start_end))) * earth_radius
    return ct_distance

Along Track Distance¶

Distance along a great circle arc, closest to a point is defined as:

dat = acos(\frac{cos(δ13)}{cos(δxt)}) * R

(4)

δ13 (delta_13) is (angular) distance from start point to third point
δxt (delta_xt) is (angular) cross-track distance
R is the earth’s radius

d_at = np.arccos(\frac{np.cos(delta_13)}{np.cos(d_xt/R)} * R

(5)

ATD=acos(\frac{cos(dist_AD)}{cos(XTD)})

(6)

For very short distances (is less susceptible to rounding error):

ATD=asin(\sqrt{\frac{sin(dist_AD))^2 - (sin(XTD))^2 )}{cos(XTD)}}

(7)

def along_track_distance(start_point=None, end_point=None, new_point=None):
    crosst_distance = cross_track_distance(start_point, end_point, new_point)

    _, _, distance_m_start_new = geodesic.inv(location_df.loc[start_point, "longitude"],
                                              location_df.loc[start_point, "latitude"],
                                              location_df.loc[new_point, "longitude"],
                                              location_df.loc[new_point, "latitude"])  
    angular_distance_start_new = distance_m_start_new / earth_radius
    at_distance = np.arccos(np.cos(angular_distance_start_new) / np.cos(crosst_distance / earth_radius)) * earth_radius
    return at_distance

Generate Points at Intermediate Points along an Arc/Path¶

# Distance point along great circle path
def point_along_path(start_point=None, end_point=None, distance_m=None):
    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"])  
    origin = geopy.Point(location_df.loc[start_point, "latitude"],
                         location_df.loc[start_point, "longitude"])
    distance_to_move = geopy.distance.distance(
                            kilometers=distance_m / 1000)  # distance to move towards the next point
    final_position = distance_to_move.destination(origin, bearing=fwd_bearing)
    return (final_position.latitude, final_position.longitude)

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, reverse_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

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)

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

def plot_cross_track(start_point=None, end_point=None, new_point=None,
                     lon_west=-130, lon_east=-60,
                     lat_south=20, lat_north=60):
    # Set up world map plot
    fig = plt.subplots(figsize=(16, 10))
    projection_map = ccrs.PlateCarree()
    ax = plt.axes(projection=projection_map)
    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")

    # Cross-Track and Along-Track Distances
    ct_distance = cross_track_distance(start_point, end_point, new_point)
    print(f"Cross Track Distance: \n{ct_distance} meters ({ct_distance/1000} km)")

    at_distance = along_track_distance(start_point, end_point, new_point)
    print(f"Along Track Distance: \n{at_distance} meters ({at_distance/1000} km)\n")

    closest_point = point_along_path(start_point, end_point, at_distance)
    print(f"Closest Point To Point Along Great Circle Path:\n{closest_point}") 

    # Plot Latitude/Longitude Location
    great_circle_arc_pts = 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"],
                                      10)
    longitudes = [x[1] for x in great_circle_arc_pts] # longitude
    latitudes = [x[0] for x in great_circle_arc_pts] # latitude
    plt.plot(longitudes, latitudes, c="purple")
    plt.scatter(longitudes, latitudes, c="purple")

    cross_track_arc = arc_points(closest_point[0],
                                 closest_point[1],
                                 location_df.loc[new_point, "latitude"],
                                 location_df.loc[new_point, "longitude"],
                                   10)
    longitudes = [x[1] for x in cross_track_arc] # longitude
    latitudes = [x[0] for x in cross_track_arc] # latitude
    plt.plot(longitudes, latitudes, c="green")
    plt.scatter(longitudes, latitudes, c="green")

    # plot closest_point in red
    plt.scatter(closest_point[1], closest_point[0], c="red")

    plt.title(f"Closest Point {closest_point} from {start_point.title()}->{end_point.title()} to {new_point.title()}, Cross-Track Distance = {ct_distance/1000:4f} km")
    plt.show()

Positive Cross-Track Distance: Point lies in the hemisphere to the left of the great circle¶

plot_cross_track(start_point="boulder", end_point="greenwich", new_point="greenwich")

Cross Track Distance: 
0.0 meters (0.0 km)
Along Track Distance: 
7561763.794332366 meters (7561.763794332366 km)

Closest Point To Point Along Great Circle Path:
(51.49340000000001, 0.009800000000012687)

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

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

plot_cross_track(start_point="boulder", end_point="boston", new_point="cape canaveral")

Cross Track Distance: 
1593669.526094791 meters (1593.669526094791 km)
Along Track Distance: 
2076501.5510165778 meters (2076.5015510165777 km)

Closest Point To Point Along Great Circle Path:
(42.75525245755491, -80.62124342116076)

plot_cross_track(start_point="boulder", end_point="boston", new_point="arecibo",
                lat_south=15)

Cross Track Distance: 
2577081.0417989404 meters (2577.0810417989405 km)
Along Track Distance: 
3669432.5407487787 meters (3669.4325407487786 km)

Closest Point To Point Along Great Circle Path:
(41.142677689865174, -61.4898025092551)

Negative Cross-Track Distance: Point lies in the hemiphere to the right of the great circle¶

plot_cross_track(start_point="boulder", end_point="boston", new_point="redwoods")

Cross Track Distance: 
-744048.9243466797 meters (-744.0489243466797 km)
Along Track Distance: 
1409025.1954944504 meters (1409.0251954944504 km)

Closest Point To Point Along Great Circle Path:
(42.46116756301668, -88.74894120764374)

plot_cross_track(start_point="boulder", end_point="boston", new_point="greenwich",
                lon_east=15)

Cross Track Distance: 
-3381043.8402817464 meters (-3381.0438402817463 km)
Along Track Distance: 
7144535.346708467 meters (7144.535346708467 km)

Closest Point To Point Along Great Circle Path:
(28.39735920969611, -26.451361846821584)

Determine if a point lies on a great circle arc and path¶

With and without tolerances (in meters):

def is_point_on_arc(start_point=None, end_point=None,
                check_point=None, tolerance=0):
    # tolerance in meters
    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"])
    check_lon = np.deg2rad(location_df.loc[check_point, "longitude"])
    check_lat = np.deg2rad(location_df.loc[check_point, "latitude"])

    # Verify not meridian (longitude passes through the poles)
    if np.sin(lon1 - lon2) == 0:
        print("Invalid inputs: start/end points are meridians")
        return np.nan
    
    # 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 np.nan
    # account for tolerance based on cross-track distance from arc
    ct_distance = cross_track_distance(start_point, end_point, check_point)
    print(f"Cross-Track Distance = {ct_distance} meters")
    if np.abs(ct_distance) <= tolerance:
        return True

    # determine expected latitude
    num = np.sin(lat1)*np.cos(lat2)*np.sin(check_lon-lon2)-np.sin(lat2)*np.cos(lat1)*np.sin(check_lon-lon1)
    den = np.cos(lat1)*np.cos(lat2)*np.sin(lon1-lon2)
    new_lat = np.arctan(num/den)
    expected_lat = np.rad2deg(new_lat)
    return check_lat == expected_lat

Check if a point lies on a great circle arc¶

is_point_on_arc("boulder", "boston", "rockford", tolerance=0)

Cross-Track Distance = 18201.48035911659 meters

np.False_

plot_cross_track(start_point="boulder", end_point="boston", new_point="rockford")

Cross Track Distance: 
18201.48035911659 meters (18.20148035911659 km)
Along Track Distance: 
1378654.5186233742 meters (1378.654518623374 km)

Closest Point To Point Along Great Circle Path:
(42.434120910748035, -89.11630028269337)

plot_cross_track(start_point="boulder", end_point="boston", new_point="rockford", 
                 lon_west=-95, lon_east=-85,
                 lat_south=40, lat_north=45)

Cross Track Distance: 
18201.48035911659 meters (18.20148035911659 km)
Along Track Distance: 
1378654.5186233742 meters (1378.654518623374 km)

Closest Point To Point Along Great Circle Path:
(42.434120910748035, -89.11630028269337)

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

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

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

# increase tolerance to capture point
print("tolerance = 0")
print(is_point_on_arc("boulder", "boston", "rockford", tolerance=0))
print("\ntolerance >= cross-track distance")
print(is_point_on_arc("boulder", "boston", "rockford", tolerance=18202))

tolerance = 0
Cross-Track Distance = 18201.48035911659 meters
False

tolerance >= cross-track distance
Cross-Track Distance = 18201.48035911659 meters
True

Summary¶

Calculating and plotting the cross-track and along-trackd distance of a great circle arc/path and a point.

What’s next?¶

Determine when a great circle path crosses a given parallel and the maximum and minimum latitude coordinates of a great circle path.

Resources and references¶

Tutorials

Great Circle Arcs and Path

Tutorials

Great Circles and Parallels