GRIB Plotting

from datetime import datetime
from urllib.request import urlopen

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pygrib
import xarray as xr

# Download the GRIB file from the webserver and save to a local file
fname = 'forecast.bin'
now = datetime.utcnow()
with open(fname, 'wb') as cache_file:
    fobj = urlopen(f'https://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/hrrr.{now:%Y%m%d}/conus/hrrr.t00z.wrfnatf00.grib2')
    cache_file.write(fobj.read())

# Open grib and get the first message
grb = pygrib.open(fname)
msg = grb[1]

# Set up the projection information using information from the GRIB message
globe = ccrs.Globe(semimajor_axis=msg['radius'])
data_proj = ccrs.LambertConformal(globe=globe, central_latitude=msg['LaDInDegrees'],
                                  central_longitude=msg['LoVInDegrees'], standard_parallels=(msg['Latin1InDegrees'], msg['Latin2InDegrees']))

# Transform the corner point lat/lon to an x/y in the projected coordinates
x0, y0 = data_proj.transform_point(msg['longitudeOfFirstGridPointInDegrees'],
                                   msg['latitudeOfFirstGridPointInDegrees'],
                                   src_crs=ccrs.PlateCarree())

# Generate the grid x/y from that corner using the number of point and the spacing
x = x0 + np.arange(msg['Nx']) * msg['DxInMetres']
y = y0 + np.arange(msg['Ny']) * msg['DyInMetres']

# Pull valid time out to use as a coordinate
t = msg.validDate

# Grab lat/lons from pygrib
lats, lons = msg.latlons()

# Create an xarray DataArray from the coordinates and information in the message
data = xr.DataArray(data=msg.values, dims=["y", "x"],
                    coords=dict(x=x, y=y, time=t, latitude=(('y', 'x'), lats), longitude=(('y', 'x'), lons)),
                    attrs=dict(parameter=msg['parameter'], unit=msg['units']))

# Convert a lon/lat into the x/y coordinates
lat_pt, lon_pt = (40, -105)
x_pt, y_pt = data_proj.transform_point(lon_pt, lat_pt, ccrs.PlateCarree())

# Use that point to select
my_value = data.sel(x=x_pt, y=y_pt, method='nearest')

# Plot on a map to verify we got this right
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=data_proj)

# Plot the data
ax.pcolormesh(data.x, data.y, data)

# Plot original point in lon/lat space
ax.plot(lon_pt, lat_pt, 'bs', markersize=20, transform=ccrs.PlateCarree())

# Plot point after transforming manually
ax.plot(x_pt, y_pt, 'ro')

ax.coastlines()