🗺️ Brian’s Maps

Notice: Carpenter Workshop must be installed.

This notebook demonstrates how to plot GRIB2 data retrieved from Herbie on a map using my custom Carpenter Workshop tools.This is not meant to replace an other Cartopy tutorials, so if you are unfamilar with Cartopy, I suggest you first refer to the Project Pythia Cartopy tutorial.

import matplotlib.pyplot as plt
import numpy as np

from herbie import Herbie
from toolbox import EasyMap, pc, ccrs

Download some HRRR data. Let’s look at composite reflectivity

ds = Herbie("2022-12-10 12:00").xarray("REFC:entire")

✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2022-Dec-10 12:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws

---

What does the data look like when we plot it?

ds.refc.plot()

<matplotlib.collections.QuadMesh at 0x19d6cda2470>

Fun fact: this reflectivity field uses -10 as the missing value. Lets just replace any negative number with nan. (Can anyone tell me why there are reflectivities less than zero?)

ds["refc"] = ds.refc.where(ds.refc > 0)  # sets all negative values to nan

ds.refc.plot()

<matplotlib.collections.QuadMesh at 0x19d6f412fe0>

Applying the mask sure makes the reflectivity data look better.

Ok, but we really want to see this data on a map. Let’s use the Cartopy Workshop EasyMap().

ax = EasyMap().ax

Notice that it gives us a Platte Carree map of the world. We can use this to plot our HRRR data onto

ax = EasyMap().ax
ds.refc.plot(
    x="longitude",
    y="latitude",
    ax=ax,
    transform=pc,
    cbar_kwargs={"shrink": 0.5, "orientation": "horizontal", "pad": 0.01},
)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x19d00211150>

We can add other features, like states and a border around the HRRR domain

ax = EasyMap().STATES().OCEAN().LAND().DOMAIN(ds).ax
ds.refc.plot(
    x="longitude",
    y="latitude",
    ax=ax,
    transform=pc,
    cbar_kwargs={"shrink": 0.5, "orientation": "horizontal", "pad": 0.01},
)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x19d00262080>

Dark mode even looks cool.

ax = EasyMap(dark=True).STATES().OCEAN().LAND().DOMAIN(ds).ax
ds.refc.plot(
    x="longitude",
    y="latitude",
    ax=ax,
    transform=pc,
    cbar_kwargs={"shrink": 0.5, "orientation": "horizontal", "pad": 0.01},
)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x19d003fbdc0>

To plot on different map projects, we need to know what coordinate system the model uses. Herbie has a custom xarray accessor that uses Metpy and Pygrib to parse the projection informtion.

crs = ds.herbie.crs

ax = EasyMap("50m", crs=crs).STATES().OCEAN().LAND().DOMAIN(ds).ax
ds.refc.plot(
    x="longitude",
    y="latitude",
    ax=ax,
    transform=pc,
    cbar_kwargs={"shrink": 0.5, "orientation": "horizontal", "pad": 0.01},
)

ax.adjust_extent();  # make the map extent slightly larger

Now lets zoom on a specific region and show a Stamen map background

crs = ds.herbie.crs

ax = EasyMap("50m", crs=crs, dpi=150).STATES().STAMEN(zoom=6).ax
ds.refc.plot(
    x="longitude",
    y="latitude",
    ax=ax,
    transform=pc,
    cmap="gist_ncar",
    vmin=-10,
    vmax=50,
    cbar_kwargs={"shrink": 0.5, "orientation": "horizontal", "pad": 0.01},
)

ax.set_extent([-120, -105, 35, 45], crs=pc)