Ship on waves gif

Great Circles and Parallels¶

Overview¶

A valid great circle path (that is not a path around the equator) will cross a maximum and minimum latitude.

Determine the maximum latitude on a Great Circle Path
Determine the minimum latitude on a Great Great path
Determine when a great circle path crosses parallels (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       # 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 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.head()

location_df.index = location_df["name"]

Maximum Latitude on a Great Circle Path¶

We have previously determined an equation to derive a great circle path from intermediate points from two points on a great circle arc.

Without additional calculations, we can use a list of points along the great circle path to find the maximum location of the maximum and minimum.

By default, the equation below will determine 360 points along longitude, so the output will only have a resolution of 1 degree. However, by defining the longitude with more points, the resolution increases.

# Generate Latitude Coordinates based on Longitude Coordinates
def generate_latitude_along_gc(start_point=None, end_point=None, number_of_lon_pts=360):
    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"])

    # 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

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=360)
print(f"Max Latitude (within 1 degree): {max(lat_lon, key=lambda x:x[0])}")

Max Latitude (within 1 degree): (np.float64(42.750406941471915), np.float64(-81.0))

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=720)
print(f"Max Latitude (within 0.5 degree): {max(lat_lon, key=lambda x:x[0])}")

Max Latitude (within 0.5 degree): (np.float64(42.751388471834524), np.float64(-80.5))

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=1080)
print(f"Max Latitude (within 0.3 degree): {max(lat_lon, key=lambda x:x[0])}")

Max Latitude (within 0.3 degree): (np.float64(42.751302958796096), np.float64(-80.66666666666667))

Plot Maximum¶

def plot_coordinate_max_min(great_circle_pts=None,
                            max_coord=None, min_coord=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.STATES, edgecolor="black")
    
    # Plot Great Circle Latitude/Longitude Location
    longitudes = [x[1] for x in great_circle_pts] # longitude
    latitudes = [x[0] for x in great_circle_pts] # latitude
    plt.plot(longitudes, latitudes)

    # Overly Max/Min Coordinates
    if max_coord is not None:
        plt.scatter([max_coord[1]], [max_coord[0]], s=100, c="red")
    if min_coord is not None:
        plt.scatter([min_coord[1]], [min_coord[0]], s=100, c="green")
    
    # Setup Axis Limits and Title/Labels
    plt.title(title)
    plt.show()

gc_lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=360)
max_lat_lon = max(gc_lat_lon, key=lambda x:x[0])
plot_coordinate_max_min(great_circle_pts=gc_lat_lon,
                            max_coord=max_lat_lon,
                            title=f"Max Latitude located at {max_lat_lon}")

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

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

gc_lat_lon = generate_latitude_along_gc("boulder", "houston", number_of_lon_pts=360)
max_lat_lon = max(gc_lat_lon, key=lambda x:x[0])
plot_coordinate_max_min(great_circle_pts=gc_lat_lon,
                            max_coord=max_lat_lon,
                            title=f"Max Latitude located at {max_lat_lon}")

Maximumn Latitude from Clairaut’s Formula¶

Clairaut’s Formula (Clairaut’s equation or Clairaut’s relation) is a differential equation which defines the relationship between the latitude, φ, and the true course (bearing, θ) where:

sin(θ) * cos(φ) = \text{constant}

(1)

So, for any two points (A and B) along the great circle:

sin(θA) * cos(φA) = sin(θB) * cos(φB)

(2)

So, to solve for the maximum latitude the true course should be when 90 and 270 degrees on the unit sphere where for any bearing/latitude along the great circle:

\text{max latitude} = acos(|sin(θ) * cos(φ)|)

(3)

For the purpose of this example, we will use pyproj geodesic to determine the bearing based on a great circle arc, but consult previous sections if you want to determine bearing mathetically based on the unit sphere instead of the ellipsoid.

Important Note

Clairaut’s Formula works from unit sphere, and as a result, is subject to errors (about 3%, about +/- 11 degrees).

Ed Williams: Clairaut’s Formula

def clairaut_formula_max(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"])
    # Clairaut
    start_lat = np.deg2rad(location_df.loc[start_point, "longitude"])
    max_lat = np.arccos(np.abs(np.sin(fwd_bearing) * np.cos(start_lat)))
    return np.rad2deg(max_lat)

max_lat = clairaut_formula_max("boulder", "boston")
print(f"Max latitude from Boulder to Boston: {max_lat}")

Max latitude from Boulder to Boston: 75.50718325253314

Minimum Latitude on a Great Circle Path¶

Antipodal Point of Max (TODO)¶

def antipodal(latitude=None, longitude=None):
    anti_lat = -1 * latitude
    if longitude > 0:
        anti_lon = longitude - 180
    else:
        anti_lon = longitude + 180
    return (anti_lat, anti_lon)

Like finding maximum from a list of great circle path, the smallest latitude can be found by analysing the list for the smallest latitude point.

gc_lat_lon = generate_latitude_along_gc("boulder", "houston", number_of_lon_pts=360)
max_lat_lon = max(gc_lat_lon, key=lambda x:x[0])
print(max_lat_lon)
antipodal(max_lat_lon[0], max_lat_lon[1])

(np.float64(59.410929434369436), np.float64(-166.0))

(np.float64(-59.410929434369436), np.float64(14.0))

Minimum Latitude along Great Circle Path¶

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=360)
print(f"Min Latitude (within 1 degree): {min(lat_lon, key=lambda x:x[0])}")

Min Latitude (within 1 degree): (np.float64(-42.75040694147194), np.float64(99.0))

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=720)
print(f"Min Latitude (within 0.5 degree): {min(lat_lon, key=lambda x:x[0])}")

Min Latitude (within 0.5 degree): (np.float64(-42.75138847183453), np.float64(99.5))

lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=1080)
print(f"Min Latitude (within 0.3 degree): {min(lat_lon, key=lambda x:x[0])}")

Min Latitude (within 0.3 degree): (np.float64(-42.7513029587961), np.float64(99.33333333333331))

gc_lat_lon = generate_latitude_along_gc("boulder", "boston", number_of_lon_pts=360)
min_lat_lon = min(gc_lat_lon, key=lambda x:x[0])
plot_coordinate_max_min(great_circle_pts=gc_lat_lon,
                            min_coord=min_lat_lon,
                            title=f"Min Latitude located at {min_lat_lon}")

Maximumn Latitude from Clairaut’s Formula¶

To solve for the minimum, the true course should be when 0 and 180 degrees on the unit sphere where for any bearing/latitude along the great circle:

\text{min latitude} = asin(|sin(θ) * cos(φ)|)

(4)

The southernmost point is the antipode to the northernmost (max) latitude.

def clairaut_formula_min(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"])
    # Clairaut Formula
    start_lat = np.deg2rad(location_df.loc[start_point, "longitude"])
    min_lat = np.arcsin(np.abs(np.cos(fwd_bearing) * np.sin(start_lat)))
    return np.rad2deg(min_lat)

min_lat = clairaut_formula_min("boulder", "boston")
print(f"Min latitude along great circle path from Boulder to Boston: {min_lat}")

Min latitude along great circle path from Boulder to Boston: 17.49699780715814

Determine when great circle path cross parallels (TODO)¶

Determine the longitude when a great circle crosses a given latitude parrellel.

Summary¶

Determine the coordinates when a great circle path crosses a specific parallel as well as the maximumn and minimum latitude coordinates.

What’s next?¶

Intersections of Great Circles.

Resources and references¶

Tutorials

Great Circles and a Point

Tutorials

Intersections of Great Circles