Great Circle

Great Circle#

Overview#

This section covers great circle functions from NCL:

Functions#

area_poly_sphere#

NCL’s area_poly_sphere calculates the area enclosed by an arbitrary polygon on the sphere

Important Note

Coordinates should be within the valid latitude/longitude range (-90° to 90° and -180° to 180°) and be in clockwise order

Due to the shape of the Earth, the radius varies, but can be assumed to be a unit sphere with a radius of 6370997 m (based on the Clarke 1866 Authalic Sphere[1] model)

Grab and Go#

from pyproj import Geod

# Points in clockise order: Boulder, Boston, Houston
latitudes = [40.0150, 42.3601, 29.5518]  # degrees
longitudes = [-105.2705, -71.0589, -95.0982]  # degrees

geod = Geod(ellps="sphere")  # radius = 6370997 m
poly_area_m, _ = geod.polygon_area_perimeter(longitudes, latitudes)
poly_area_km2 = abs(poly_area_m) * 1e-6
poly_area_km2

1940855.984630871

css2c#

NCL’s css2c converts spherical (latitude/longitude) coordinates to Cartesian coordinates on a unit sphere

Grab and Go#

from astropy.coordinates.representation import UnitSphericalRepresentation
from astropy import units

lat = 40.0150
lon = -105.2705

spherical_coords = UnitSphericalRepresentation(lat=lat * units.deg, lon=lon * units.deg)
cart_coords = spherical_coords.to_cartesian()
print(f"X = {cart_coords.x.value}")
print(f"Y = {cart_coords.y.value}")
print(f"Z = {cart_coords.z.value}")

X = -0.20171369272651396
Y = -0.7388354627678497
Z = 0.6429881376224998

csc2s#

NCL’s csc2s converts Cartesian coordinates to spherical (latitude/longitude) coordinates on a unit sphere

Grab and Go#

from astropy.coordinates.representation import (
    CartesianRepresentation,
    SphericalRepresentation,
)
import numpy as np

x = -0.20171369272651396
y = -0.7388354627678497
z = 0.6429881376224998

cart_coords = CartesianRepresentation(x=x, y=y, z=z)
spherical_coords = cart_coords.represent_as(SphericalRepresentation)

# convert latitude/longitude from radians to degrees
lat_deg = np.rad2deg(spherical_coords.lat.value)
lon_deg = (
    np.rad2deg(spherical_coords.lon.value) + 180
) % 360 - 180  # keep longitude between -180 to 180

print(f"Latitude = {lat_deg}")
print(f"Longitude = {lon_deg}")

Latitude = 40.015
Longitude = -105.27049999999997

gc_onarc#

NCL’s gc_onarc determines if a point on the globe lies on a specified great circle arc as long as the angle between the two points along the great circle arc is not exactly 180 degrees (diametrically opposite, or antipodal).

Grab and Go#

import numpy as np


# Convert latitude and longitude points to Cartesian Points (see: css2c)
def latlon_to_cart(lat, lon):
    from astropy.coordinates.representation import UnitSphericalRepresentation
    from astropy import units

    spherical_coords = UnitSphericalRepresentation(
        lat=lat * units.deg, lon=lon * units.deg
    )
    cart_coords = spherical_coords.to_cartesian()
    return np.array([cart_coords.x, cart_coords.y, cart_coords.z])


pt_within = latlon_to_cart(40.0150, -105.2705)  # Boulder
vertex_a = latlon_to_cart(50.0150, -105.2705)  # Point exactly 10 degrees above Boulder
vertex_b = latlon_to_cart(30.0150, -105.2705)  # Point exactly 10 degrees below Boulder

# Determine if point lies along great circle arc
from uxarray.grid.arcs import point_within_gca

print(
    f"Boulder lies within the great circle arc = {point_within_gca(pt_within, vertex_a, vertex_b)}"
)

Boulder lies within the great circle arc = True

gc_qarea#

NCL’s gc_qarea calculates the area of a (four-sided) quadrilateral patch on the unit sphere

Important Note

Coordinates should be within the valid latitude/longitude range (-90° to 90° and -180° to 180°) and be in clockwise or counter-clockwise order

Grab and Go#

Area on the unit sphere#

NCL’s gc_qarea function finds the area of a quadrilateral patch on the unit sphere, a sphere with radius of 1

# Unit Sphere radius = 1
from pyproj import Geod

# Points in clockise order
latitudes = [90.0, 0.0, -90.0, 0.0]
longitudes = [0.0, -90.0, 0.0, 90.0]

