GC: css2c

environment.yml

@@ -3,6 +3,7 @@ channels:
   - conda-forge
 dependencies:
   # Content
+  - astropy
   - cartopy
   - geocat-comp
   - geocat-datafiles

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,50 @@
     "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": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "X = -0.20171369272651396\n",
+      "Y = -0.7388354627678497\n",
+      "Z = 0.6429881376224998\n"
+     ]
+    }
+   ],
+   "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": {},

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)"
    ]
   },
   {
@@ -118,24 +119,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 +191,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 +233,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": {