THREDDS-served GRIB Data: HRRR (Regional Grid) Model¶

This notebook reads in realtime HRRR from a THREDDS server.¶

Overview:¶

Determine latest available model run
Use dictionaries to deal with coordinate dimension names that change from run-to-run
Subset data and map regions
Create a time series of a variable at a point

Imports¶

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import pyplot as plt
from matplotlib.dates import DateFormatter, AutoDateLocator,HourLocator,DayLocator
import numpy as np
import xarray as xr
import pandas as pd
from datetime import datetime, timedelta
from pyproj import Proj
import seaborn as sns
sns.set()

Access real-time NCEP HRRR model output from Unidata’s THREDDS server¶

Parse the current time and choose a recent hour when we might expect the model run to be complete. There is normally a ~90 minute lag, but we’ll give it a bit more than that.

now = datetime.now()
year = now.year
month = now.month
day = now.day
hour = now.hour
minute = now.minute
print (year, month, day, hour,minute)

2023 10 24 18 20

# The previous hour's data is typically available shortly after the half-hour
if (minute < 45): 
    hrDelta = 2
else:
    hrDelta = 1
runTime = now - timedelta(hours=hrDelta)
runTimeStr = runTime.strftime('%Y%m%d %H00 UTC')
modelDate = runTime.strftime('%Y%m%d')
modelHour = runTime.strftime('%H')
modelMin = runTime.strftime('%M')
modelDay = runTime.strftime('%D')
titleTime = modelDay + ' ' +  modelHour + '00 UTC'
print (modelDay)
print (modelHour)
print (runTimeStr)
print(titleTime)

10/24/23
16
20231024 1600 UTC
10/24/23 1600 UTC

Open the Dataset

ds = xr.open_dataset(f'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_{modelDate}_{modelHour}00.grib2')

Leverage MetPy’s parse_cf function, which among other things will make it easy to extract the dataset’s map projection.

ds = ds.metpy.parse_cf()

Note that there now exists a metpy_crs coordinate variable, which we’ll soon use to determine and use the Dataset’s map projection.

Select one of the variables from the Dataset; in this case, we’ll choose 2-meter temperature.

var = ds.Temperature_height_above_ground

Examine the grid

GRIB files that are served by THREDDS are automatically converted into NetCDF. A quirk in this translation can result in certain dimensions (typically those that are time or vertically-based) varying from run-to-run. For example, the time coordinate dimension for a given variable might be time in one run, but time1 (or time2, time3, etc.) in another run. This poses a challenge any time we want to subset a DataArray. For example, if a particular dataset had vertical levels in units of hPa, the dimension name for the geopotential height grid might be isobaric. If we wanted to select only the 500 and 700 hPa levels, we might write the following line of code:

Z = ds.sel(isobaric=[500, 700])

This would fail if we pointed to model output from a different time, and the vertical coordinate name for geopotential height was instead isobaric1.

Via strategic use of Python dictionaries, we can programmatically deal with this!

First, let’s use some MetPy accessors that determine the dimension names for times and vertical levels:

timeDim, vertDim = var.metpy.time.name, var.metpy.vertical.name
timeDim, vertDim

('time2', 'height_above_ground3')

We next construct dictionaries corresponding to each desired subset of the relevant dimensions. Here, we specify the first time, and the one and only vertical level:

idxTime = 0 # First time
idxVert = 0 # First (and in this case, only) vertical level

timeDict = {timeDim: idxTime}
vertDict = {vertDim: idxVert}

timeDict, vertDict

({'time2': 0}, {'height_above_ground3': 0})

Produce a quick visualization at the given time and level. Pass the two dictionaries into Xarray’s isel function.

grid = var.isel(timeDict).isel(vertDict)
grid.plot()

<matplotlib.collections.QuadMesh at 0x14b79b674250>

../../_images/062eceaca4e88695c069d3c314905b9e7bd3fa4ae10920d2f6e3bf2ed7bfd1be.png

Get the Dataset’s map projection.

proj_data = var.metpy.cartopy_crs
proj_data;

Plot the data over the full domain

Create a title string

validTime = var.isel(timeDict)[timeDim].values # Returns validTime as a NumPy Datetime 64 object
validTime = pd.Timestamp(validTime) # Converts validTime to a Pandas Timestamp object
timeStr = validTime.strftime("%Y-%m-%d %H%M UTC") # We can now use Datetime's strftime method to create a time string.
tl1 = "HRRR 2m temperature ($^\circ$C)"
tl2 = str('Valid at: '+ timeStr)
title_line = (tl1 + '\n' + tl2 + '\n')

Convert grid to desired units

grid = grid.metpy.convert_units('degC')

Create DataArray objects for the horizontal coordinate dimensions. In this case, the model is output on a Lambert Conic Conformal projection, which means the horizontal dimensions are named x and y and represent physical distance in meters from the projection’s central point.

x, y = ds.x, ds.y

res = '50m'
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1,projection=proj_data)

ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature (cfeature.LAND.with_scale(res))
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5, zorder=4)
ax.add_feature (cfeature.STATES.with_scale(res))
# Usually avoid plotting counties once the regional extent is large
#ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );
CF = ax.contourf(x,y,grid,cmap='coolwarm',alpha=0.99,transform=proj_data)
cbar = plt.colorbar(CF,fraction=0.046, pad=0.03,shrink=0.5)
cbar.ax.tick_params(labelsize=10)
cbar.ax.set_ylabel(f'Temperature ($^\circ$C)',fontsize=10);
title = ax.set_title(title_line,fontsize=16)

../../_images/f478325c81fc14777b94f520850498df7d8774a358d43d60dd8951e9d5792b88.png

