How can I get geopotential height and wind data at all pressure levels in HRRR for a specific location?

[blaylockbk]

Was asked this question in an email.

[blaylockbk]

Getting data from a specific location isn't super efficient with GRIB2 data and Heribe, but it can be done. It requires downloading the full HRRR grid for each variable you want.

Download subsets of GRIB2 files.

from herbie import Herbie

H = Herbie('2023-01-01', product="prs", fxx=1)

Note: Hopefully you have a good internet connection to help the download speed.

# Download subset of file with all geopotential heights
ds_hgt = H.xarray(
    searchString="HGT:.* mb:",
    verbose=True,
    remove_grib=False,  # Set to false to retain the file locally. (We don't want to redownload it next time)
)
ds_hgt

# Download subset of file with all U and V winds
ds_wind = H.xarray(
    searchString="[U|V]GRD:.* mb:",
    verbose=True,
    remove_grib=False,  # Set to false to retain the file locally.
)
)
ds_wnd

Note that we set remove_grib=False to retain the file locally; we don't want to redownload it next time we want to load the data.

Get data at a specific point.

Herbie has an accessor to get data at a point nearest a lat/lon location.

I uncovered a small bug while working through this myself I'll fix, but the below does work correctly.

point_wnd = ds_wnd.herbie.nearest_points(
    points=(-79.6, 31.5)  # This is a tuple (longitude, latitude)
)

This returns an xarray dataset for the winds nearest the requested point for all the vertical pressure levels.

If you look at the longitude and latitude data from the grid retrieved, you can see that the point is in fact very close to the point we requested.

# Note: Subtract 360 degrees from longitude to put it in -180 - 180 degrees units
point_wnd.longitude.item()-360, point_wnd.latitude.item()

OUT: (-79.61449010405101, 31.49582934470579)

Similarly, you can query the data at multiple points.

point_wnd = ds_wnd.herbie.nearest_points(
    points=[(-79.6, 31.5), (-100, 45), (-95, 35.3)],
)

Now see that the "point" dimension is "3"

If you want to plot these values time, you should generate a list of these "point_wnd" datasets at different times or forecast lead times and then concatenate the datasets together. Then you can do some neat statistics and plotting. I'll let you work through that yourself if you want.

[blaylockbk]

UPDATE:
Since in this case the variables are all on pressure levels, you can actually do the download in one go and get all three variables (HGT, UGRD, and VGRD) in one Dataset.

H = Herbie('2023-01-01', product="prs", fxx=1)
ds = H.xarray(":(?:HGT|UGRD|VGRD):.* mb:", verbose=True, remove_grib=False)
ds

On my windows laptop with slowish WiFi connection, this took ~3 minutes. There is overwhead with each partial range get (how Herbie does the subsetting) so you might experiment if downloading the full GRIB2 file before subsetting gives you better performance. That is, do H.download() before you do H.xarray(...).