Great Circle#
Warning
This is not meant to be a standalone notebook. This notebook is part of the process we have for adding entries to the NCL Index and is not meant to be used as tutorial or example code.
Functions covered#
NCL code#
; area_poly_sphere
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml
; Boulder, Boston, Houston (3 Sides)
lat0 = (/40.0150, 42.3601, 29.5518/)
lon0 = (/-105.2705, -71.0589, -95.0982/)
print(gc_clkwise(lat0, lon0))
; (0) True
print(area_poly_sphere(lat0, lon0, 6370.997))
; (0) 1940856
; Roughly the Four Corners of Colorado (4 Sides)
lat1 = (/41.00488, 41.00203, 37.00540, 37.00051/)
lon1 = (/-109.05001, -102.05348, -103.04633, -109.04720/)
print(gc_clkwise(lat1, lon1))
; (0) True
print(area_poly_sphere(lat1, lon1, 6370.997))
; (0) 250007.6
; Caltech, Alberta, Greenwich, Paris, Madrid (5 Sides)
lat2 = (/34.1377, 53.9333, 51.4934, 48.8575, 40.4167/)
lon2 = (/-118.1253, -116.5765, 0.0098, 2.3514, -3.7033/)
print(gc_clkwise(lat2, lon2))
; (0) True
print(area_poly_sphere(lat2, lon2, 6370.997))
; (0) 1.16348e+07
;; Edge Cases:
; Crossing the Equator (0 degrees Latitude)
; Libreville, Moanda, San Tome and Principe
lat3 = (/0.4078, -5.9230, 0.1864/)
lon3 = (/9.4403, 12.3636, 6.6131/)
print(gc_clkwise(lat3, lon3))
; (0) True
print(area_poly_sphere(lat3, lon3, 6370.997))
; (0) 114894.8
; Crossing the Prime Meridian (0 Degrees Longitude)
; Dublin, Amsterdam, Fishguard
lat4 = (/53.3498, 52.3676, 51.9939/)
lon4 = (/-6.2603, 4.9041, -4.9760/)
print(gc_clkwise(lat4, lon4))
; (0) True
print(area_poly_sphere(lat4, lon4, 6370.997))
; (0) 54450.39
; Half of the World
lat5 = (/90.0, 0.0, -90.0, 0.0/)
lon5 = (/0.0, -90.0, 0.0, 90.0/)
print(gc_clkwise(lat5, lon5))
; (0) True
print(area_poly_sphere(lat5, lon5, 6370.997))
; (0) 2.55032e+08
; Single Point -> Invalid NCL
lat6 = (/40, 40, 40, 40/)
lon6 = (/-105, -105, -105, -105/)
print(gc_clkwise(lat6, lon6))
; (0) True
print(area_poly_sphere(lat6, lon6, 6370.997))
; (0) -1.27516e+08
; Single Degree
lat7 = (/40, 40, 41, 41/)
lon7 = (/-105, -106, -106, -105/)
print(gc_clkwise(lat7, lon7))
; (0) True
print(area_poly_sphere(lat7, lon7, 6370.997))
;(0) 9401.705
; css2c
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml
; ncl -n css2c.ncl >> css2c_output.txt
print("Latitude (Degree), Longitude (Degree), Cartesian X, Cartesian Y, Cartesian Z")
do lat=-90,90
do lon=-180,180
begin
cart = css2c(lat, lon)
print (lat + "," + lon + "," + cart(0,0) + "," + cart(1,0) + "," + cart(2,0))
end
end do
end do
; csc2s
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml
; ncl -n csc2s.ncl >> csc2s_output.txt
print("Input Latitude (Degree), Input Longitude (Degree), Cartesian X, Cartesian Y, Cartesian Z, Output Latitude (Degree), Output Longitude (Degree)")
do lat=-90,90
do lon=-180,180
begin
cart = css2c(lat, lon)
; determine a list of xyz coordinates based on lat/lon
x = cart(0,0)
y = cart(1,0)
z = cart(2,0)
sph = csc2s(x, y, z)
print(lat + "," + lon + "," + x + "," + y + "," + z + "," + sph(0,0) + "," + sph(1,0))
end
end do
end do
; gc_onarc
; ncl -n gc_onarc.ncl >> gc_onarc_output.txt
; Arc above/below latitude with point in Boulder with same longitude
; True
is_on_arc = gc_onarc(40.0150, -105.2705, (/50.0150, 30.0150/), (/-105.2705, -105.2705/))
print(is_on_arc)
; Arc from Denver to Boston, point with Boulder
; False
is_on_arc = gc_onarc(40.0150, -105.2705, (/39.73915, 42.35843/), (/-104.9847, -71.05977/))
print(is_on_arc)
; All three points in the same position
; True
is_on_arc = gc_onarc(40.0150, -105.2705, (/40.0150, 40.0150/), (/-105.2705, -105.2705/))
print(is_on_arc)
; Anti-meridian
; False
is_on_arc = gc_onarc(24.86658349036834, 151.60735730237568, (/55.182644873824785, 49.05673786129486/), (/30.6703372467811, 179.59874545906715/))
print(is_on_arc)
; gc_qarea
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/gc_qarea.shtml
; Roughly the Four Corners of Colorado (4 Sides)
lat = (/41.00488, 41.00203, 37.00540, 37.00051/)
lon = (/-109.05001, -102.05348, -103.04633, -109.04720/)
print(gc_clkwise(lat, lon))
; (0) True
print(gc_qarea(lat, lon)*6370.997^2)
; (0) 250007.7
; Boulder, Boston, Cape Canaveral, Houston
lat = (/40.0150, 42.3601, 28.3968, 29.5518/)
lon = (/-105.2705, -71.0589, -80.6057, -95.0982/)
print(gc_clkwise(lat, lon))
; (0) True
print(gc_qarea(lat, lon)*6370.997^2)
; (0) 3154871
;; Edge Cases
; Crossing the Equator (0 degrees Latitude)
; Boulder, Boston, Montevideo, Quito
lat = (/40.0150, 42.3601, -34.5420, -0.1312/)
lon = (/-105.2705, -71.0589, -56.1103, -78.3045/)
print(gc_clkwise(lat, lon))
; (0) True
print(gc_qarea(lat, lon)*6370.997^2)
; (0) 1.507316e+07
; Crossing the Prime Meridian (0 Degrees Longitude)
; Dublin, Norwich, London, Cardiff
lat = (/53.2100, 52.3743, 51.3026, 51.2854/)
lon = (/-6.1537, 1.1734, -0.739, -3.1045/)
print(gc_clkwise(lat, lon))
; (0) True
print(gc_qarea(lat, lon)*6370.997^2)
; (0) 54846.59
; Half of the World
lat = (/90.0, 0.0, -90.0, 0.0/)
lon = (/0.0, -90.0, 0.0, 90.0/)
print(gc_clkwise(lat, lon))
; (0) True
print(gc_qarea(lat, lon)*6370.997^2)
; (0) 2.55032e+08
; gc_tarea
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/gc_tarea.shtml
; 1/8th the surface of Earth (North)
lat = (/0.0, 0.0, 90.0/)
lon = (/0.0, 90.0, 0.0/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 6.3758e+07
; 1/8th the surface of Earth (South)
lat = (/0.0, 0.0, -90.0/)
lon = (/0.0, 90.0, 0.0/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 6.3758e+07
; Triangle Across United (Redwoods, Boston, Houston)
lat = (/41.4017, 42.3601, 29.5518/)
lon = (/-124.0417, -71.0589, -95.0982/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 3782549
;; Edge Cases
; Boulder, Boston, Montevideo (Crosses 0 degrees Latitude)
lat = (/40.0150, 42.3601, -34.5420/)
lon = (/-105.2705, -71.0589, -56.1103/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 1.493506e+07
; Boulder, Cairo, Houston (Crosses 0 degrees Longitude)
lat = (/40.0150, 30.240, 29.5518/)
lon = (/-105.2705, 31.149, -95.0982/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 1.114434e+07
; Equator and Pole (two 90 degree spherical angles)
lat = (/0, 0, 90/)
lon = (/-105, 0, 0/)
tri_area = gc_tarea(lat, lon)
print(tri_area*6370.997^2)
; (0) 7.438434e+07
; gc_latlon
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/gc_latlon.shtml
;; Calculate great circle distances between two points
; Distance between Equator and North Pole
lat1 = 0
lon1 = 0
lat2 = 90
lon2 = 0
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
; (0) 10007.89
; Distance between Sacramento and Albany
lat1 = 38.54623
lon1 = -121.42660
lat2 = 42.66575
lon2 = -73.79901
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
; (0) 3993.532
; Distance between Boulder and Boston
lat1 = 40.0150
lon1 = -105.2705
lat2 = 42.3601
lon2 = -71.0589
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
;(0) 2855.432
;; Edge Cases
; Distance to Nairobi and Cairo (Crosses 0 degrees Latitude)
lat1 = -1.1711
lon1 = 36.4902
lat2 = 30.240
lon2 = 31.149
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
;(0) 3538.456
; Distance to Dublin and London (Crosses 0 degrees Longitude)
lat1 = 53.21
lon1 = -5.1537
lat2 = 51.3026
lon2 = 0.739
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
;(0) 453.5349
; Distance from North to South Pole
lat1 = 90
lon1 = 0
lat2 = -90
lon2 = 0
gc_distance_km = gc_latlon(lat1,lon1,lat2,lon2,2,4)
print(gc_distance_km)
;(0) 20015.78
;; Interpolate points along great circle arc between two points
; Distance between Equator and North Pole
lat1 = 0
lon1 = 0
lat2 = 90
lon2 = 0
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; units : degrees
; gclat : ( 0, 10, 20, 30, 40, 50, 60, 70, 80, 90 )
; gclon : ( 0, 6.186176e-16, 1.276937e-15, 2.025549e-15, 2.943859e-15, 4.181094e-15, 6.076649e-15, 9.639125e-15, 1.989687e-14, 0 )
; spacing : 10
; (0) 90
; Distance between Sacramento and Albany
lat1 = 38.54623
lon1 = -121.42660
lat2 = 42.66575
lon2 = -73.79901
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; spacing : 3.990384
; gclat : ( 38.54623, 39.95334, 41.15266, 42.12738, 42.86257, 43.34622, 43.56999, 43.52992, 43.22673, 42.66575 )
; gclon : ( 238.5734, 243.3962, 248.4062, 253.585, 258.9052, 264.3314, 269.8211, 275.3275, 280.8027, 286.201 )
; units : degrees
; (0) 35.91346
; Distance between Boulder and Boston
lat1 = 40.0150
lon1 = -105.2705
lat2 = 42.3601
lon2 = -71.0589
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; spacing : 2.853181
; gclat : ( 40.015, 40.76997, 41.40982, 41.92961, 42.32515, 42.59314, 42.73126, 42.73831, 42.61423, 42.3601 )
; gclon : ( 254.7295, 258.3425, 262.0321, 265.7881, 269.5983, 273.4489, 277.3247, 281.2095, 285.087, 288.9411 )
; units : degrees
; (0) 25.67863
;; Edge Cases
; Distance to Nairobi and Cairo (Crosses 0 degrees Latitude)
lat1 = -1.1711
lon1 = 36.4902
lat2 = 30.240
lon2 = 31.149
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; spacing : 3.535667
; gclat : ( -1.1711, 2.3232, 5.817294, 9.310863, 12.80358, 16.29507, 19.78496, 23.27277, 26.758, 30.24 )
; gclon : ( 36.4902, 35.95082, 35.40877, 34.85993, 34.3, 33.72431, 33.12772, 32.50438, 31.84749, 31.149 )
; units : degrees
; (0) 31.821
; Distance to Dublin and London (Crosses 0 degrees Longitude)
lat1 = 53.21
lon1 = -5.1537
lat2 = 51.3026
lon2 = 0.739
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; spacing : 0.4531774
; gclat : ( 53.21, 53.01284, 52.81184, 52.60707, 52.39858, 52.18644, 51.97071, 51.75145, 51.52873, 51.3026 )
; gclon : ( 354.8463, 355.5261, 356.1996, 356.8669, 357.5279, 358.1826, 358.8311, 359.4733, 0.1092852, 0.739 )
; units : degrees
; (0) 4.078597
; Distance from North to South Pole
lat1 = 90
lon1 = 0
lat2 = -90
lon2 = 0
gc_distance = gc_latlon(lat1,lon1,lat2,lon2,10,2)
print(gc_distance)
; spacing : 20
; gclat : ( 90, 70, 50, 30, 10, -10, -30, -50, -70, -90 )
; gclon : ( 0, 3.508355e-15, 3.508355e-15, 3.508355e-15, 3.508355e-15, 3.508355e-15, 3.508355e-15, 3.508355e-15, 3.508355e-15, 0 )
; units : degrees
; (0) 180
Python Functionality#
area_poly_sphere#
import math
ncl_results = {}
poly_name = [
"Boulder, Boston, Houston",
"Four Corners of Colorado",
"Caltech, Alberta, Greenwich, Paris, Madrid",
"Crossing the Equator",
"Crossing the Prime Meridian",
"Half of the World",
"Single Point -> Invalid NCL",
"Single Degree",
]
ncl_lat = [
[40.0150, 42.3601, 29.5518],
[41.00488, 41.00203, 37.00540, 37.00051],
[34.1377, 53.9333, 51.4934, 48.8575, 40.4167],
[0.4078, -5.9230, 0.1864],
[53.3498, 52.3676, 51.9939],
[90.0, 0.0, -90.0, 0.0],
[40, 40, 40, 40],
[40, 40, 41, 41],
]
ncl_lon = [
[-105.2705, -71.0589, -95.0982],
[-109.05001, -102.05348, -103.04633, -109.04720],
[-118.1253, -116.5765, 0.0098, 2.3514, -3.7033],
[9.4403, 12.3636, 6.6131],
[-6.2603, 4.9041, -4.9760],
[0.0, -90.0, 0.0, 90.0],
[-105, -105, -105, -105],
[-105, -106, -106, -105],
]
ncl_results[poly_name[0]] = 1940856
ncl_results[poly_name[1]] = 250007.6
ncl_results[poly_name[2]] = 11634800
ncl_results[poly_name[3]] = 114894.8
ncl_results[poly_name[4]] = 54450.39
ncl_results[poly_name[5]] = 255032000
ncl_results[poly_name[6]] = -127516000
ncl_results[poly_name[7]] = 9401.705
from pyproj import Geod
python_results = {}
geod = Geod(ellps="sphere") # radius = 6370997 m
for i in range(len(poly_name)):
poly_area_m, _ = geod.polygon_area_perimeter(ncl_lon[i], ncl_lat[i])
poly_area_km2 = abs(poly_area_m) * 1e-6
python_results[poly_name[i]] = poly_area_km2
css2c#
import geocat.datafiles as gcd
import numpy as np
from astropy.coordinates.representation import UnitSphericalRepresentation
from astropy import units
css2c_data = gcd.get('applications_files/ncl_outputs/css2c_output.txt')
css2c_data = np.loadtxt(css2c_data, delimiter=',', skiprows=6)
lat_lon = tuple(map(tuple, (css2c_data[::, 0:2])))
cart_values = tuple(css2c_data[::, 2:])
ncl_css2c = dict(zip(lat_lon, cart_values))
Downloading file 'applications_files/ncl_outputs/css2c_output.txt' from 'https://github.com/NCAR/geocat-datafiles/raw/main/applications_files/ncl_outputs/css2c_output.txt' to '/home/runner/.cache/geocat'.
## Collect Latitude and Longtiude Pairs
lat_lon = []
for lat in range(-90, 90 + 1):
for lon in range(-180, 180 + 1):
lat_lon.append((lat, lon))
## Calculate Cartesian coordinates
astropy_css2c = {}
for i, pair in enumerate(lat_lon):
lat, lon = pair
spherical_coords = UnitSphericalRepresentation(
lat=lat * units.deg, lon=lon * units.deg
)
cart_coords = spherical_coords.to_cartesian()
astropy_css2c[pair] = cart_coords.xyz.value
csc2s#
import geocat.datafiles as gcd
from astropy.coordinates.representation import (
CartesianRepresentation,
SphericalRepresentation,
)
import numpy as np
csc2s_data = gcd.get('applications_files/ncl_outputs/csc2s_output.txt')
csc2s_data = np.loadtxt(csc2s_data, delimiter=',', skiprows=6)
input_lat_lon = tuple(map(tuple, csc2s_data[::, 0:2]))
cart_values = tuple(map(tuple, (csc2s_data[::, 2:5])))
output_lat_lon = tuple(map(tuple, (csc2s_data[::, 5:])))
ncl_csc2s = dict(zip(input_lat_lon, cart_values))
ncl_csc2s_input_output = dict(zip(input_lat_lon, output_lat_lon))
Downloading file 'applications_files/ncl_outputs/csc2s_output.txt' from 'https://github.com/NCAR/geocat-datafiles/raw/main/applications_files/ncl_outputs/csc2s_output.txt' to '/home/runner/.cache/geocat'.
## Calculate Spherical coordinates
def spherical_cart(x, y, z):
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
return (lat_deg, lon_deg)
astropy_csc2s = {}
for xyz in cart_values:
x, y, z = xyz
lat_lon = spherical_cart(x, y, z)
astropy_csc2s[lat_lon] = tuple(xyz)
gc_onarc#
ncl_gc_onarc = {}
# Arc above/below latitude with point in Boulder with same longitude = True
ncl_lat = (40.0150, 50.0150, 30.0150)
ncl_lon = (-105.2705, -105.2705, -105.2705)
ncl_gc_onarc[(ncl_lat, ncl_lon)] = True
# Arc from Denver to Boston, point with Boulder
ncl_lat = (40.0150, 39.73915, 42.35843)
ncl_lon = (-105.2705, -104.9847, -71.05977)
ncl_gc_onarc[(ncl_lat, ncl_lon)] = False
# All three points in the same position
ncl_lat = (40.0150, 40.0150, 40.0150)
ncl_long = (-105.2705, -105.2705, -105.2705)
ncl_gc_onarc[(ncl_lat, ncl_long)] = True
# Anti-meridian
ncl_lat = (24.86658349036834, 55.182644873824785, 49.05673786129486)
ncl_lon = (151.60735730237568, 30.6703372467811, 179.59874545906715)
ncl_gc_onarc[(ncl_lat, ncl_lon)] = False
import numpy as np
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])
python_onarc = {}
from uxarray.grid.arcs import point_within_gca
for ncl_latlon in ncl_gc_onarc.keys():
pt_within = latlon_to_cart(ncl_latlon[0][0], ncl_latlon[1][0])
vertex_a = latlon_to_cart(ncl_latlon[0][1], ncl_latlon[1][1])
vertex_b = latlon_to_cart(ncl_latlon[0][2], ncl_latlon[1][2])
python_onarc[ncl_latlon] = point_within_gca(pt_within, vertex_a, vertex_b)
gc_qarea#
ncl_gc_qarea = {}
# Roughly the Four Corners of Colorado
ncl_lat = (41.00488, 41.00203, 37.00540, 37.00051)
ncl_lon = (-109.05001, -102.05348, -103.04633, -109.04720)
ncl_gc_qarea[(ncl_lat, ncl_lon)] = 250007.7
# Boulder, Boston, Cape Canaveral, Houston
ncl_lat = (40.0150, 42.3601, 28.3968, 29.5518)
ncl_lon = (-105.2705, -71.0589, -80.6057, -95.0982)
ncl_gc_qarea[(ncl_lat, ncl_lon)] = 3154871
# Crossing the Equator (0 degrees Latitude) = Boulder, Boston, Montevideo, Quito
ncl_lat = (40.0150, 42.3601, -34.5420, -0.1312)
ncl_long = (-105.2705, -71.0589, -56.1103, -78.3045)
ncl_gc_qarea[(ncl_lat, ncl_long)] = 15073160 # 1.507316e+07
# Crossing the Prime Meridian (0 Degrees Longitude)= Dublin, Norwich, London, Cardiff
ncl_lat = (53.2100, 52.3743, 51.3026, 51.2854)
ncl_lon = (-6.1537, 1.1734, -0.739, -3.1045)
ncl_gc_qarea[(ncl_lat, ncl_lon)] = 54846.59
# Half of the World
ncl_lat = (90.0, 0.0, -90.0, 0.0)
ncl_lon = (0.0, -90.0, 0.0, 90.0)
ncl_gc_qarea[(ncl_lat, ncl_lon)] = 255032000 # 2.55032e+08
from pyproj import Geod
python_gc_qarea = {}
geod = Geod(ellps="sphere") # radius = 6370997 m
for lat_lon in ncl_gc_qarea.keys():
lat, lon = lat_lon
poly_area_m, _ = geod.polygon_area_perimeter(lon, lat)
poly_area_km2 = abs(poly_area_m) * 1e-6
python_gc_qarea[lat_lon] = poly_area_km2
gc_tarea#
ncl_gc_tarea = {}
# 1/8th the surface of Earth (North)
ncl_lat = (0.0, 0.0, 90.0)
ncl_long = (0.0, 90.0, 0.0)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 63758000 # 6.3758e+07
# 1/8th the surface of Earth (South)
ncl_lat = (0.0, 0.0, -90.0)
ncl_long = (0.0, 90.0, 0.0)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 63758000 # 6.3758e+07
# Triangle Across United (Redwoods, Boston, Houston)
ncl_lat = (41.4017, 42.3601, 29.5518)
ncl_long = (-124.0417, -71.0589, -95.0982)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 3782549
# Boulder, Boston, Montevideo (Crosses 0 degrees Latitude)
ncl_lat = (40.0150, 42.3601, -34.5420)
ncl_long = (-105.2705, -71.0589, -56.1103)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 14935060 # 1.493506e+07
# Boulder, Cairo, Houston (Crosses 0 degrees Longitude)
ncl_lat = (40.0150, 30.240, 29.5518)
ncl_long = (-105.2705, 31.149, -95.0982)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 11144340 # 1.114434e+07
# BEquator and Pole (two 90 degree spherical angles)
ncl_lat = (0, 0, 90)
ncl_long = (-105, 0, 0)
ncl_gc_tarea[(ncl_lat, ncl_long)] = 74384340 # 7.438434e+07
from pyproj import Geod
python_gc_tarea = {}
geod = Geod(ellps="sphere") # radius = 6370997 m
for lat_lon in ncl_gc_tarea.keys():
lat, lon = lat_lon
poly_area_m, _ = geod.polygon_area_perimeter(lon, lat)
poly_area_km2 = abs(poly_area_m) * 1e-6
python_gc_tarea[lat_lon] = poly_area_km2
gc_latlon#
Great Circle Distance#
ncl_gc_latlon_dist = {}
# Distance between Equator and North Pole
lat1 = 0
lon1 = 0
lat2 = 90
lon2 = 0
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 10007.89
# Distance between Sacramento and Albany
lat1 = 38.54623
lon1 = -121.42660
lat2 = 42.66575
lon2 = -73.79901
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 3993.532
# Distance between Boulder and Boston
lat1 = 40.0150
lon1 = -105.2705
lat2 = 42.3601
lon2 = -71.0589
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 2855.432
# Distance bewtween Nairobi and Cairo (Crosses 0 degrees Latitude)
lat1 = -1.1711
lon1 = 36.4902
lat2 = 30.240
lon2 = 31.149
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 3538.456
# Distance between Dublin and London (Crosses 0 degrees Longitude)
lat1 = 53.21
lon1 = -5.1537
lat2 = 51.3026
lon2 = 0.739
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 453.5349
# Distance between North to South Pole
lat1 = 90
lon1 = 0
lat2 = -90
lon2 = 0
ncl_gc_latlon_dist[(lat1, lon1, lat2, lon2)] = 20015.78
from pyproj import Geod
python_gc_latlon_dist = {}
for latlon in ncl_gc_latlon_dist.keys():
geodesic = Geod(ellps="sphere")
_, _, distance_m = geodesic.inv(latlon[1], latlon[0], latlon[3], latlon[2])
distance_km = distance_m / 1000
python_gc_latlon_dist[latlon] = distance_km
Interpolated Points#
ncl_gc_latlon_npts = {}
# Distance between Equator and North Pole
lat1 = 0
lon1 = 0
lat2 = 90
lon2 = 0
lon = [
0,
6.186176e-16,
1.276937e-15,
2.025549e-15,
2.943859e-15,
4.181094e-15,
6.076649e-15,
9.639125e-15,
1.989687e-14,
0,
]
lat = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
# Distance between Sacramento and Albany
lat1 = 38.54623
lon1 = -121.42660
lat2 = 42.66575
lon2 = -73.79901
lat = [
38.54623,
39.95334,
41.15266,
42.12738,
42.86257,
43.34622,
43.56999,
43.52992,
43.22673,
42.66575,
]
lon = [
238.5734,
243.3962,
248.4062,
253.585,
258.9052,
264.3314,
269.8211,
275.3275,
280.8027,
286.201,
]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
# Distance between Boulder and Boston
lat1 = 40.0150
lon1 = -105.2705
lat2 = 42.3601
lon2 = -71.0589
lat = [
40.015,
40.76997,
41.40982,
41.92961,
42.32515,
42.59314,
42.73126,
42.73831,
42.61423,
42.3601,
]
lon = [
254.7295,
258.3425,
262.0321,
265.7881,
269.5983,
273.4489,
277.3247,
281.2095,
285.087,
288.9411,
]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
# Distance between Nairobi and Cairo (Crosses 0 degrees Latitude)
lat1 = -1.1711
lon1 = 36.4902
lat2 = 30.240
lon2 = 31.149
lat = [
-1.1711,
2.3232,
5.817294,
9.310863,
12.80358,
16.29507,
19.78496,
23.27277,
26.758,
30.24,
]
lon = [
36.4902,
35.95082,
35.40877,
34.85993,
34.3,
33.72431,
33.12772,
32.50438,
31.84749,
31.149,
]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
# Distance between Dublin and London (Crosses 0 degrees Longitude)
lat1 = 53.21
lon1 = -5.1537
lat2 = 51.3026
lon2 = 0.739
lat = [
53.21,
53.01284,
52.81184,
52.60707,
52.39858,
52.18644,
51.97071,
51.75145,
51.52873,
51.3026,
]
lon = [
354.8463,
355.5261,
356.1996,
356.8669,
357.5279,
358.1826,
358.8311,
359.4733,
0.1092852,
0.739,
]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
# Distance between North to South Pole
lat1 = 90
lon1 = 0
lat2 = -90
lon2 = 0
lat = [90, 70, 50, 30, 10, -10, -30, -50, -70, -90]
lon = [
0,
3.508355e-15,
3.508355e-15,
3.508355e-15,
3.508355e-15,
3.508355e-15,
3.508355e-15,
3.508355e-15,
3.508355e-15,
0,
]
ncl_gc_latlon_npts[(lat1, lon1, lat2, lon2)] = list(zip(lat, lon))
from pyproj import Geod
python_gc_latlon_npts = {}
for latlon in ncl_gc_latlon_dist.keys():
geodesic = Geod(ellps="sphere")
pts = geodesic.npts(latlon[1], latlon[0], latlon[3], latlon[2], 8)
python_gc_latlon_npts[latlon] = [[pt[1], pt[0] % 360] for pt in pts]
# add starting and ending points of GC arc
lat1, lon1 = latlon[0], latlon[1]
lat2, lon2 = latlon[2], latlon[3]
python_gc_latlon_npts[latlon] = (
[[lat1, lon1]] + python_gc_latlon_npts[latlon] + [[lat2, lon2]]
)
Comparison#
area_poly_sphere#
Important Note
NCL does not return a valid value for a single point ("Single Point -> Invalid NCL") where Python does, so this error is ignored belowfor key in ncl_results.keys():
if key in python_results.keys() and key in python_results.keys():
try:
assert math.isclose(
ncl_results[key], python_results[key], rel_tol=1e-5
) # within 4 decimal points
print(
f"VALID:\n\t{key}: \n\tncl:\t\t{ncl_results[key]}\n\tpython:\t\t{python_results[key]}\n"
)
except Exception:
print(
f"INVALID:\n\t{key}: \n\t\tncl:\t\t{ncl_results[key]}\n\t\tpython:\t\t{python_results[key]}\n"
)
if key != "Single Point -> Invalid NCL":
assert False
VALID:
Boulder, Boston, Houston:
ncl: 1940856
python: 1940855.984630871
VALID:
Four Corners of Colorado:
ncl: 250007.6
python: 250007.81588628562
VALID:
Caltech, Alberta, Greenwich, Paris, Madrid:
ncl: 11634800
python: 11634798.025760762
VALID:
Crossing the Equator:
ncl: 114894.8
python: 114894.73868013734
VALID:
Crossing the Prime Meridian:
ncl: 54450.39
python: 54450.41832867822
VALID:
Half of the World:
ncl: 255032000
python: 255031995.77390912
INVALID:
Single Point -> Invalid NCL:
ncl: -127516000
python: 0.0
VALID:
Single Degree:
ncl: 9401.705
python: 9401.704877506347
css2c#
for key in ncl_css2c.keys():
if key in astropy_css2c.keys():
assert abs(ncl_css2c[key][0] - astropy_css2c[key][0]) < 0.000005
assert abs(ncl_css2c[key][1] - astropy_css2c[key][1]) < 0.000005
assert abs(ncl_css2c[key][2] - astropy_css2c[key][2]) < 0.000005
csc2s#
Important Note
To generate the Cartesian coordinates to test against, the NCL script for this receipt converts a range of latitude/longitude to Cartesian coordinates (with the `css2c` function). The Carestian coordinates are then converted back into latitude/longitude with the `csc2s` function. This allows the receipt to test `csc2s` across a full range of coordinates. However, NCL coordinates representing the poles (+90/-90) and the antimeridian (+180/-180) produced through this process return as an equivalent, but different value. For example, an input at the pole (-90, -179) produces an output of (-90, 1) and an input of (-90,13) produces an output (-90,-167).ncl 0> cart = css2c(-90, 87)
ncl 1> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))
(0,0) -90
(1,0) -92.99999
The same applies for the antimerdian where, for example, an input of (-89,-180) produces an output of (-89,180)
ncl 4> cart = css2c(89,180)
ncl 5> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))
(0,0) 89.00005
(1,0) -180
# Verify Latitude/Longitude Inputs match the Latitude/Longtiude Outputs
for key in ncl_csc2s_input_output.keys():
try:
assert ncl_csc2s_input_output[key][0] == key[0]
assert ncl_csc2s_input_output[key][1] == key[1]
except Exception:
if (
abs(ncl_csc2s_input_output[key][0]) != 90
and abs(ncl_csc2s_input_output[key][1]) != 180
):
print(Exception)
# Expected places where input lat/lon will not match output lat/lon in NCL
# NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)
# NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)
assert False
# Verify conversions from cartesian coordinates to latitude/longtiude
for i, key in enumerate(ncl_csc2s.keys()):
if i % 3 == 0: # test every third point to prevent CellTimeoutError
try:
assert abs(key[0] - list(astropy_csc2s.keys())[i][0]) < 0.0005
assert abs(key[1] - list(astropy_csc2s.keys())[i][1]) < 0.0005
except Exception:
if not math.isclose(abs(key[0]), 90) and not math.isclose(abs(key[1]), 180):
print(abs(key[0] - list(astropy_csc2s.keys())[i][0]))
print(abs(key[1] - list(astropy_csc2s.keys())[i][1]))
# Expected places where input lat/lon will not match output lat/lon in NCL
# NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)
# NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)
assert False
gc_onarc#
Important Note
The tolerance for the function to check if a point lies on a great circle arc varies between NCL and UXarray which can generate some potential mismatches when comparing between each functionfor latlon in ncl_gc_onarc.keys():
assert ncl_gc_onarc[latlon] is python_onarc[latlon]
gc_qarea#
import math
for latlon in ncl_gc_qarea.keys():
assert math.isclose(
ncl_gc_qarea[latlon], python_gc_qarea[latlon], rel_tol=1e-2
) # within 2 decimal points
gc_tarea#
import math
for latlon in ncl_gc_tarea.keys():
assert math.isclose(
ncl_gc_tarea[latlon], python_gc_tarea[latlon], rel_tol=1e-2
) # within 2 decimal points
gc_latlon#
Great Circle Distance#
import math
for latlon in ncl_gc_latlon_dist.keys():
assert math.isclose(
ncl_gc_latlon_dist[latlon], python_gc_latlon_dist[latlon], rel_tol=1e-2
) # within 2 decimal points
Interpolated Points#
import math
for latlon in ncl_gc_latlon_npts.keys():
for i, pts in enumerate(ncl_gc_latlon_npts[latlon]):
assert (
ncl_gc_latlon_npts[latlon][i][0]
- (python_gc_latlon_npts[latlon][i][0] % 360)
< 0.01
)
assert (
ncl_gc_latlon_npts[latlon][i][1]
- (python_gc_latlon_npts[latlon][i][1] % 360)
< 0.01
)