GC: css2c (#152)

--------- Co-authored-by: Julia Kent <46687291+jukent@users.noreply.github.com>
NCAR · Dec 10, 2024 · 12653be · 12653be
1 parent 796b294
commit 12653be
Show file tree

Hide file tree

Showing 5 changed files with 151 additions and 31 deletions.
diff --git a/environment.yml b/environment.yml
@@ -3,6 +3,7 @@ channels:
   - conda-forge
 dependencies:
   # Content
+  - astropy
   - cartopy
   - geocat-comp
   - geocat-datafiles

diff --git a/ncl/ncl_entries/great_circle.ipynb b/ncl/ncl_entries/great_circle.ipynb
@@ -15,15 +15,16 @@
     "\n",
     "This section covers great circle functions from NCL:\n",
     "\n",
-    "- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)"
+    "- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)\n",
+    "- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## area_poly_sphere\n",
-    "NCL's `area_poly_sphere`  calculates the area enclosed by an arbitrary polygon on the sphere\n",
+    "NCL's `area_poly_sphere` calculates the area enclosed by an arbitrary polygon on the sphere\n",
     "\n",
     "<div class=\"admonition alert alert-info\">\n",
     "    <p class=\"admonition-title\" style=\"font-weight:bold\">Important Note</p>\n",
@@ -58,6 +59,40 @@
     "poly_area_km2"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## css2c\n",
+    "NCL's `css2c` converts spherical (latitude/longitude) coordinates to Cartesian coordinates on a unit sphere"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Grab and Go"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from astropy.coordinates.representation import UnitSphericalRepresentation\n",
+    "from astropy import units\n",
+    "\n",
+    "lat = 40.0150\n",
+    "lon = -105.2705\n",
+    "\n",
+    "spherical_coords = UnitSphericalRepresentation(lat=lat * units.deg, lon=lon * units.deg)\n",
+    "cart_coords = spherical_coords.to_cartesian()\n",
+    "print(f\"X = {cart_coords.x.value}\")\n",
+    "print(f\"Y = {cart_coords.y.value}\")\n",
+    "print(f\"Z = {cart_coords.z.value}\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -71,6 +106,7 @@
    "source": [
     "## Python Resources\n",
     "- [pyroj.geod() great circle computations](https://pyproj4.github.io/pyproj/stable/api/geod.html)\n",
+    "- [Astropy Coordinate Systems](https://docs.astropy.org/en/stable/coordinates/representations.html)\n",
     "\n",
     "## Additional Reading\n",
     "- [Aviation Formulary for working with great circles](https://www.edwilliams.org/avform147.htm)"

diff --git a/ncl/ncl_index/ncl-index-table.csv b/ncl/ncl_index/ncl-index-table.csv
@@ -46,4 +46,5 @@ NCL Function,Description,Python Equivalent,Notes
 `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>`__
 `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>`__
 `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>`__
+`css2c <https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml>`__,"Converts spherical coordinates (lat/lon) to Cartesian coordinates on a unit sphere","``astropy.UnitSphericalRepresentation()``",`example notebook <../ncl_entries/great_circle.ipynb#css2c>`__
 `satvpr_temp_fao56 <https://www.ncl.ucar.edu/Document/Functions/Crop/satvpr_temp_fao56.shtml>`__," Compute saturation vapor pressure using temperature as described in FAO 56","``geocat.comp.saturation_vapor_pressure()``",`example notebook <../ncl_entries/meteorology.ipynb#satvpr-temp-fao56>`__
diff --git a/ncl/ncl_raw/great_circle.ncl b/ncl/ncl_raw/great_circle.ncl
@@ -68,3 +68,18 @@ 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
diff --git a/ncl/receipts/great_circle.ipynb b/ncl/receipts/great_circle.ipynb
@@ -24,7 +24,8 @@
    "metadata": {},
    "source": [
     "## Functions covered\n",
-    "- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)"
+    "- [area_poly_sphere](https://www.ncl.ucar.edu/Document/Functions/Built-in/area_poly_sphere.shtml)\n",
+    "- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)"
    ]
   },
   {
@@ -53,10 +54,18 @@
     "## Python Functionality"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "6dc1d6c4-ac26-4346-9aa3-2c73182a7511",
+   "metadata": {},
+   "source": [
+    "### area_poly_sphere"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "95225429-62d5-4d38-b170-850526c28107",
+   "id": "2ba99b37-a743-4776-bc7b-d5a08b977642",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -73,24 +82,8 @@
     "    \"Half of the World\",\n",
     "    \"Single Point -> Invalid NCL\",\n",
     "    \"Single Degree\",\n",
-    "]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "6dc1d6c4-ac26-4346-9aa3-2c73182a7511",
-   "metadata": {},
-   "source": [
-    "### area_poly_sphere"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "2ba99b37-a743-4776-bc7b-d5a08b977642",
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "]\n",
+    "\n",
     "ncl_lat = [\n",
     "    [40.0150, 42.3601, 29.5518],\n",
     "    [41.00488, 41.00203, 37.00540, 37.00051],\n",
@@ -118,24 +111,68 @@
     "ncl_results[poly_name[4]] = 54450.39\n",
     "ncl_results[poly_name[5]] = 255032000\n",
     "ncl_results[poly_name[6]] = -127516000\n",
-    "ncl_results[poly_name[7]] = 9401.705"
+    "ncl_results[poly_name[7]] = 9401.705\n",
+    "\n",
+    "from pyproj import Geod\n",
+    "\n",
+    "python_results = {}\n",
+    "geod = Geod(ellps=\"sphere\")  # radius = 6370997 m\n",
+    "for i in range(len(poly_name)):\n",
+    "    poly_area_m, _ = geod.polygon_area_perimeter(ncl_lon[i], ncl_lat[i])\n",
+    "    poly_area_km2 = abs(poly_area_m) * 1e-6\n",
+    "    python_results[poly_name[i]] = poly_area_km2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c0e65c3e-2377-47ec-b94c-e7eb753966d9",
+   "metadata": {},
+   "source": [
+    "### css2c"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "35803c5d-bf83-4a35-b61c-50466d9d5095",
+   "id": "35abde81-5843-4504-8e32-a137ee1aa094",
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pyproj import Geod\n",
+    "import geocat.datafiles as gdf\n",
+    "import numpy as np\n",
+    "from astropy.coordinates.representation import UnitSphericalRepresentation\n",
+    "from astropy import units\n",
     "\n",
-    "python_results = {}\n",
-    "geod = Geod(ellps=\"sphere\")  # radius = 6370997 m\n",
-    "for i in range(len(poly_name)):\n",
-    "    poly_area_m, _ = geod.polygon_area_perimeter(ncl_lon[i], ncl_lat[i])\n",
-    "    poly_area_km2 = abs(poly_area_m) * 1e-6\n",
-    "    python_results[poly_name[i]] = poly_area_km2"
+    "css2c_data = gdf.get('applications_files/ncl_outputs/css2c_output.txt')\n",
+    "css2c_data = np.loadtxt(css2c_data, delimiter=',', skiprows=6)\n",
+    "\n",
+    "lat_lon = tuple(map(tuple, (css2c_data[::, 0:2])))\n",
+    "cart_values = css2c_data[::, 2:]\n",
+    "ncl_css2c = dict(zip(lat_lon, cart_values))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee104ea1-e287-4635-b404-5b06ccfb6949",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Collect Latitude and Longtiude Pairs\n",
+    "lat_lon = []\n",
+    "for lat in range(-90, 90 + 1):\n",
+    "    for lon in range(-180, 180 + 1):\n",
+    "        lat_lon.append((lat, lon))\n",
+    "\n",
+    "## Calculate Cartesian coordinates\n",
+    "astropy_css2c = {}\n",
+    "for i, pair in enumerate(lat_lon):\n",
+    "    lat, lon = pair\n",
+    "    spherical_coords = UnitSphericalRepresentation(\n",
+    "        lat=lat * units.deg, lon=lon * units.deg\n",
+    "    )\n",
+    "    cart_coords = spherical_coords.to_cartesian()\n",
+    "    astropy_css2c[pair] = cart_coords.xyz.value"
    ]
   },
   {
@@ -146,6 +183,14 @@
     "## Comparison"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "384790de-53f9-45d7-8bc1-fef86df21b57",
+   "metadata": {},
+   "source": [
+    "### area_poly_sphere"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "cb61d4f0-bd10-4034-bfe8-b0c6e9137ec3",
@@ -180,6 +225,28 @@
     "            if key != \"Single Point -> Invalid NCL\":\n",
     "                assert False"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "57251c09-ef8c-438f-b3cf-cf798e7f4028",
+   "metadata": {},
+   "source": [
+    "### css2c"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d5738d2d-5abe-4ed4-8e90-c42881c2bfb0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for key in ncl_css2c.keys():\n",
+    "    if key in astropy_css2c.keys():\n",
+    "        assert abs(ncl_css2c[key][0] - astropy_css2c[key][0]) < 0.000005\n",
+    "        assert abs(ncl_css2c[key][1] - astropy_css2c[key][1]) < 0.000005\n",
+    "        assert abs(ncl_css2c[key][2] - astropy_css2c[key][2]) < 0.000005"
+   ]
   }
  ],
  "metadata": {