Great Circle#
Overview#
This section covers great circle functions from NCL:
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 orderDue 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 orderGrab 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°)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°)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 notfrom 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)