Add rotated_ll to supported CF CRS grid types

check_pygrib_vs_herbie_crs_extraction.ipynb

@@ -15,19 +15,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
     "from herbie import Herbie\n",
     "import pygrib\n",
+    "import pandas as pd\n",
     "from pyproj import CRS\n",
     "from herbie.crs import get_cf_crs"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
@@ -49,7 +50,6 @@
       "  pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}\n",
       "  Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}\n",
       "  equal= True\n",
-      "👨🏻‍🏭 Created directory: [/home/blaylock/data/graphcast/20250101]\n",
       "\n",
       "Model: GRAPHCAST\n",
       "  pygrib {'a': 4326.0, 'b': 4326.0, 'proj': 'longlat'}\n",
@@ -81,16 +81,6 @@
       "  Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}\n",
       "  equal= True\n",
       "\n",
-      "Model: HAFSA\n",
-      "  pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}\n",
-      "  Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}\n",
-      "  equal= True\n",
-      "\n",
-      "Model: HREF\n",
-      "  pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}\n",
-      "  Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}\n",
-      "  equal= True\n",
-      "\n",
       "Model: NAM\n",
       "  pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}\n",
       "  Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}\n",
@@ -104,6 +94,11 @@
       "Model: RTMA\n",
       "  pygrib {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}\n",
       "  Herbie {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}\n",
+      "  equal= True\n",
+      "\n",
+      "Model: HRDPS\n",
+      "  pygrib {'a': 6371229, 'b': 6371229, 'lon_0': 245.305142, 'o_lat_p': 36.08852, 'o_lon_p': 0.0, 'o_proj': 'longlat', 'proj': 'ob_tran'}\n",
+      "  Herbie {'a': 6371229, 'b': 6371229, 'lon_0': 245.305142, 'o_lat_p': 36.08852, 'o_lon_p': 0.0, 'o_proj': 'longlat', 'proj': 'ob_tran'}\n",
       "  equal= True\n"
      ]
     }
@@ -119,24 +114,31 @@
     "    (dict(model=\"aifs\"), \":2t:\"),\n",
     "    (dict(model=\"gefs\", member=1), \":TMP:2 m above\"),\n",
     "    (dict(model=\"gefs\", product=\"wave\", member=1), \":WIND:surf\"),\n",
-    "    (dict(model=\"hafsa\", product=\"storm.atm\", storm=\"07s\"), \":TMP:2 m above\"),\n",
-    "    (dict(model=\"href\", products=\"mean\", domain=\"conus\", fxx=1), \":TMP:2 m above\"),\n",
+    "    # FAILING\n",
+    "    # (dict(model=\"hafsa\", product=\"storm.atm\", storm=\"07s\"), \":TMP:2 m above\"),\n",
+    "    # FAILING\n",
+    "    # (dict(model=\"href\", products=\"mean\", domain=\"conus\", fxx=1), \":TMP:2 m above\"),\n",
     "    (dict(model=\"nam\"), \":TMP:2 m above\"),\n",
     "    (dict(model=\"urma\"), \":TMP:2 m above\"),\n",
     "    (dict(model=\"rtma\"), \":TMP:2 m above\"),\n",
+    "    (dict(model=\"hrdps\", fxx=0, product=\"continental/2.5km\", variable=\"TMP\", level=\"AGL-2m\"), None),\n",
+    "\n",
     "]:\n",
     "    if model[\"model\"] in (\"hafsa\", \"href\"):  # only on nomads\n",
     "        date = \"2025-01-16 06:00\"\n",
+    "    elif model[\"model\"] in (\"hrdps\",): # no historical data\n",
+    "        date = pd.Timestamp(\"now\").floor(\"6h\") - pd.Timedelta(\"6h\")\n",
     "    else:\n",
-    "        date = \"2025-01-01\"\n",
+    "        date = \"2025-07-05T06:00:00\"\n",
     "\n",
-    "    ds = Herbie(date, verbose=False, **model).xarray(\n",
-    "        search, remove_grib=False, _use_pygrib_for_crs=False\n",
-    "    )\n",
+    "    # ds = Herbie(date, verbose=False, **model).xarray(\n",
+    "    #    search, remove_grib=False, _use_pygrib_for_crs=False\n",
+    "    # )\n",
+    "    ds = Herbie(date, verbose=False, **model).download()\n",
     "\n",
     "    print()\n",
     "    print(f\"Model: {model['model'].upper()}\")\n",
-    "    with pygrib.open(str(ds.local_grib)) as grb:\n",
+    "    with pygrib.open(str(ds)) as grb:\n",
     "        msg = grb.message(1)\n",
     "        projparams_dict_pygrib = dict(sorted(msg.projparams.items()))\n",
     "        print(\"  pygrib\", projparams_dict_pygrib)\n",
