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
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 = 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 gdf
from astropy.coordinates.representation import (
CartesianRepresentation,
SphericalRepresentation,
)
import numpy as np
csc2s_data = gdf.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)
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