Skip to content

Commit efc8fd7

Browse files
cyschneckcyschneck
authored
GC: area_poly_sphere (#130)
1 parent b15325b commit efc8fd7

File tree

6 files changed

+397
-0
lines changed

6 files changed

+397
-0
lines changed

‎ncl/ncl_entries/great_circle.ipynb

+111
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Great Circle\n"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"## Overview\n",
15+
"\n",
16+
"This section covers great circle functions from NCL:\n",
17+
"\n",
18+
"- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)"
19+
]
20+
},
21+
{
22+
"cell_type": "markdown",
23+
"metadata": {},
24+
"source": [
25+
"## area_poly_sphere\n",
26+
"NCL's `area_poly_sphere` calculates the area enclosed by an arbitrary polygon on the sphere\n",
27+
"\n",
28+
"<div class=\"admonition alert alert-info\">\n",
29+
" <p class=\"admonition-title\" style=\"font-weight:bold\">Important Note</p>\n",
30+
" Coordinates should be within the valid latitude/longitude range (-90° to 90° and -180° to 180°) and be in clockwise order\n",
31+
"</div>\n",
32+
"\n",
33+
"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{footcite}`usgs_1987` model)"
34+
]
35+
},
36+
{
37+
"cell_type": "markdown",
38+
"metadata": {},
39+
"source": [
40+
"### Grab and Go"
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": null,
46+
"metadata": {},
47+
"outputs": [],
48+
"source": [
49+
"from pyproj import Geod\n",
50+
"\n",
51+
"# Points in clockise order: Boulder, Boston, Houston\n",
52+
"latitudes = [40.0150, 42.3601, 29.5518] # degrees\n",
53+
"longitudes = [-105.2705, -71.0589, -95.0982] # degrees\n",
54+
"\n",
55+
"geod = Geod(ellps=\"sphere\") # radius = 6370997 m\n",
56+
"poly_area_m, _ = geod.polygon_area_perimeter(longitudes, latitudes)\n",
57+
"poly_area_km2 = abs(poly_area_m) * 1e-6\n",
58+
"poly_area_km2"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"metadata": {},
64+
"source": [
65+
"---"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"## Python Resources\n",
73+
"- [pyroj.geod() great circle computations](https://pyproj4.github.io/pyproj/stable/api/geod.html)\n",
74+
"\n",
75+
"## Additional Reading\n",
76+
"- [Aviation Formulary for working with great circles](https://www.edwilliams.org/avform147.htm)"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"## References:\n",
84+
"\n",
85+
"```{footbibliography}\n",
86+
"```"
87+
]
88+
}
89+
],
90+
"metadata": {
91+
"kernelspec": {
92+
"display_name": "Python 3 (ipykernel)",
93+
"language": "python",
94+
"name": "python3"
95+
},
96+
"language_info": {
97+
"codemirror_mode": {
98+
"name": "ipython",
99+
"version": 3
100+
},
101+
"file_extension": ".py",
102+
"mimetype": "text/x-python",
103+
"name": "python",
104+
"nbconvert_exporter": "python",
105+
"pygments_lexer": "ipython3",
106+
"version": "3.12.7"
107+
}
108+
},
109+
"nbformat": 4,
110+
"nbformat_minor": 4
111+
}

‎ncl/ncl_entries/ncl_entries.rst

+1
Original file line numberDiff line numberDiff line change
@@ -38,3 +38,4 @@ and links to additional resources and references.
3838

3939
- `Climatology <climatology.ipynb>`_
4040
- `Meteorology <meteorology.ipynb>`_
41+
- `Great circle <great_circle.ipynb>`_

‎ncl/ncl_index/ncl-index-table.csv

