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 herbie.toolbox import EasyMap, pc
from herbie import paint
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"
[5]:
<xarray.Dataset> Dimensions: (latitude: 801, longitude: 1001) Coordinates: time datetime64[ns] 2023-09-09 step timedelta64[ns] 12:00:00 heightAboveGround float64 10.0 * latitude (latitude) float64 12.12 12.13 12.15 ... 28.09 28.11 * longitude (longitude) float64 292.3 292.4 292.4 ... 312.3 312.3 valid_time datetime64[ns] 2023-09-09T12:00:00 Data variables: si10 (latitude, longitude) float32 nan nan nan ... nan nan gribfile_projection object None Attributes: GRIB_edition: 2 GRIB_centre: kwbc GRIB_centreDescription: US National Weather Service - NCEP GRIB_subCentre: 0 Conventions: CF-1.7 institution: US National Weather Service - NCEP model: hafsa product: storm.atm description: Hurricane Analysis and Forecast System (HAFS-A) ... remote_grib: https://nomads.ncep.noaa.gov/pub/data/nccf/com/h... local_grib: /home/blaylock/data/hafsa/20230909/subset_711281... search: WIND
[6]:
ax = (
EasyMap("50m", crs=ds.herbie.crs, figsize=[8, 8], theme="dark")
.OCEAN()
.LAND()
.STATES()
.ax
)
p = ax.pcolormesh(
ds.longitude, ds.latitude, ds.si10, transform=pc, **paint.NWSWindSpeed.kwargs2
)
plt.colorbar(p, ax=ax, orientation="horizontal", pad=0.05, **paint.NWSWindSpeed.cbar_kwargs2)
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')

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], theme="dark")
.OCEAN()
.LAND()
.STATES()
.ax
)
p = ax.pcolormesh(
dsb.longitude, dsb.latitude, dsb.si10, transform=pc, **paint.NWSWindSpeed.kwargs2
)
plt.colorbar(p, ax=ax, orientation="horizontal", pad=0.05, **paint.NWSWindSpeed.cbar_kwargs2)
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')

And we can show the difference between the two models
[9]:
ax = (
EasyMap("50m", crs=dsb.herbie.crs, figsize=[8, 8], theme="dark")
.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')

[10]:
ds
[10]:
<xarray.Dataset> Dimensions: (latitude: 801, longitude: 1001) Coordinates: time datetime64[ns] 2023-09-09 step timedelta64[ns] 12:00:00 heightAboveGround float64 10.0 * latitude (latitude) float64 12.12 12.13 12.15 ... 28.09 28.11 * longitude (longitude) float64 292.3 292.4 292.4 ... 312.3 312.3 valid_time datetime64[ns] 2023-09-09T12:00:00 Data variables: si10 (latitude, longitude) float32 nan nan nan ... nan nan gribfile_projection object None Attributes: GRIB_edition: 2 GRIB_centre: kwbc GRIB_centreDescription: US National Weather Service - NCEP GRIB_subCentre: 0 Conventions: CF-1.7 institution: US National Weather Service - NCEP model: hafsa product: storm.atm description: Hurricane Analysis and Forecast System (HAFS-A) ... remote_grib: https://nomads.ncep.noaa.gov/pub/data/nccf/com/h... local_grib: /home/blaylock/data/hafsa/20230909/subset_711281... search: WIND
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>

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'])