GFS GRIB filter example

.gitignore

@@ -1,4 +1,6 @@
 *.pyc
 .vscode/
 *.nc
-*.png
+*.png
+*.grb2
+*.idx

environment.yml

@@ -10,6 +10,7 @@ dependencies:
   - cartopy
   - typer
   - netCDF4
+  - cfgrib
   - pip:
     - pydap==3.2.2
     - xarray==2022.11.0

notebook/get_gfs_filter.ipynb

@@ -0,0 +1,181 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# GFS filter service\n",
+    "\n",
+    "This notebooks contains an example of downloading GFS data using the filter service, that works with the original GRIB files."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import datetime as dt\n",
+    "from pathlib import Path\n",
+    "import urllib.request as request\n",
+    "from contextlib import closing\n",
+    "from pathlib import Path\n",
+    "import shutil\n",
+    "import numpy as np\n",
+    "\n",
+    "range1 = lambda start, end, step=1: np.arange(start, end + step, step)\n",
+    "\n",
+    "\n",
+    "def get_gfs_grib_filter(\n",
+    "    date: dt.datetime,\n",
+    "    run: int,\n",
+    "    vars: set,\n",
+    "    levels: set,\n",
+    "    lat: tuple[float, float] = (-90, 90),\n",
+    "    lon: tuple[float, float] = (-180, 180),\n",
+    "    future_window: int = 8,\n",
+    "    outdir: str = \".\",\n",
+    "):\n",
+    "    \"\"\"Download meteorological data from GFS HTTP server using the grib filter\n",
+    "    service. This allows subsetting coordinates, time steps, variables and levels\n",
+    "\n",
+    "    References\n",
+    "    ----------\n",
+    "    [1] https://nomads.ncep.noaa.gov/\n",
+    "    [2] https://www.nco.ncep.noaa.gov/pmb/products/gfs/\n",
+    "    \"\"\"\n",
+    "    base_url = \"https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl\"\n",
+    "    fout = Path(outdir) / f\"gfs_{date:%Y%m%d}_{run:02d}.grb2\"\n",
+    "\n",
+    "    # Transform lon from -180, 180 to 0, 360. The URL of the requests\n",
+    "    # changes whether it is the full region or a sub-region\n",
+    "    lon0, lon1 = lon[0] % 360, lon[1] % 360\n",
+    "\n",
+    "    if lon0 == lon1 == 180 or lon0 > lon1:\n",
+    "        lon0 = 0\n",
+    "        lon1 = 360\n",
+    "        prefix = \"\"\n",
+    "    else:\n",
+    "        prefix = \"subregion=&\"\n",
+    "\n",
+    "    dir = f\"%2Fgfs.{date:%Y%m%d}%2F{run:02d}%2Fatmos\"\n",
+    "\n",
+    "    coords_ = (\n",
+    "        f\"{prefix}leftlon={lon0}&rightlon={lon1}&toplat={lat[1]}&bottomlat={lat[0]}\"\n",
+    "    )\n",
+    "    levels_ = \"&\".join([f\"lev_{level}=on\" for level in levels])\n",
+    "    vars_ = \"&\".join([f\"var_{var}=on\" for var in vars])\n",
+    "\n",
+    "    # Until hour 120 (5 days), time resolution is 1 hour. From hour 120 to 384 (16 days),\n",
+    "    # time resolution is 3 hours, so we need to combine two ranges in the last case\n",
+    "    last_hour = future_window * 24\n",
+    "    if last_hour <= 120:\n",
+    "        time_range = range1(0, last_hour)\n",
+    "    else:\n",
+    "        time_range = np.concatenate((range1(0, 120), range1(123, last_hour, 3)))\n",
+    "\n",
+    "    try:\n",
+    "        with open(fout, \"wb\") as f:\n",
+    "            for time in time_range:\n",
+    "                file = f\"gfs.t{run:02d}z.pgrb2.0p25.f{time:03d}\"\n",
+    "                url = f\"{base_url}?file={file}&{levels_}&{vars_}&{coords_}&dir={dir}\"\n",
+    "\n",
+    "                # Print URL for debugging\n",
+    "                if time == 0:\n",
+    "                    print(url)\n",
+    "\n",
+    "                with closing(request.urlopen(url)) as r:\n",
+    "                    shutil.copyfileobj(r, f)\n",
+    "    except Exception as err:\n",
+    "        fout.unlink(missing_ok=True)\n",
+    "        raise IOError(f\"{err}\")\n",
+    "\n",
+    "    print(f'File \"{fout.name}\" downloaded')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are 4 runs per day (00, 06, 12, 18). Information about available variables can be found [here](https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.f000.shtml). Example: download wind speed forecasts (u and v components) 10 meters above ground for today (24 values, one per hour, starting at 00)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f000&lev_10_m_above_ground=on&var_VGRD=on&var_UGRD=on&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%2Fgfs.20221109%2F00%2Fatmos\n",
+      "File \"gfs_20221109_00.grb2\" downloaded\n"
+     ]
+    }
+   ],
+   "source": [
+    "date = dt.date.today()\n",
+    "run = 0\n",
+    "vars = {\"UGRD\", \"VGRD\"}\n",
+    "levels = {\"10_m_above_ground\"}\n",
+    "\n",
+    "get_gfs_grib_filter(date, run, vars, levels, future_window=1)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Depending on the variables and levels that you download, you can open the GRIB file with xarray and export it to NetCDF. See [this](https://github.com/albertotb/get-gfs/issues/9#issuecomment-886061535) comment for more information."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Ignoring index file 'gfs_20221109_00.grb2.923a8.idx' older than GRIB file\n"
+     ]
+    }
+   ],
+   "source": [
+    "import xarray as xr\n",
+    "\n",
+    "with xr.open_dataset(f\"gfs_{date:%Y%m%d}_{run:02d}.grb2\", engine=\"cfgrib\") as ds:\n",
+    "    ds.to_netcdf(f\"gfs_{date:%Y%m%d}_{run:02d}.nc\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.10.6 ('get-gfs')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.0"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "e6d217ddd860d14589860f0e06cb2192a9028cf6c110b539310ed5bf6d43e16d"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}