+1
Original file line numberDiff line numberDiff line change
@@ -45,3 +45,4 @@ NCL Function,Description,Python Equivalent,Notes
4545
`round <https://www.ncl.ucar.edu/Document/Functions/Built-in/round.shtml>`__,"Rounds a float or double variable to the nearest whole number","``round()`` or ``numpy.round()``",`example notebook <../ncl_entries/general_applied_math.ipynb#decimalplaces-round>`__
4646
`daylight_fao56 <https://www.ncl.ucar.edu/Document/Functions/Crop/daylight_fao56.shtml>`__," Compute maximum number of daylight hours as described in FAO 56","``geocat.comp.meteorology.max_daylight()``",`example notebook <../ncl_entries/meteorology.ipynb#daylight-fao56>`__
4747
`dewtemp_trh <https://www.ncl.ucar.edu/Document/Functions/Built-in/dewtemp_trh.shtml>`__,"Calculates the dew point temperature given temperature and relative humidity","``geocat.comp.dewtemp()``",`example notebook <../ncl_entries/meteorology.ipynb#dewtemp-trh>`__
48+
`area_poly_sphere <https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml>`__,"Calculates the area enclosed by an arbitrary polygon on the sphere","``pyproj.Geod()`` and ``shapely.geometry.polygon.Polygon()``",`example notebook <../ncl_entries/great_circle.ipynb#area-poly-sphere>`__

‎ncl/ncl_raw/great_circle.ncl

+70
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
; area_poly_sphere
2+
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml
3+
4+
; Boulder, Boston, Houston (3 Sides)
5+
lat0 = (/40.0150, 42.3601, 29.5518/)
6+
lon0 = (/-105.2705, -71.0589, -95.0982/)
7+
print(gc_clkwise(lat0, lon0))
8+
; (0) True
9+
print(area_poly_sphere(lat0, lon0, 6370.997))
10+
; (0) 1940856
11+
12+
; Roughly the Four Corners of Colorado (4 Sides)
13+
lat1 = (/41.00488, 41.00203, 37.00540, 37.00051/)
14+
lon1 = (/-109.05001, -102.05348, -103.04633, -109.04720/)
15+
print(gc_clkwise(lat1, lon1))
16+
; (0) True
17+
print(area_poly_sphere(lat1, lon1, 6370.997))
18+
; (0) 250007.6
19+
20+
; Caltech, Alberta, Greenwich, Paris, Madrid (5 Sides)
21+
lat2 = (/34.1377, 53.9333, 51.4934, 48.8575, 40.4167/)
22+
lon2 = (/-118.1253, -116.5765, 0.0098, 2.3514, -3.7033/)
23+
print(gc_clkwise(lat2, lon2))
24+
; (0) True
25+
print(area_poly_sphere(lat2, lon2, 6370.997))
26+
; (0) 1.16348e+07
27+
28+
;; Edge Cases:
29+
30+
; Crossing the Equator (0 degrees Latitude)
31+
; Libreville, Moanda, San Tome and Principe
32+
lat3 = (/0.4078, -5.9230, 0.1864/)
33+
lon3 = (/9.4403, 12.3636, 6.6131/)
34+
print(gc_clkwise(lat3, lon3))
35+
; (0) True
36+
print(area_poly_sphere(lat3, lon3, 6370.997))
37+
; (0) 114894.8
38+
39+
; Crossing the Prime Meridian (0 Degrees Longitude)
40+
; Dublin, Amsterdam, Fishguard
41+
lat4 = (/53.3498, 52.3676, 51.9939/)
42+
lon4 = (/-6.2603, 4.9041, -4.9760/)
43+
print(gc_clkwise(lat4, lon4))
44+
; (0) True
45+
print(area_poly_sphere(lat4, lon4, 6370.997))
46+
; (0) 54450.39
47+
48+
; Half of the World
49+
lat5 = (/90.0, 0.0, -90.0, 0.0/)
50+
lon5 = (/0.0, -90.0, 0.0, 90.0/)
51+
print(gc_clkwise(lat5, lon5))
52+
; (0) True
53+
print(area_poly_sphere(lat5, lon5, 6370.997))
54+
; (0) 2.55032e+08
55+
56+
; Single Point -> Invalid NCL
57+
lat6 = (/40, 40, 40, 40/)
58+
lon6 = (/-105, -105, -105, -105/)
59+
print(gc_clkwise(lat6, lon6))
60+
; (0) True
61+
print(area_poly_sphere(lat6, lon6, 6370.997))
62+
; (0) -1.27516e+08
63+
64+
; Single Degree
65+
lat7 = (/40, 40, 41, 41/)
66+
lon7 = (/-105, -106, -106, -105/)
67+
print(gc_clkwise(lat7, lon7))
68+
; (0) True
69+
print(area_poly_sphere(lat7, lon7, 6370.997))
70+
;(0) 9401.705

