GC: csc2s

ncl/ncl_entries/great_circle.ipynb

@@ -16,7 +16,8 @@
     "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)\n",
-    "- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)"
+    "- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)\n",
+    "- [csc2s](https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml)"
    ]
   },
   {
@@ -93,6 +94,50 @@
     "print(f\"Z = {cart_coords.z.value}\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## csc2s\n",
+    "NCL's `csc2s` converts Cartesian coordinates to spherical (latitude/longitude) 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 (\n",
+    "    CartesianRepresentation,\n",
+    "    SphericalRepresentation,\n",
+    ")\n",
+    "import numpy as np\n",
+    "\n",
+    "x = -0.20171369272651396\n",
+    "y = -0.7388354627678497\n",
+    "z = 0.6429881376224998\n",
+    "\n",
+    "cart_coords = CartesianRepresentation(x=x, y=y, z=z)\n",
+    "spherical_coords = cart_coords.represent_as(SphericalRepresentation)\n",
+    "\n",
+    "# convert latitude/longitude from radians to degrees\n",
+    "lat_deg = np.rad2deg(spherical_coords.lat.value)\n",
+    "lon_deg = (\n",
+    "    np.rad2deg(spherical_coords.lon.value) + 180\n",
+    ") % 360 - 180  # keep longitude between -180 to 180\n",
+    "\n",
+    "print(f\"Latitude = {lat_deg}\")\n",
+    "print(f\"Longitude = {lon_deg}\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

ncl/ncl_index/ncl-index-table.csv

@@ -51,3 +51,4 @@ NCL Function,Description,Python Equivalent,Notes
 `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>`__
 `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>`__
 `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>`__
+`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

@@ -83,3 +83,23 @@ do lat=-90,90
 		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

ncl/receipts/great_circle.ipynb

@@ -25,7 +25,8 @@
    "source": [
     "## Functions covered\n",
     "- [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)"
+    "- [css2c](https://www.ncl.ucar.edu/Document/Functions/Built-in/css2c.shtml)\n",
+    "- [csc2s](https://www.ncl.ucar.edu/Document/Functions/Built-in/csc2s.shtml)"
    ]
   },
   {
@@ -147,7 +148,7 @@
     "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",
+    "cart_values = tuple(css2c_data[::, 2:])\n",
     "ncl_css2c = dict(zip(lat_lon, cart_values))"
    ]
   },
@@ -175,6 +176,64 @@
     "    astropy_css2c[pair] = cart_coords.xyz.value"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "399d047d-f22c-41cf-996d-d84e1f5b096c",
+   "metadata": {},
+   "source": [
+    "### csc2s"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ae18411-506d-455c-867e-50273bfff7e2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geocat.datafiles as gdf\n",
+    "from astropy.coordinates.representation import (\n",
+    "    CartesianRepresentation,\n",
+    "    SphericalRepresentation,\n",
+    ")\n",
+    "import numpy as np\n",
+    "\n",
+    "csc2s_data = gdf.get('applications_files/ncl_outputs/csc2s_output.txt')\n",
+    "csc2s_data = np.loadtxt(csc2s_data, delimiter=',', skiprows=6)\n",
+    "\n",
+    "input_lat_lon = tuple(map(tuple, csc2s_data[::, 0:2]))\n",
+    "cart_values = tuple(map(tuple, (csc2s_data[::, 2:5])))\n",
+    "output_lat_lon = tuple(map(tuple, (csc2s_data[::, 5:])))\n",
+    "ncl_csc2s = dict(zip(input_lat_lon, cart_values))\n",
+    "ncl_csc2s_input_output = dict(zip(input_lat_lon, output_lat_lon))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "043a0ec3-e0a2-4c6e-952d-1d459237a1f4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Calculate Spherical coordinates\n",
+    "def spherical_cart(x, y, z):\n",
+    "    cart_coords = CartesianRepresentation(x=x, y=y, z=z)\n",
+    "    spherical_coords = cart_coords.represent_as(SphericalRepresentation)\n",
+    "    # convert latitude/longitude from radians to degrees\n",
+    "    lat_deg = np.rad2deg(spherical_coords.lat.value)\n",
+    "    lon_deg = (\n",
+    "        np.rad2deg(spherical_coords.lon.value) + 180\n",
+    "    ) % 360 - 180  # keep longitude between -180 to 180\n",
+    "    return (lat_deg, lon_deg)\n",
+    "\n",
+    "\n",
+    "astropy_csc2s = {}\n",
+    "for xyz in cart_values:\n",
+    "    x, y, z = xyz\n",
+    "    lat_lon = spherical_cart(x, y, z)\n",
+    "    astropy_csc2s[lat_lon] = tuple(xyz)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "3237a0bffc6827fc",
@@ -247,6 +306,87 @@
     "        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"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "90d7474d-60fb-4c22-8191-a120560174af",
+   "metadata": {},
+   "source": [
+    "### csc2s"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa9fb6d4-550b-4d61-85df-51b268a96256",
+   "metadata": {},
+   "source": [
+    "<div class=\"admonition alert alert-info\">\n",
+    "    <p class=\"admonition-title\" style=\"font-weight:bold\">Important Note</p>\n",
+    "    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",
+    "    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",
+    "\n",
+    "```\n",
+    "ncl 0> cart = css2c(-90, 87)\n",
+    "ncl 1> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))\n",
+    "(0,0)\t-90\n",
+    "(1,0)\t-92.99999\n",
+    "```\n",
+    "The same applies for the antimerdian where, for example, an input of (-89,-180) produces an output of (-89,180)\n",
+    "```\n",
+    "ncl 4> cart = css2c(89,180)                          \n",
+    "ncl 5> print(csc2s(cart(0,0), cart(1,0), cart(2,0)))\n",
+    "(0,0)\t89.00005\n",
+    "(1,0)\t-180\n",
+    "```\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2e651572-aeb0-458f-9590-2ee6d2008235",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Verify Latitude/Longitude Inputs match the Latitude/Longtiude Outputs\n",
+    "for key in ncl_csc2s_input_output.keys():\n",
+    "    try:\n",
+    "        assert ncl_csc2s_input_output[key][0] == key[0]\n",
+    "        assert ncl_csc2s_input_output[key][1] == key[1]\n",
+    "    except Exception:\n",
+    "        if (\n",
+    "            abs(ncl_csc2s_input_output[key][0]) != 90\n",
+    "            and abs(ncl_csc2s_input_output[key][1]) != 180\n",
+    "        ):\n",
+    "            print(Exception)\n",
+    "            # Expected places where input lat/lon will not match output lat/lon in NCL\n",
+    "            # NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)\n",
+    "            # NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)\n",
+    "            assert False"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "53360a09-f1f3-4b0f-b7b4-9501aa5c92e1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Verify conversions from cartesian coordinates to latitude/longtiude\n",
+    "for i, key in enumerate(ncl_csc2s.keys()):\n",
+    "    if i % 3 == 0:  # test every third point to prevent CellTimeoutError\n",
+    "        try:\n",
+    "            assert abs(key[0] - list(astropy_csc2s.keys())[i][0]) < 0.0005\n",
+    "            assert abs(key[1] - list(astropy_csc2s.keys())[i][1]) < 0.0005\n",
+    "        except Exception:\n",
+    "            if not math.isclose(abs(key[0]), 90) and not math.isclose(abs(key[1]), 180):\n",
+    "                print(abs(key[0] - list(astropy_csc2s.keys())[i][0]))\n",
+    "                print(abs(key[1] - list(astropy_csc2s.keys())[i][1]))\n",
+    "                # Expected places where input lat/lon will not match output lat/lon in NCL\n",
+    "                # NCL produces flipped longitude value for +/-90 latitude, example: (90,-179)->(90,1)\n",
+    "                # NCL produces flipped longitude value for all latitude values when longitude is 180, example: (79,-180)->(79,180)\n",
+    "                assert False"
+   ]
   }
  ],
  "metadata": {
@@ -265,7 +405,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.7"
+   "version": "3.12.8"
   }
  },
  "nbformat": 4,