@@ -145,7 +147,7 @@
     "        sorted(get_cf_crs(ds, _return_projparams=True).items())\n",
     "    )\n",
     "    print(\"  Herbie\", projparams_dict_herbie)\n",
-    "    print(\"  equal=\", projparams_dict_herbie == projparams_dict_pygrib)\n"
+    "    print(\"  equal=\", projparams_dict_herbie == projparams_dict_pygrib)"
    ]
   },
   {
@@ -5787,7 +5789,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "herbie-dev",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -5801,9 +5803,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.1"
+   "version": "3.12.9"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

pyproject.toml

@@ -36,7 +36,7 @@ dependencies = [
     "pandas>=2.1",
     "pyproj>=3.7.0",
     "requests>=2.23.3",
-    "toml>=0.10.2", # TODO: Drop in favor of tomllib when Python >=3.11 is required.
+    "toml>=0.10.2",     # TODO: Drop in favor of tomllib when Python >=3.11 is required.
     "xarray>=2025.1.1",
 ]
 dynamic = ["version"]
@@ -55,6 +55,7 @@ build-backend = "hatchling.build"
 [dependency-groups]
 dev = [
     "ipykernel>=6.29.5",
+    "jupyter>=1.1.1",
     "pytest>=8.3.5",
     "pytest-cov>=6.1.1",
     "ruff>=0.11.6",

src/herbie/crs.py

@@ -6,9 +6,8 @@
 
 from typing import Any
 
-from pyproj import CRS
-
 import xarray as xr
+from pyproj import CRS
 
 
 def get_cf_crs(
@@ -82,6 +81,15 @@ def get_cf_crs(
         projparams["a"] = a
         projparams["b"] = b
 
+    elif da.GRIB_gridType == "rotated_ll":
+        projparams = {"proj": "ob_tran"}
+        projparams["o_proj"] = "longlat"
+        projparams["a"] = a
+        projparams["b"] = b
+        projparams["lon_0"] = da.GRIB_longitudeOfSouthernPoleInDegrees
+        projparams["o_lon_p"] = 0.0
+        projparams["o_lat_p"] = -da.GRIB_latitudeOfSouthernPoleInDegrees
+
     elif da.GRIB_gridType == "regular_gg":
         projparams = {"proj": "longlat"}
         projparams["a"] = a

tests/test_hrdps.py

@@ -0,0 +1,76 @@
+## Nikhil Shankar
+## August 10, 2025
+
+"""
+Tests for downloading HRDPS model.
+
+DRAGONS:
+
+At ~00Z the HRDPS servers are cleared of the previous day's forecasts.
+The next forecast arrives at ~06Z. During the interim period these tests will fail.
+"""
+
+from datetime import datetime
+
+import pandas as pd
+
+from herbie import Herbie, config
+
+now = datetime.now()
+latest = pd.Timestamp("now").floor("6h") - pd.Timedelta("6h")
+
+save_dir = config["default"]["save_dir"] / "Herbie-Tests-Data/"
+
+# Remove all previous test data
+for i in (save_dir / "hrdps").rglob("*"):
+    if i.is_file():
+        i.unlink()
+
+
+def test_hrdps_download():
+    H = Herbie(
+        latest,
+        model="hrdps",
+        product="continental/2.5km",
+        variable="TMP",
+        level="AGL-2m",
+        overwrite=True,
+        save_dir=save_dir,
+    )
+    f = H.download()
+    assert H.get_localFilePath().exists()
+    f.unlink()
+
+
+def test_hrdps_xarray():
+    H = Herbie(
+        latest,
+        model="hrdps",
+        product="continental/2.5km",
+        variable="TMP",
+        level="AGL-2m",
+        overwrite=True,
+        save_dir=save_dir,
+    )
+    H.xarray(remove_grib=False)
+    assert H.get_localFilePath().exists()
+    H.get_localFilePath().unlink()
+
+
+def test_hrdps_to_netcdf():
+    """Check that a xarray Dataset can be written to a NetCDF file.
+
+    It is important that I have haven't put any python objects in the
+    xarray Dataset attributes.
+    """
+    H = Herbie(
+        latest,
+        model="hrdps",
+        product="continental/2.5km",
+        variable="TMP",
+        level="AGL-2m",
+        overwrite=True,
+        save_dir=save_dir,
+    )
+    ds = H.xarray(remove_grib=False)
+    ds.to_netcdf(save_dir / "test_hrdps_to_netcdf.nc")