Subset the dataset to a smaller region¶

lonW, lonE, latS, latN = [-80,-71.3,40.2,45]
mapExtent = lonW, lonE, latS, latN

Use the Pyproj library to perform regional subsets¶

Create a pyproj object that corresponds to the dataset’s native projection

pFull = Proj(proj_data)

Define a function that finds the index of the closest gridpoint in x/y coordinates to a specific location in lon-lat coordinates.

def find_closest(array, value):
    idx = (np.abs(array-value)).argmin()
    return idx

Find the x- and y-coordinate indicies that are closest to the values of the four corners.

# Lats/Lons of corner points
SWLon, SWLat = (mapExtent[0]), (mapExtent[2])
NWLon, NWLat = (mapExtent[0], mapExtent[3])
SELon, SELat = (mapExtent[1], mapExtent[2])
NELon, NELat = (mapExtent[1], mapExtent[3])

# Corresponding corner points' coordinates in the dataset's native CRS
xSW, ySW = pFull(SWLon, SWLat)
xNW, yNW = pFull(NWLon, NWLat)
xSE, ySE = pFull(SELon, SELat)
xNE, yNE = pFull(NELon, NELat)

# Find the indicies that correspond to these corner points coordinates
xSWds, ySWds = find_closest(x, xSW), find_closest(y, ySW)
xNWds, yNWds = find_closest(x, xNW), find_closest(y, yNW)
xSEds, ySEds = find_closest(x, xSE), find_closest(y, ySE)
xNEds, yNEds = find_closest(x, xNE), find_closest(y, yNE)

# Get the minimum and maximum values of the x- and y- coordinates
xMin, xMax = np.min([xNWds, xSWds]), np.max([xNEds, xSEds]) + 1
yMin, yMax = np.min([ySEds, ySWds]), np.max([yNEds, yNWds]) + 1

In order for the data region to fill the plotted region, pad the x and y ranges (trial and error, depends on the subregion)

xPad = 90
yPad = 40

xRange = np.arange(xMin - xPad, xMax + xPad)
yRange = np.arange(yMin - yPad, yMax + yPad)

Read in the temperature data over the selected region. Also subset the time and vertical coordinates. Convert the units.

dsSub = ds.isel(x = xRange, y = yRange)
xSub, ySub = dsSub.x, dsSub.y
varSub = dsSub.Temperature_height_above_ground
gridSub = varSub.isel(timeDict).isel(vertDict)
gridSub = gridSub.metpy.convert_units('degC')

Specify the output map projection (Lambert Conformal, centered on New York State)

lon1 = -75
lat1 = 42.5
slat = 42.5

proj_map= ccrs.LambertConformal(central_longitude=lon1,
                             central_latitude=lat1,
                             standard_parallels=[slat,slat])

Plot the regional map. Omit plotting the land/water features since they tend to obscure the nuances of the temperature field, and also slows down the creation of the figure.

res = '10m'
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1,projection=proj_map)
ax.set_extent(mapExtent, crs=ccrs.PlateCarree())
#ax.set_facecolor(cfeature.COLORS['water'])
#ax.add_feature (cfeature.LAND.with_scale(res))
ax.add_feature(cfeature.COASTLINE.with_scale(res))
#ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeature.STATES.with_scale(res))
# Usually avoid plotting counties once the regional extent is large
#ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );
CF = ax.contourf(xSub,ySub,gridSub, cmap='coolwarm',alpha=0.5,transform=proj_data)
cbar = plt.colorbar(CF,fraction=0.046, pad=0.03,shrink=0.5)
cbar.ax.tick_params(labelsize=10)
cbar.ax.set_ylabel(f'Temperature ($^\circ$C)',fontsize=10);
title = ax.set_title(title_line,fontsize=16)

../../_images/e8c568cd86b51bf3081b45ca8c549deeac00eef8f6a543044f1869b9e4f169e4.png

Create a time series at the gridpoint closest to a given location.¶

siteName = "ALB"
siteLat, siteLon = (42.75, -73.8) # ALB Airport
siteX, siteY = pFull(siteLon, siteLat)
siteXidx, siteYidx = find_closest(x, siteX), find_closest(y, siteY)

Subset the full grid, and convert the units. For the time series, we want all times, so we don’t subset the time dimension.

timeSer = var.isel(x = siteXidx, y = siteYidx).isel(vertDict)
timeSer = timeSer.metpy.convert_units('degC')
times = timeSer[timeDim]

fig = plt.figure(figsize=(12,10))
fig.suptitle(f'HRRR {validTime} Initialization',fontweight='bold',fontsize=15)
ax = fig.add_subplot(111)

ax.set_xlabel('Date/Time')
ax.set_ylabel('2 m Temperature ($^\circ$C)')
ax.set_title(f'Time series at gridpoint closest to {siteName}: Latitude {siteLat} , Longitude {siteLon}' )


# Improve on the default ticking

ax.plot(times, timeSer, linewidth=0.8,color='tab:blue',marker='o',linestyle='--')
ax.xaxis.set_major_locator(HourLocator(interval=1))
hoursFmt = DateFormatter('%d/%H')

ax.xaxis.set_major_formatter(hoursFmt)
ax.tick_params(axis='x', rotation=45, labelsize=8)
ax.set_xlim(times.values[0],times.values[-1]);

                                         

../../_images/e4ca2a5f0b8308bd21aca5861ad0e4c0ac6cc70e53048ddc6eeb533b91f4441f.png

Things to try:¶

Try other variables, locations, and regional models
Create a suite of forecast products, looping over forecast hours, model levels, etc.

ATM433/533 Fall 2023

THREDDS-served GRIB Data: HRRR (Regional Grid) Model

Contents