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:241: 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:2155, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2152     # we do this instead of `self.figure.draw_without_rendering`
   2153     # so that we can inject the orientation
   2154     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2155         self.figure.draw(renderer)
   2156 if bbox_inches:
   2157     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:524, in GeoAxes.draw(self, renderer, **kwargs)
    519         self.imshow(img, extent=extent, origin=origin,
    520                     transform=factory.crs, *factory_args[1:],
    521                     **factory_kwargs)
    522 self._done_img_factory = True
--> 524 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:3216, in _AxesBase.draw(self, renderer)
   3213 if artists_rasterized:
   3214     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3216 mimage._draw_list_compositing_images(
   3217     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3219 renderer.close_group('axes')
   3220 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:185, in FeatureArtist.draw(self, renderer)
    180     geoms = self._feature.geometries()
    181 else:
    182     # For efficiency on local maps with high resolution features (e.g
    183     # from Natural Earth), only create paths for geometries that are
    184     # in view.
--> 185     geoms = self._feature.intersecting_geometries(extent)
    187 stylised_paths = {}
    188 # Make an empty placeholder style dictionary for when styler is not
    189 # used.  Freeze it so that we can use it as a dict key.  We will need
    190 # 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:309, in NaturalEarthFeature.intersecting_geometries(self, extent)
    302 """
    303 Returns an iterator of shapely geometries that intersect with
    304 the given extent.
    305 The extent is assumed to be in the CRS of the feature.
    306 If extent is None, the method returns all geometries for this dataset.
    307 """
    308 self.scaler.scale_from_extent(extent)
--> 309 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:294, in NaturalEarthFeature.geometries(self)
    290 if key not in _NATURAL_EARTH_GEOM_CACHE:
    291     path = shapereader.natural_earth(resolution=self.scale,
    292                                      category=self.category,
    293                                      name=self.name)
--> 294     geometries = tuple(shapereader.Reader(path).geometries())
    295     _NATURAL_EARTH_GEOM_CACHE[key] = geometries
    296 else:

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

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1483, in Reader.iterShapes(self, bbox)
   1479 if self.numShapes:
   1480     # Iterate exactly the number of shapes from shx header
   1481     for i in xrange(self.numShapes):
   1482         # MAYBE: check if more left of file or exit early? 
-> 1483         shape = self.__shape(oid=i, bbox=bbox)
   1484         if shape:
   1485             yield shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1340, in Reader.__shape(self, oid, bbox)
   1338 # Read points - produces a list of [x,y] values
   1339 if nPoints:
-> 1340     flat = unpack("<%sd" % (2 * nPoints), f.read(16*nPoints))
   1341     record.points = list(izip(*(iter(flat),) * 2))
   1342 # Read z extremes and values

error: unpack requires a buffer of 736 bytes

<Figure size 1600x1000 with 2 Axes>

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)

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:241: 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="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:241: 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:241: 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:241: 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