diff --git a/environment.yml b/environment.yml index d757600f..cd286364 100644 --- 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 index 42d4cbf1..df8a6f9c 100644 --- a/ncl/ncl_entries/great_circle.ipynb +++ b/ncl/ncl_entries/great_circle.ipynb @@ -15,7 +15,8 @@ "\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)" ] }, { @@ -23,7 +24,7 @@ "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", "
Important Note
\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": {}, diff --git a/ncl/receipts/great_circle.ipynb b/ncl/receipts/great_circle.ipynb index a89522a8..cc7392c8 100644 --- 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)" ] }, { @@ -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": {