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
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 gdf
import numpy as np
from astropy.coordinates.representation import UnitSphericalRepresentation
from astropy import units
css2c_data = gdf.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 = 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
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