HAFS#

This tutorial shows you how to download HAFS data from NOMADS with the Herbie Python package.

The HAFS model comes in two flavors

HAFS-A; In Herbie, set model="hafsa"
HAFS-B; In Herbie, set model="hafsb"

What’s the difference? According to AndyHazelton

They use different physics (HAFS-A GFDL microphysics, HAFS-B Thompson, and HAFS-B uses the tc-pbl option @XiaominChen7 and I worked on).

You must supply the storm=<stormName_or_ID> to Herbie.

[1]:

from herbie import Herbie
from toolbox import EasyMap, pc
from paint.standard2 import cm_tmp, cm_wind

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

You can get a Herbie for a storm by its ID or Name

[21]:

H = Herbie("2023-9-9", model="hafsa", fxx=12, storm="lee")
H.storm, H.storm_name

✅ Found ┊ model=hafsa ┊ product=storm.atm ┊ 2023-Sep-09 00:00 UTC F12 ┊ GRIB2 @ nomads ┊ IDX @ nomads

[21]:

('13l', 'lee')

[22]:

H = Herbie("2023-9-9", model="hafsa", fxx=12, storm="11e")
H.storm, H.storm_name

✅ Found ┊ model=hafsa ┊ product=storm.atm ┊ 2023-Sep-09 00:00 UTC F12 ┊ GRIB2 @ nomads ┊ IDX @ nomads

[22]:

('11e', 'jova')

Ok, let’s get and plot some HAFS data!

[2]:

H = Herbie("2023-9-9", model="hafsa", fxx=12, storm="lee")

✅ Found ┊ model=hafsa ┊ product=storm.atm ┊ 2023-Sep-09 00:00 UTC F12 ┊ GRIB2 @ nomads ┊ IDX @ nomads

[3]:

H.inventory("WIND")

[3]:

	grib_message	start_byte	end_byte	range	reference_time	valid_time	variable	level	forecast_time	search_this
698	699	220591614	220881145	220591614-220881145	2023-09-09	2023-09-09 12:00:00	WIND	10 m above ground	11-12 hour max fcst	:WIND:10 m above ground:11-12 hour max fcst

[4]:

H.PRODUCTS

[4]:

{'storm.atm': '13L-Lee',
 'storm.sat': '13L-Lee',
 'parent.atm': '13L-Lee',
 'parent.sat': '13L-Lee',
 'parent.swath': '13L-Lee',
 'ww3': '13L-Lee'}

[5]:

ds = H.xarray("WIND")
ds

curl -s --range 220591614-220881145 "https://nomads.ncep.noaa.gov/pub/data/nccf/com/hafs/prod/hfsa.20230909/00/13l.2023090900.hfsa.storm.atm.f012.grb2" > "/home/blaylock/data/hafsa/20230909/subset_71128152__13l.2023090900.hfsa.storm.atm.f012.grb2"

[6]:

ax = (
    EasyMap("50m", crs=ds.herbie.crs, figsize=[8, 8], dark=True)
    .OCEAN()
    .LAND()
    .STATES()
    .ax
)
p = ax.pcolormesh(
    ds.longitude, ds.latitude, ds.si10, transform=pc, **cm_wind().cmap_kwargs
)
plt.colorbar(p, ax=ax, orientation="horizontal", pad=0.05, **cm_wind().cbar_kwargs)

ax.set_title(
    f"{ds.model.upper()}: {H.product_description} ({H.product})\nValid: {ds.valid_time.dt.strftime('%H:%M UTC %d %b %Y').item()}",
    loc="left",
)
ax.set_title(ds.si10.GRIB_name, loc="right")

[6]:

Text(1.0, 1.0, '10 metre wind speed')

../../../_images/user_guide_tutorial_model_notebooks_hafs_10_1.png

You can also get the HAFS-B data

[7]:

H = Herbie("2023-9-9", model="hafsb", fxx=12, storm="lee")

✅ Found ┊ model=hafsb ┊ product=storm.atm ┊ 2023-Sep-09 00:00 UTC F12 ┊ GRIB2 @ nomads ┊ IDX @ nomads

[8]:

dsb = H.xarray("WIND")
ax = (
    EasyMap("50m", crs=dsb.herbie.crs, figsize=[8, 8], dark=True)
    .OCEAN()
    .LAND()
    .STATES()
    .ax
)
p = ax.pcolormesh(
    dsb.longitude, dsb.latitude, dsb.si10, transform=pc, **cm_wind().cmap_kwargs
)
plt.colorbar(p, ax=ax, orientation="horizontal", pad=0.05, **cm_wind().cbar_kwargs)

