diff --git a/ncl/ncl_entries/great_circle.ipynb b/ncl/ncl_entries/great_circle.ipynb index b3c57426..92fa221c 100644 --- a/ncl/ncl_entries/great_circle.ipynb +++ b/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)" ] }, { @@ -43,9 +44,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "1940855.984630871" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "from pyproj import Geod\n", "\n", @@ -76,9 +88,19 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "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", @@ -93,6 +115,59 @@ "print(f\"Z = {cart_coords.z.value}\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## csc2s\n", + "NCL's `csc2s` converts Cartesian coordaintes to spherical (latitude/longitude) coordinates on a unit sphere" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Grab and Go" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Latitude = 40.015\n", + "Longitude = -105.27049999999997\n" + ] + } + ], + "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": {}, diff --git a/ncl/ncl_raw/great_circle.ncl b/ncl/ncl_raw/great_circle.ncl index b3093f69..1fdd9c05 100644 --- a/ncl/ncl_raw/great_circle.ncl +++ b/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("Cartesian X, Cartesian Y, Cartesian Z, Latitude (Degree), 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(x + "," + y + "," + z + "," + sph(0,0) + "," + sph(1,0)) + end + end do +end do diff --git a/ncl/receipts/great_circle.ipynb b/ncl/receipts/great_circle.ipynb index 68300d1c..6a4ee398 100644 --- a/ncl/receipts/great_circle.ipynb +++ b/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)" ] }, { @@ -64,7 +65,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "2ba99b37-a743-4776-bc7b-d5a08b977642", "metadata": {}, "outputs": [], @@ -133,7 +134,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "35abde81-5843-4504-8e32-a137ee1aa094", "metadata": {}, "outputs": [], @@ -147,13 +148,13 @@ "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))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "ee104ea1-e287-4635-b404-5b06ccfb6949", "metadata": {}, "outputs": [], @@ -175,6 +176,62 @@ " astropy_css2c[pair] = cart_coords.xyz.value" ] }, + { + "cell_type": "markdown", + "id": "399d047d-f22c-41cf-996d-d84e1f5b096c", + "metadata": {}, + "source": [ + "### csc2s" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "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(\"../ncl_raw/csc2s_output.txt\", delimiter=',', skiprows=6)\n", + "\n", + "cart_values = csc2s_data[::, 0:3]\n", + "lat_lon = tuple(map(tuple, (csc2s_data[::, 3:])))\n", + "ncl_csc2s = dict(zip(lat_lon, cart_values))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "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[tuple(xyz)] = lat_lon" + ] + }, { "cell_type": "markdown", "id": "3237a0bffc6827fc", @@ -204,10 +261,57 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "74362fd9-0e9f-4cf9-91da-08cd81be625c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "VALID:\n", + "\tBoulder, Boston, Houston: \n", + "\tncl:\t\t1940856\n", + "\tpython:\t\t1940855.984630871\n", + "\n", + "VALID:\n", + "\tFour Corners of Colorado: \n", + "\tncl:\t\t250007.6\n", + "\tpython:\t\t250007.81588628562\n", + "\n", + "VALID:\n", + "\tCaltech, Alberta, Greenwich, Paris, Madrid: \n", + "\tncl:\t\t11634800\n", + "\tpython:\t\t11634798.025760762\n", + "\n", + "VALID:\n", + "\tCrossing the Equator: \n", + "\tncl:\t\t114894.8\n", + "\tpython:\t\t114894.73868013734\n", + "\n", + "VALID:\n", + "\tCrossing the Prime Meridian: \n", + "\tncl:\t\t54450.39\n", + "\tpython:\t\t54450.41832867822\n", + "\n", + "VALID:\n", + "\tHalf of the World: \n", + "\tncl:\t\t255032000\n", + "\tpython:\t\t255031995.77390912\n", + "\n", + "INVALID:\n", + "\tSingle Point -> Invalid NCL: \n", + "\t\tncl:\t\t-127516000\n", + "\t\tpython:\t\t0.0\n", + "\n", + "VALID:\n", + "\tSingle Degree: \n", + "\tncl:\t\t9401.705\n", + "\tpython:\t\t9401.704877506347\n", + "\n" + ] + } + ], "source": [ "for key in ncl_results.keys():\n", " if key in python_results.keys() and key in python_results.keys():\n", @@ -236,7 +340,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "d5738d2d-5abe-4ed4-8e90-c42881c2bfb0", "metadata": {}, "outputs": [], @@ -247,6 +351,26 @@ " 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": "code", + "execution_count": 8, + "id": "53360a09-f1f3-4b0f-b7b4-9501aa5c92e1", + "metadata": {}, + "outputs": [], + "source": [ + "for i, key in enumerate(ncl_csc2s.keys()):\n", + " assert abs(key[0] - astropy_csc2s[list(astropy_csc2s.keys())[i]][0]) < 0.0005\n", + " assert abs(key[1] - astropy_csc2s[list(astropy_csc2s.keys())[i]][1]) < 0.0005" + ] } ], "metadata": {