Skip to content

Commit b7b4fbf

Browse files
cyschneckcyschneck
authored
GC: csc2s (#162)
1 parent 3008321 commit b7b4fbf

File tree

4 files changed

+210
-4
lines changed

4 files changed

+210
-4
lines changed

‎ncl/ncl_entries/great_circle.ipynb

+46-1
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,8 @@
1616
"This section covers great circle functions from NCL:\n",
1717
"\n",
1818
"- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)\n",
19-
"- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)"
19+
"- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)\n",
20+
"- [csc2s](https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml)"
2021
]
2122
},
2223
{
@@ -93,6 +94,50 @@
9394
"print(f\"Z = {cart_coords.z.value}\")"
9495
]
9596
},
97+
{
98+
"cell_type": "markdown",
99+
"metadata": {},
100+
"source": [
101+
"## csc2s\n",
102+
"NCL's `csc2s` converts Cartesian coordinates to spherical (latitude/longitude) coordinates on a unit sphere"
103+
]
104+
},
105+
{
106+
"cell_type": "markdown",
107+
"metadata": {},
108+
"source": [
109+
"### Grab and Go"
110+
]
111+
},
112+
{
113+
"cell_type": "code",
114+
"execution_count": null,
115+
"metadata": {},
116+
"outputs": [],
117+
"source": [
118+
"from astropy.coordinates.representation import (\n",
119+
" CartesianRepresentation,\n",
120+
" SphericalRepresentation,\n",
121+
")\n",
122+
"import numpy as np\n",
123+
"\n",
124+
"x = -0.20171369272651396\n",
125+
"y = -0.7388354627678497\n",
126+
"z = 0.6429881376224998\n",
127+
"\n",
128+
"cart_coords = CartesianRepresentation(x=x, y=y, z=z)\n",
129+
"spherical_coords = cart_coords.represent_as(SphericalRepresentation)\n",
130+
"\n",
131+
"# convert latitude/longitude from radians to degrees\n",
132+
"lat_deg = np.rad2deg(spherical_coords.lat.value)\n",
133+
"lon_deg = (\n",
134+
" np.rad2deg(spherical_coords.lon.value) + 180\n",
135+
") % 360 - 180 # keep longitude between -180 to 180\n",
136+
"\n",
137+
"print(f\"Latitude = {lat_deg}\")\n",
138+
"print(f\"Longitude = {lon_deg}\")"
139+
]
140+
},
96141
{
97142
"cell_type": "markdown",
98143
"metadata": {},

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

+1
Original file line numberDiff line numberDiff line change
@@ -51,3 +51,4 @@ NCL Function,Description,Python Equivalent,Notes
5151
`satvpr_tdew_fao56 <https://www.ncl.ucar.edu/Document/Functions/Crop/satvpr_tdew_fao56.shtml>`__,"Compute actual saturation vapor pressure as described in FAO 56","``geocat.comp.actual_saturation_vapor_pressure()``",`example notebook <../ncl_entries/meteorology.ipynb#satvpr-tdew-fao56>`__
5252
`satvpr_slope_fao56 <https://www.ncl.ucar.edu/Document/Functions/Crop/satvpr_slope_fao56.shtml>`__," Compute the slope of the saturation vapor pressure curve as described in FAO 56","``geocat.comp.saturation_vapor_pressure_slope()``",`example notebook <../ncl_entries/meteorology.ipynb#satvpr-slope-fao56>`__
5353
`coriolis_param <https://www.ncl.ucar.edu/Document/Functions/Contributed/coriolis_param.shtml>`__,"Calculate the Coriolis parameter","``metpy.calc.coriolis_parameter()``",`example notebook <../ncl_entries/meteorology.ipynb#coriolis-param>`__
54+
`csc2s <https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml>`__,"Converts Cartesian coordinates on a unit sphere to spherical coordinates (lat/lon)","``astropy.coordinates.representation``",`example notebook <../ncl_entries/great_circle.ipynb#csc2s>`__

‎ncl/ncl_raw/great_circle.ncl

+20
Original file line numberDiff line numberDiff line change
@@ -83,3 +83,23 @@ do lat=-90,90
8383
end
8484
end do
8585
end do
86+
87+
; csc2s
88+
; Adapted from https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml
89+
90+
; ncl -n csc2s.ncl >> csc2s_output.txt
91+
92+
print("Input Latitude (Degree), Input Longitude (Degree), Cartesian X, Cartesian Y, Cartesian Z, Output Latitude (Degree), Output Longitude (Degree)")
93+
do lat=-90,90
94+
do lon=-180,180
95+
begin
96+
cart = css2c(lat, lon)
97+
; determine a list of xyz coordinates based on lat/lon
98+
x = cart(0,0)
99+
y = cart(1,0)
100+
z = cart(2,0)
101+
sph = csc2s(x, y, z)
102+
print(lat + "," + lon + "," + x + "," + y + "," + z + "," + sph(0,0) + "," + sph(1,0))
103+
end
104+
end do
105+
end do

‎ncl/receipts/great_circle.ipynb

+143-3
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,8 @@
2525
"source": [
2626
"## Functions covered\n",
2727
"- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)\n",
28-
"- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)"
28+
"- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)\n",
29+
"- [csc2s](https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml)"
2930
]
3031
},
3132
{
@@ -147,7 +148,7 @@
147148
"css2c_data = np.loadtxt(css2c_data, delimiter=',', skiprows=6)\n",
148149
"\n",
149150
"lat_lon = tuple(map(tuple, (css2c_data[::, 0:2])))\n",
150-
"cart_values = css2c_data[::, 2:]\n",
151+
"cart_values = tuple(css2c_data[::, 2:])\n",
151152
"ncl_css2c = dict(zip(lat_lon, cart_values))"
152153
]
153154
},
@@ -175,6 +176,64 @@
175176
" astropy_css2c[pair] = cart_coords.xyz.value"
176177
]
177178
},
179+
{
180+
"cell_type": "markdown",
181+
"id": "399d047d-f22c-41cf-996d-d84e1f5b096c",
182+
"metadata": {},
183+
"source": [
184+
"### csc2s"
185+
]
186+
},
187+
{
188+
"cell_type": "code",
189+
"execution_count": null,
190+
"id": "5ae18411-506d-455c-867e-50273bfff7e2",
191+
"metadata": {},
192+
"outputs": [],
193+
"source": [
194+
"import geocat.datafiles as gdf\n",
195+
"from astropy.coordinates.representation import (\n",
196+
" CartesianRepresentation,\n",
197+
" SphericalRepresentation,\n",
198+
")\n",
199+
"import numpy as np\n",
200+
"\n",
201+
"csc2s_data = gdf.get('applications_files/ncl_outputs/csc2s_output.txt')\n",
202+
"csc2s_data = np.loadtxt(csc2s_data, delimiter=',', skiprows=6)\n",
203+
"\n",
204+
"input_lat_lon = tuple(map(tuple, csc2s_data[::, 0:2]))\n",
205+
"cart_values = tuple(map(tuple, (csc2s_data[::, 2:5])))\n",
206+
"output_lat_lon = tuple(map(tuple, (csc2s_data[::, 5:])))\n",
207+
"ncl_csc2s = dict(zip(input_lat_lon, cart_values))\n",
208+
"ncl_csc2s_input_output = dict(zip(input_lat_lon, output_lat_lon))"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"id": "043a0ec3-e0a2-4c6e-952d-1d459237a1f4",
215+
"metadata": {},
216+
"outputs": [],
217+
"source": [
218+
"## Calculate Spherical coordinates\n",
219+
"def spherical_cart(x, y, z):\n",
220+
" cart_coords = CartesianRepresentation(x=x, y=y, z=z)\n",
221+
" spherical_coords = cart_coords.represent_as(SphericalRepresentation)\n",
222+
" # convert latitude/longitude from radians to degrees\n",
223+
" lat_deg = np.rad2deg(spherical_coords.lat.value)\n",
224+
" lon_deg = (\n",
225+
" np.rad2deg(spherical_coords.lon.value) + 180\n",
226+
" ) % 360 - 180 # keep longitude between -180 to 180\n",
227+
" return (lat_deg, lon_deg)\n",
228+
"\n",
229+
"\n",
230+
"astropy_csc2s = {}\n",
231+
"for xyz in cart_values:\n",
232+
" x, y, z = xyz\n",
233+
" lat_lon = spherical_cart(x, y, z)\n",
234+
" astropy_csc2s[lat_lon] = tuple(xyz)"
235+
]
236+
},
178237
{
179238
"cell_type": "markdown",
180239
"id": "3237a0bffc6827fc",
@@ -247,6 +306,87 @@
247306
" assert abs(ncl_css2c[key][1] - astropy_css2c[key][1]) < 0.000005\n",
248307
" assert abs(ncl_css2c[key][2] - astropy_css2c[key][2]) < 0.000005"
249308
]
309+
},
310+
{
311+
"cell_type": "markdown",
312+
"id": "90d7474d-60fb-4c22-8191-a120560174af",
313+
"metadata": {},
314+
"source": [
315+
"### csc2s"
316+
]
317+
},
318+
{
319+
"cell_type": "markdown",
320+
"id": "fa9fb6d4-550b-4d61-85df-51b268a96256",
321+
"metadata": {},
322+
"source": [
323+
"<div class=\"admonition alert alert-info\">\n",
324+
" <p class=\"admonition-title\" style=\"font-weight:bold\">Important Note</p>\n",
325+
" 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. \n",
326+
" 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).\n",
327+
"\n",
328+
"```\n",
329+
"ncl 0> cart = css2c(-90, 87)\n",
330+
"ncl 1> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))\n",
331+
"(0,0)\t-90\n",
332+
"(1,0)\t-92.99999\n",
333+
"```\n",
334+
"The same applies for the antimerdian where, for example, an input of (-89,-180) produces an output of (-89,180)\n",
335+
"```\n",
336+
"ncl 4> cart = css2c(89,180) \n",
337+
"ncl 5> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))\n",
338+
"(0,0)\t89.00005\n",
339+
"(1,0)\t-180\n",
340+
"```\n",
341+
"</div>"
342+
]
343+
},
344+
{
345+
"cell_type": "code",
346+
"execution_count": null,
347+
"id": "2e651572-aeb0-458f-9590-2ee6d2008235",
348+
"metadata": {},
349+
"outputs": [],
350+
"source": [
351+
"# Verify Latitude/Longitude Inputs match the Latitude/Longtiude Outputs\n",
352+
"for key in ncl_csc2s_input_output.keys():\n",
353+
" try:\n",
354+
" assert ncl_csc2s_input_output[key][0] == key[0]\n",
355+
" assert ncl_csc2s_input_output[key][1] == key[1]\n",
356+
" except Exception:\n",
357+
" if (\n",
358+
" abs(ncl_csc2s_input_output[key][0]) != 90\n",
359+
" and abs(ncl_csc2s_input_output[key][1]) != 180\n",
360+
" ):\n",
361+
" print(Exception)\n",
362+
" # Expected places where input lat/lon will not match output lat/lon in NCL\n",
363+
" # NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)\n",
364+
" # NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)\n",
365+
" assert False"
366+
]
367+
},
368+
{
369+
"cell_type": "code",
370+
"execution_count": null,
371+
"id": "53360a09-f1f3-4b0f-b7b4-9501aa5c92e1",
372+
"metadata": {},
373+
"outputs": [],
374+
"source": [
375+
"# Verify conversions from cartesian coordinates to latitude/longtiude\n",
376+
"for i, key in enumerate(ncl_csc2s.keys()):\n",
377+
" if i % 3 == 0: # test every third point to prevent CellTimeoutError\n",
378+
" try:\n",
379+
" assert abs(key[0] - list(astropy_csc2s.keys())[i][0]) < 0.0005\n",
380+
" assert abs(key[1] - list(astropy_csc2s.keys())[i][1]) < 0.0005\n",
381+
" except Exception:\n",
382+
" if not math.isclose(abs(key[0]), 90) and not math.isclose(abs(key[1]), 180):\n",
383+
" print(abs(key[0] - list(astropy_csc2s.keys())[i][0]))\n",
384+
" print(abs(key[1] - list(astropy_csc2s.keys())[i][1]))\n",
385+
" # Expected places where input lat/lon will not match output lat/lon in NCL\n",
386+
" # NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)\n",
387+
" # NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)\n",
388+
" assert False"
389+
]
250390
}
251391
],
252392
"metadata": {
@@ -265,7 +405,7 @@
265405
"name": "python",
266406
"nbconvert_exporter": "python",
267407
"pygments_lexer": "ipython3",
268-
"version": "3.12.7"
408+
"version": "3.12.8"
269409
}
270410
},
271411
"nbformat": 4,

0 commit comments

Comments
 (0)