‎ncl/receipts/great_circle.ipynb

+206
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,206 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "e6eddaf181eacfd3",
6+
"metadata": {},
7+
"source": [
8+
"# Great Circle"
9+
]
10+
},
11+
{
12+
"cell_type": "markdown",
13+
"id": "e6be6389ef38b00b",
14+
"metadata": {},
15+
"source": [
16+
"```{warning} This is not meant to be a standalone notebook.\n",
17+
"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.\n",
18+
"```"
19+
]
20+
},
21+
{
22+
"cell_type": "markdown",
23+
"id": "536ffded4355d9c6",
24+
"metadata": {},
25+
"source": [
26+
"## Functions covered\n",
27+
"- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)"
28+
]
29+
},
30+
{
31+
"cell_type": "markdown",
32+
"id": "4a129c971c083695",
33+
"metadata": {},
34+
"source": [
35+
"## NCL code"
36+
]
37+
},
38+
{
39+
"cell_type": "markdown",
40+
"id": "3d70616a8934f0fb",
41+
"metadata": {},
42+
"source": [
43+
"```{literalinclude} ../ncl_raw/great_circle.ncl\n",
44+
"\n",
45+
"```"
46+
]
47+
},
48+
{
49+
"cell_type": "markdown",
50+
"id": "d918dec004b6456b",
51+
"metadata": {},
52+
"source": [
53+
"## Python Functionality"
54+
]
55+
},
56+
{
57+
"cell_type": "code",
58+
"execution_count": null,
59+
"id": "95225429-62d5-4d38-b170-850526c28107",
60+
"metadata": {},
61+
"outputs": [],
62+
"source": [
63+
"import math\n",
64+
"\n",
65+
"ncl_results = {}\n",
66+
"\n",
67+
"poly_name = [\n",
68+
" \"Boulder, Boston, Houston\",\n",
69+
" \"Four Corners of Colorado\",\n",
70+
" \"Caltech, Alberta, Greenwich, Paris, Madrid\",\n",
71+
" \"Crossing the Equator\",\n",
72+
" \"Crossing the Prime Meridian\",\n",
73+
" \"Half of the World\",\n",
74+
" \"Single Point -> Invalid NCL\",\n",
75+
" \"Single Degree\",\n",
76+
"]"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"id": "6dc1d6c4-ac26-4346-9aa3-2c73182a7511",
82+
"metadata": {},
83+
"source": [
84+
"### area_poly_sphere"
85+
]
86+
},
87+
{
88+
"cell_type": "code",
89+
"execution_count": null,
90+
"id": "2ba99b37-a743-4776-bc7b-d5a08b977642",
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"ncl_lat = [\n",
95+
" [40.0150, 42.3601, 29.5518],\n",
96+
" [41.00488, 41.00203, 37.00540, 37.00051],\n",
97+
" [34.1377, 53.9333, 51.4934, 48.8575, 40.4167],\n",
98+
" [0.4078, -5.9230, 0.1864],\n",
99+
" [53.3498, 52.3676, 51.9939],\n",
100+
" [90.0, 0.0, -90.0, 0.0],\n",
101+
" [40, 40, 40, 40],\n",
102+
" [40, 40, 41, 41],\n",
103+
"]\n",
104+
"ncl_lon = [\n",
105+
" [-105.2705, -71.0589, -95.0982],\n",
106+
" [-109.05001, -102.05348, -103.04633, -109.04720],\n",
107+
" [-118.1253, -116.5765, 0.0098, 2.3514, -3.7033],\n",
108+
" [9.4403, 12.3636, 6.6131],\n",
109+
" [-6.2603, 4.9041, -4.9760],\n",
110+
" [0.0, -90.0, 0.0, 90.0],\n",
111+
" [-105, -105, -105, -105],\n",
112+
" [-105, -106, -106, -105],\n",
113+
"]\n",
114+
"ncl_results[poly_name[0]] = 1940856\n",
115+
"ncl_results[poly_name[1]] = 250007.6\n",
116+
"ncl_results[poly_name[2]] = 11634800\n",
117+
"ncl_results[poly_name[3]] = 114894.8\n",
118+
"ncl_results[poly_name[4]] = 54450.39\n",
119+
"ncl_results[poly_name[5]] = 255032000\n",
120+
"ncl_results[poly_name[6]] = -127516000\n",
121+
"ncl_results[poly_name[7]] = 9401.705"
122+
]
123+
},
124+
{
125+
"cell_type": "code",
126+
"execution_count": null,
127+
"id": "35803c5d-bf83-4a35-b61c-50466d9d5095",
128+
"metadata": {},
129+
"outputs": [],
130+
"source": [
131+
"from pyproj import Geod\n",
132+
"\n",
133+
"python_results = {}\n",
134+
"geod = Geod(ellps=\"sphere\") # radius = 6370997 m\n",
135+
"for i in range(len(poly_name)):\n",
136+
" poly_area_m, _ = geod.polygon_area_perimeter(ncl_lon[i], ncl_lat[i])\n",
137+
" poly_area_km2 = abs(poly_area_m) * 1e-6\n",
138+
" python_results[poly_name[i]] = poly_area_km2"
139+
]
140+
},
141+
{
142+
"cell_type": "markdown",
143+
"id": "3237a0bffc6827fc",
144+
"metadata": {},
145+
"source": [
146+
"## Comparison"
147+
]
148+
},
149+
{
150+
"cell_type": "markdown",
151+
"id": "cb61d4f0-bd10-4034-bfe8-b0c6e9137ec3",
152+
"metadata": {},
153+
"source": [
154+
"<div class=\"admonition alert alert-info\">\n",
155+
" <p class=\"admonition-title\" style=\"font-weight:bold\">Important Note</p>\n",
156+
" NCL does not return a valid value for a single point (\"Single Point -> Invalid NCL\") where Python does, so this error is ignored below\n",
157+
"</div>"
158+
]
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": null,
163+
"id": "74362fd9-0e9f-4cf9-91da-08cd81be625c",
164+
"metadata": {},
165+
"outputs": [],
166+
"source": [
167+
"for key in ncl_results.keys():\n",
168+
" if key in python_results.keys() and key in python_results.keys():\n",
169+
" try:\n",
170+
" assert math.isclose(\n",
171+
" ncl_results[key], python_results[key], rel_tol=1e-5\n",
172+
" ) # within 4 decimal points\n",
173+
" print(\n",
174+
" f\"VALID:\\n\\t{key}: \\n\\tncl:\\t\\t{ncl_results[key]}\\n\\tpython:\\t\\t{python_results[key]}\\n\"\n",
175+
" )\n",
176+
" except Exception:\n",
177+
" print(\n",
178+
" f\"INVALID:\\n\\t{key}: \\n\\t\\tncl:\\t\\t{ncl_results[key]}\\n\\t\\tpython:\\t\\t{python_results[key]}\\n\"\n",
179+
" )\n",
180+
" if key != \"Single Point -> Invalid NCL\":\n",
181+
" assert False"
182+
]
183+
}
184+
],
185+
"metadata": {
186+
"kernelspec": {
187+
"display_name": "Python 3 (ipykernel)",
188+
"language": "python",
189+
"name": "python3"
190+
},
191+
"language_info": {
192+
"codemirror_mode": {
193+
"name": "ipython",
194+
"version": 3
195+
},
196+
"file_extension": ".py",
197+
"mimetype": "text/x-python",
198+
"name": "python",
199+
"nbconvert_exporter": "python",
200+
"pygments_lexer": "ipython3",
201+
"version": "3.12.7"
202+
}
203+
},
204+
"nbformat": 4,
205+
"nbformat_minor": 5
206+
}

0 commit comments

Comments
 (0)