# Adjust the radius of the spherical datum to describe the unit sphere
radius = 1

geod = Geod(a=radius)
poly_area, _ = geod.polygon_area_perimeter(longitudes, latitudes)
poly_area = abs(poly_area)
poly_area

6.283185307179586
Area on a spherical datum with Earth radius#

pyproj includes additional ellipsoid options, but the ellipsoid sphere treats the Earth as a sphere with an equal radius of 6370997 meters

# Normal Sphere: radius = 6370997 m
from pyproj import Geod

# Points in clockise order: Roughly Four Corners of Colorado
latitudes = [41.00488, 41.00203, 37.00540, 37.00051]  # degrees
longitudes = [-109.05001, -102.05348, -103.04633, -109.04720]  # degrees

geod = Geod(ellps="sphere")  # radius = 6370997 m
poly_area_m, _ = geod.polygon_area_perimeter(longitudes, latitudes)
poly_area_km2 = abs(poly_area_m) * 1e-6
poly_area_km2

250007.81588628562

gc_tarea#

NCL’s gc_tarea calculates the area of a triangular patch on the unit sphere

Important Note

Coordinates should be within the valid latitude/longitude range (-90° to 90° and -180° to 180°)

Grab and Go#

Area on the unit sphere#

NCL’s gc_tarea function finds the area of a triangular patch on the unit sphere, a sphere with radius of 1

# Unit Sphere radius = 1
from pyproj import Geod

# Latitude and Longitude points for one eighth surface area of a unit sphere
latitudes = [0.0, 0.0, 90.0]  # degrees
longitudes = [0.0, 90.0, 0.0]  # degrees

# Adjust the radius of the spherical datum to describe the unit sphere
radius = 1

geod = Geod(a=radius)
poly_area, _ = geod.polygon_area_perimeter(longitudes, latitudes)
poly_area = abs(poly_area)
poly_area

1.5707963267948966
Area on a spherical datum with Earth radius#

pyproj includes additional ellipsoid options, but the ellipsoid sphere treats the Earth as a sphere with an equal radius of 6370997 meters

# Normal Sphere: radius = 6370997 m
from pyproj import Geod

# Latitude and Longitude points for one eighth surface area of Earth
latitudes = [0.0, 0.0, 90.0]  # degrees
longitudes = [0.0, 90.0, 0.0]  # degrees

geod = Geod(ellps="sphere")  # radius = 6370997 m
poly_area_m, _ = geod.polygon_area_perimeter(longitudes, latitudes)
poly_area_km2 = abs(poly_area_m) * 1e-6
poly_area_km2

63757998.94347728

gc_latlon#

NCL’s gc_latlon calculates the great circle distance (true surface distance) between two points on the globe and interpolates points along the great circle

Important Note

Coordinates should be within the valid latitude/longitude range (-90° to 90° and -180° to 180°)

Grab and Go#

Great Circle Distance between Two Points#
from pyproj import Geod

# Latitude and Longitude points
lat1, lon1 = 40.0150, -105.2705  # Boulder
lat2, lon2 = 42.3601, -71.0589  # Boston

geodesic = Geod(ellps="sphere")
_, _, distance_m = geodesic.inv(lon1, lat1, lon2, lat2)

print(f"Distance = {distance_m / 1000} km")

Distance = 2855.332232201915 km
Interpolate Points along Great Circle Arc#

Important Note

NCL includes the starting and ending points along the great circle arc, while pyproj does not 
from pyproj import Geod

# Latitude and Longitude points
lat1, lon1 = 40.0150, -105.2705  # Boulder
lat2, lon2 = 42.3601, -71.0589  # Boston

geodesic = Geod(ellps="sphere")
pts = geodesic.npts(lon1, lat1, lon2, lat2, npts=10)

# pts returned in list of longtiude/latitude pairs
for latlon in pts:
    lat, lon = latlon[1], latlon[0]
    print(f"({lat}, {lon})")

(40.64111263335264, -102.32037915941086)
(41.1907927513628, -99.31785339056732)
(41.66124166643123, -96.26819663211273)
(42.049965273527796, -93.17759637679102)
(42.3548305568477, -90.05308125904574)
(42.574117071624066, -86.90241083568893)
(42.70656062076261, -83.73392987279547)
(42.751386588473196, -80.55639250241298)
(42.708330888194006, -77.3787644195555)
(42.577647195636175, -74.21001345000079)

Python Resources#

Additional Reading#

References:#

Contents