ax.set_title(
    f"{dsb.model.upper()}: {H.product_description} ({H.product})\nValid: {dsb.valid_time.dt.strftime('%H:%M UTC %d %b %Y').item()}",
    loc="left",
)
ax.set_title(dsb.si10.GRIB_name, loc="right")

curl -s --range 217405795-217676933 "https://nomads.ncep.noaa.gov/pub/data/nccf/com/hafs/prod/hfsb.20230909/00/13l.2023090900.hfsb.storm.atm.f012.grb2" > "/home/blaylock/data/hafsb/20230909/subset_71128152__13l.2023090900.hfsb.storm.atm.f012.grb2"

[8]:

Text(1.0, 1.0, '10 metre wind speed')

../../../_images/user_guide_tutorial_model_notebooks_hafs_13_2.png

And we can show the difference between the two models

[9]:

ax = (
    EasyMap("50m", crs=dsb.herbie.crs, figsize=[8, 8], dark=True)
    .OCEAN()
    .LAND()
    .STATES()
    .ax
)
p = ax.pcolormesh(
    ds.longitude,
    ds.latitude,
    (ds.si10.values - dsb.si10.values),
    transform=pc,
    cmap="bwr",
    vmax=8,
    vmin=-8,
)
plt.colorbar(p, ax=ax, orientation="horizontal", pad=0.05, label="Wind Speed")

ax.set_title(
    f"HAVS-A minus HAVS-B\nValid: {dsb.valid_time.dt.strftime('%H:%M UTC %d %b %Y').item()}",
    loc="left",
)
ax.set_title(dsb.si10.GRIB_name, loc="right")

[9]:

Text(1.0, 1.0, '10 metre wind speed')

../../../_images/user_guide_tutorial_model_notebooks_hafs_15_1.png

[10]:

ds

There might be a better way to do this, but you can drop the NAN values around the perimeter of the domain like this…

[11]:

ds.dropna(dim="latitude", how="all").dropna(dim="longitude", how="all").si10.plot()

[11]:

<matplotlib.collections.QuadMesh at 0x7efe5d59d7d0>

../../../_images/user_guide_tutorial_model_notebooks_hafs_18_1.png

What if I don’t know what “storms” are available?#

[15]:

from herbie.models.hafs import Storms

[18]:

# Get a dict of storm ids and storm names.
Storms().id_to_name

[18]:

{'99w': 'invest',
 '14l': 'margot',
 '13l': 'lee',
 '10w': 'haikui',
 '12w': 'yun-yeung',
 '09w': 'saola',
 '11e': 'jova'}

[14]:

# herbie will also tell you if you ask for a storm that is unknown
H = Herbie("2023-9-9", model="hafsa", fxx=12, storm="test")

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[14], line 2
      1 # herbie will tell you if you put a storm that is unknown
----> 2 H = Herbie("2023-9-9", model="hafsa", fxx=12, storm="test")

File ~/GITHUB/Herbie/herbie/core.py:256, in Herbie.__init__(self, date, valid_date, model, fxx, product, priority, save_dir, overwrite, verbose, **kwargs)
    248     setattr(self, key, value)
    250 # Get details from the template of the specified model.
    251 # This attaches the details from the `models.<model>.template`
    252 # class to this Herbie object.
    253 # This line is equivalent to `model_templates.gfs.template(self)`.
    254 # I do it this way because the model name is a variable.
    255 # (see https://stackoverflow.com/a/7936588/2383070 for what I'm doing here)
--> 256 getattr(model_templates, self.model).template(self)
    258 if product is None:
    259     # The user didn't specify a product, so let's use the first
    260     # product in the model template.
    261     self.product = list(self.PRODUCTS)[0]

File ~/GITHUB/Herbie/herbie/models/hafs.py:55, in hafsa.template(self)
     53     self.storm = S.name_to_id.get(self.storm.lower())
     54     if self.storm is None:
---> 55         raise ValueError(f"`storm` should be one of {S.name_to_id.keys()}")
     57 self.storm_name = S.id_to_name.get(self.storm)
     59 self.PRODUCTS = {
     60     "storm.atm": f"{self.storm.upper()}-{self.storm_name.title()}",
     61     "storm.sat": f"{self.storm.upper()}-{self.storm_name.title()}",
   (...)
     65     "ww3": f"{self.storm.upper()}-{self.storm_name.title()}",
     66 }

ValueError: `storm` should be one of dict_keys(['invest', 'margot', 'lee', 'haikui', 'yun-yeung', 'saola', 'jova'])