📍 Pick Points#

Extract Model Grid Data at Specific Points with Herbie#

Extract weather model data at point locations | Interpolate gridded data | Sample NWP grids

Herbie’s pick_points xarray accessor extracts values from weather and climate model grids at specific latitude/longitude locations. Whether you need to sample HRRR, GFS, or other numerical weather prediction (NWP) model output at station locations, city coordinates, or arbitrary points, this tool provides efficient nearest-neighbor and distance-weighted interpolation.

What you’ll learn: Extract and Interpolate Model Data#

In this tutorial, you’ll learn how to:

Extract point values from gridded model data (HRRR, GFS, NAM, etc.)
Sample weather models at specific latitude/longitude coordinates
Work with both regular grids (lat/lon) and curvilinear grids (Lambert Conformal, polar stereographic)
Use nearest-neighbor vs distance-weighted interpolation methods
Create time series from model forecasts at station locations
Extract vertical profiles (model soundings) at points
Interpolate GRIB2 data to observation sites
Optimize performance for extracting hundreds or thousands of points
Cache spatial indices (BallTree) for faster repeated queries

Common Use Cases#

Compare model output to weather station observations
Extract forecast data for specific cities or locations
Create model soundings at radiosonde sites
Sample gridded forecasts along flight paths or trajectories
Validate NWP models against point measurements
Extract model data for machine learning training datasets

Requirements: This tutorial requires the scikit-learn package for spatial indexing:

pip install scikit-learn
# or install with all Herbie extras
pip install 'herbie-data[extras]'

How it Works: Spatial Indexing with BallTree#

Herbie uses scikit-learn’s BallTree algorithm with the haversine formula for accurate nearest-neighbor queries on Earth’s spherical surface. This approach works for any model grid projection without requiring coordinate transformations.

Advantages over other methods:#

✅ Works with curvilinear/projected grids (HRRR, NAM, WRF)
✅ Handles longitude convention differences (0-360° vs ±180°)
✅ Returns distance to nearest grid points (great-circle distance in km)
✅ Supports k-nearest neighbors for interpolation and uncertainty quantification
✅ Fast queries even for thousands of points
✅ No coordinate transformation required

Background: When to Use Pick Points#

For regular latitude-longitude grids (e.g., GFS, IFS), you could use xarray’s advanced indexing to select your points of interest. However, pick_points is necessary when:

Curvilinear grids: Models like HRRR use Lambert Conformal or other map projections
Longitude conventions differ: Model uses [0, 360) but your points are in [-180, 180)
Distance matters: You need to know how far the nearest grid point is from your location
Multiple neighbors: You want k-nearest points for interpolation or statistics
Performance: Extracting many points at once (spatial indexing is very fast)
Batch processing: Repeatedly querying the same grid (BallTree caching)

Development History: Why I Built This

Picking values at nearest-neighbor points in a curvilinear grid has been one of my longstanding challenges. In November 2019, I asked How to select the nearest lat/lon location with multi-dimension coordinates on Stack Overflow, which has had over 28k views. My initial solution found the minimum distance between points and grid cells, but it didn’t scale well for many points.

I iterated through several approaches:

Previous approach (deprecated): ds.herbie.nearest_points used MetPy’s assign_y_x to transform coordinates to the model’s projection (e.g., Lambert Conformal for HRRR), then selected nearest neighbors. This worked but required Cartopy coordinate transformations, and not all datasets have sufficient projection metadata.

Current approach: Using BallTree with haversine distance:

✅ No coordinate transformation needed
✅ Can extract k-nearest neighbors
✅ Returns actual distances (km) using haversine formula
✅ Supports inverse-distance weighted interpolation
✅ Much faster for bulk queries

I later discovered the xoak package, which also uses BallTree for nearest neighbors.

Getting Started#

To use Herbie’s xarray accessors, just import Herbie (the accessor registers automatically):

[1]:

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

import herbie
from herbie import FastHerbie, Herbie

A Very Simple Demonstration#

First, I’ll show how Herbie extracts nearest points from a simple grid. Let’s make a 3x3 xarray Dataset with latitude/longitude coordinates and then extract data nearest two points of interest.

[2]:

# Create a 3x3 xarray dataset
ds = xr.Dataset(
    {"a": (["latitude", "longitude"], [[0, 1, 2], [0, 1, 0], [0, 0, 0]])},
    coords={
        "latitude": (["latitude"], [44, 45, 46]),
        "longitude": (["longitude"], [-99, -100, -101]),
    },
)
ds

The points you want to extract must be given as a Pandas DataFrame with columns latitude and longitude given in degrees.

[3]:

# We want to pick data closest to two points
points = pd.DataFrame(
    {
        "longitude": [-100.25, -99.4],
        "latitude": [44.25, 45.4],
    }
)
points

[3]:

	longitude	latitude
0	-100.25	44.25
1	-99.40	45.40

Now we can pick the nearest grid points with the custom Herbie xarray accessor.

Note method='nearest' is the default behavior.

[4]:

# Pick the value nearest the requested points
matched = ds.herbie.pick_points(points, method="nearest")
matched

WARNING: BallTree caching disabled - tree_name not specified.
         Provide tree_name parameter to enable caching.
INFO: 🌱 Growing new BallTree...🌳 Complete in 0.01s

Alternatively, we can get the inverse-distance weighted mean of the 4 nearest points.

[5]:

# Distance weighted mean
matched_w = ds.herbie.pick_points(points, method="weighted")
matched_w

WARNING: BallTree caching disabled - tree_name not specified.
         Provide tree_name parameter to enable caching.
INFO: 🌱 Growing new BallTree...🌳 Complete in 0.00s

Understanding the Output#

The returned Dataset includes several useful coordinates:

point_latitude, point_longitude: Your requested coordinates
point_grid_distance: Distance (km) to the nearest grid point
latitude, longitude: The actual grid point coordinates
Any additional columns from your points DataFrame (e.g., point_stid)

The point dimension indexes your requested points.

That’s a lot of info to digest. Let’s visualize what we have done with the following figure:

[6]:

fig, axes = plt.subplots(1, 2, figsize=[13, 5])
for p, ax in zip(matched_w.point, axes):
    zz = matched.sel(point=p)

    # Plot grid
    ax.pcolormesh(
        ds.longitude, ds.latitude, ds.a, edgecolor="1", lw=1, cmap="Pastel2_r"
    )
    x, y = np.meshgrid(ds.longitude, ds.latitude)
    x = x.flatten()
    y = y.flatten()
    ax.scatter(x, y, facecolor="tab:red", edgecolor="tab:red", s=100, zorder=100)

    # Plot requested point
    ax.scatter(
        zz.point_longitude,
        zz.point_latitude,
        color="yellow",
        ec="k",
        marker="p",
        s=300,
        zorder=100,
    )
    for i in ds.latitude:
        for j in ds.longitude:
            z = ds.sel(latitude=i, longitude=j)
            ax.text(j, i, f"\na={z.a.item()}", ha="center", va="top")

    # Plot path to nearest point and distance
    # for i, j in zip(x, y):
    #    ax.plot([i, point.longitude.values[0]], [j, point.latitude.values[0]])
    for i in matched_w.k:
        if i == 0:
            kwargs = dict(lw=2, color="k", ls="--")
        else:
            kwargs = dict(color=".5", ls="--")
        z = matched_w.sel(k=i, point=p)
        ax.plot(
            [z.longitude, z.point_longitude], [z.latitude, z.point_latitude], **kwargs
        )
        ax.text(
            z.longitude,
            z.latitude,
            f"  {z.point_grid_distance.item():.1f} km",
            va="center",
            fontweight="bold",
        )
    ax.set_xticks(ds.longitude)
    ax.set_yticks(ds.latitude)
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_title(f"Point={p.item()}", loc="left", fontsize="xx-large")
    ax.set_title(
        f"nearest value:  $a={zz.a.item():.2f}$\nweighted mean:  $a={z.a.item():.2f}$",
        loc="right",
    )

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_11_0.png

The yellow pentagon is the requested point. The grid’s nearest-neighbor point is connected by a thick dashed line. The three other nearest neighbors used to compute the distance-weighted mean are connected by a thin dashed line.

Summary: From a simple 2D Dataset with latitude and longitude coordinates, we got the nearest value for two points of interest. We also compute the inverse-distance weighted mean from the four grids nearest our point of interest.

Pick points from HRRR data#

The above is nice, but you might say, “Yeah, but you can already select data from a grid using xarray advanced selection.” That is true, but it does not work for models with curvilienar grids, like the HRRR model.

Here is a demonstration using real HRRR data.

[7]:

H = Herbie("2024-03-01", model="hrrr")
ds = H.xarray(":(?:TMP|DPT):2 m")
ds

✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Mar-01 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local

Let’s extract the data from some points with some fruit-themed station id names.

[8]:

points = pd.DataFrame(
    {
        "latitude": np.linspace(44, 45.5, 5),
        "longitude": np.linspace(-100, -101, 5),
        "stid": ["McIntosh", "Golden", "Fuji", "Gala", "Honeycrisp"],
    }
)
points

[8]:

	latitude	longitude	stid
0	44.000	-100.00	McIntosh
1	44.375	-100.25	Golden
2	44.750	-100.50	Fuji
3	45.125	-100.75	Gala
4	45.500	-101.00	Honeycrisp

[9]:

%%time
matched = ds.herbie.pick_points(points)
matched

CPU times: user 330 ms, sys: 414 ms, total: 744 ms
Wall time: 764 ms

Notice that a BallTree object for the HRRR model was saved with the name <model_name>_<x_dim>_<y_dim>.pkl in the directory you save Herbie data. This is done so it can be loaded more quickly next time you extract points from this grid.

In the next cell I run the same command, and it uses the cached tree instead.

[10]:

%%time
matched = ds.herbie.pick_points(points)
matched

CPU times: user 523 ms, sys: 241 ms, total: 764 ms
Wall time: 760 ms

You may want to swap the dimensions to use the point_stid coordinate as the dimension instead. Doing this makes it possible to select data by station ID instead of point index.

[11]:

matched = matched.swap_dims({"point": "point_stid"})
matched.sel(point_stid="Honeycrisp")

Now let’s plot each point on a map with the value at the nearest grid point.

[12]:

from herbie.toolbox import EasyMap, ccrs, pc

ax = EasyMap(crs=ds.herbie.crs).ax
ax.pcolormesh(
    ds.longitude,
    ds.latitude,
    ds.t2m,
    cmap="Spectral_r",
    vmax=290,
    vmin=270,
    transform=pc,
)

for i in matched.point_stid:
    z = matched.sel(point_stid=i)
    ax.scatter(z.longitude, z.latitude, color="k", transform=pc)
    ax.scatter(
        z.point_longitude, z.point_latitude, color="yellow", marker=".", transform=pc
    )
    ax.text(
        z.point_longitude,
        z.point_latitude,
        f" {z.point_stid.item()}\n {z.t2m.item() - 273.15:.2f} C",
        transform=pc,
    )
ax.set_extent([-101, -99, 44, 46], crs=pc)
ax.EasyMap.INSET_GLOBE()
ax.adjust_extent()

[12]:

(-289442.889311017, -108684.3334160587, 605320.486730601, 849853.0858732163)

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_23_1.png

Here is another example with several cities

[13]:

points = pd.DataFrame(
    {
        "latitude": [40.7608, 41.8781, 39.7392, 44.9778, 47.6062],
        "longitude": [-111.8910, -87.6298, -104.9903, -93.2650, -122.3321],
        "stid": ["Salt Lake City", "Chicago", "Denver", "Minneapolis", "Seattle"],
        "elevation": [1288, 182, 1655, 264, 56],  # Optional: add metadata
    }
)

# ---------------
# Pick the points
matched = ds.herbie.pick_points(points)

# -------------
# Make the plot

from herbie.paint import NWSTemperature

ax = EasyMap(crs=ds.herbie.crs).STATES().ax
ax.pcolormesh(
    ds.longitude,
    ds.latitude,
    ds.t2m - 273.15,
    **NWSTemperature.kwargs,
    transform=pc,
)

for i in matched.point:
    z = matched.sel(point=i)
    ax.scatter(z.longitude, z.latitude, color="k", transform=pc)
    ax.scatter(
        z.point_longitude, z.point_latitude, color="yellow", marker=".", transform=pc
    )
    ax.text(
        z.point_longitude,
        z.point_latitude,
        f" {z.point_stid.item()}\n {z.t2m.item() - 273.15:.2f} C",
        transform=pc,
    )

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_25_0.png

Pick Points: Timeseries#

[17]:

# TODO: Could demonstrate using FastHerbie here, if the kernel wouldn't crash (memory??)

i = []
for date in pd.date_range("2024-01-01", periods=9, freq="3h"):
    ds = Herbie(date).xarray("(?:DPT|TMP):2 m", remove_grib=False)
    i.append(
        ds.herbie.pick_points(
            pd.DataFrame(
                {
                    "latitude": [40.77069],
                    "longitude": [-111.96503],
                    "stid": ["KSLC"],
                }
            )
        )
    )
slc_ts = xr.concat(i, dim="valid_time")

slc_ts.t2m.plot(x="valid_time", marker=".", color="tab:red")
slc_ts.d2m.plot(x="valid_time", marker=".", color="tab:green")

✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 03:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 06:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 09:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 12:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 15:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 18:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 21:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-02 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local

[17]:

[<matplotlib.lines.Line2D at 0x7d6e21642490>]

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_31_2.png

Benchmark#

Let’s benchmark picking different numbers of points to understand performance scaling when harvesting many points.

[18]:

import warnings

import matplotlib.pyplot as plt
import pandas as pd

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

[19]:

H = Herbie("2024-03-28 00:00", model="hrrr")
ds = H.xarray(r"TMP:\d* mb", remove_grib=False)

✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Mar-28 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local

Using the model’s own grid, I will generate 100 random samples to extract.

[20]:

def generate_test_points(ds, n):
    """Generate n random points from the model grid."""
    import random
    import string

    points = (
        ds[["latitude", "longitude"]]
        .to_dataframe()[["latitude", "longitude"]]
        .sample(n)
        .reset_index(drop=True)
    )
    points["stid"] = [
        "".join(random.choices(string.ascii_letters, k=8)) for _ in range(n)
    ]
    return points


points_self_100 = generate_test_points(ds, 100)
points_self_100

[20]:

	latitude	longitude	stid
0	29.718029	274.785946	SkNxqxCY
1	26.961634	255.333244	DONsGClt
2	47.380357	258.360720	zBYYWqvk
3	48.499379	236.466228	PNKRmpQX
4	37.763784	272.271894	kgiLfUwp
...	...	...	...
95	34.706803	275.869335	RrdbpPvj
96	25.423777	256.984473	ELGWpzoF
97	21.841066	285.435707	iSCwEgju
98	35.878342	285.655062	UQBJvXzT
99	45.136663	294.870079	KYvNRgCr

100 rows × 3 columns

[21]:

%%timeit
y1 = ds.herbie.pick_points(points_self_100)

542 ms ± 144 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[22]:

%%timeit
y2 = ds.herbie.pick_points(points_self_100, method="weighted")

1.84 s ± 234 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[23]:

%%timeit
# Try the deprecated `nearest_points` method which used MetPy
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    y3 = ds.herbie.nearest_points(points_self_100)

541 ms ± 76.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[24]:

# Do I get the same answer for the old and new method?
y1 = ds.herbie.pick_points(points_self_100)
y2 = ds.herbie.pick_points(points_self_100, method="weighted")
y3 = ds.herbie.nearest_points(points_self_100)

all(y1.latitude == y3.latitude), all(y1.longitude == y3.longitude)

/home/blaylock/GITHUB/Herbie/src/herbie/nearest_points.py:109: UserWarning: More than one time coordinate present for variable  "t".
  ds = ds.metpy.assign_y_x()

[24]:

(True, True)

[25]:

%%time
ax = EasyMap(crs=ds.herbie.crs).ax

ax.pcolormesh(
    ds.longitude,
    ds.latitude,
    ds.isel(isobaricInhPa=0).t,
    cmap="Spectral_r",
    transform=pc,
)

ax.scatter(y1.longitude, y1.latitude, color="k", transform=pc)
ax.scatter(
    y1.point_longitude, y1.point_latitude, color="yellow", marker=".", transform=pc
)

ax.EasyMap.INSET_GLOBE()

CPU times: user 1.48 s, sys: 257 ms, total: 1.73 s
Wall time: 2.6 s

[25]:

<GeoAxes: label='inset_axes'>

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_41_2.png

Not bad. Less than half a second to extract 100 points from the HRRR grid.

Now let’s try even more…

[26]:

%%time

n_samples = (
    list(range(1, 10, 2)) + list(range(10, 100, 20)) + list(range(100, 1_000, 200))
    # + list(range(1_000, 10_000, 2_000))
    # + [100_000, 1_000_000]
)

times_pick_points_nearest = []
times_pick_points_weighted = []
times_nearest_points = []

with warnings.catch_warnings():
    warnings.simplefilter("ignore")

    for n in n_samples:
        points = generate_test_points(ds, n)

        timer = pd.Timestamp("now")
        _ = ds.herbie.nearest_points(points)
        times_nearest_points.append((pd.Timestamp("now") - timer).total_seconds())

        timer = pd.Timestamp("now")
        _ = ds.herbie.pick_points(points, method="nearest")
        times_pick_points_nearest.append((pd.Timestamp("now") - timer).total_seconds())

        timer = pd.Timestamp("now")
        _ = ds.herbie.pick_points(points, method="weighted")
        times_pick_points_weighted.append((pd.Timestamp("now") - timer).total_seconds())

        print(f"✓ Completed {n:,} points")

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(
    n_samples,
    times_pick_points_nearest,
    marker=".",
    lw=2,
    color="tab:blue",
    label="pick_points (nearest)",
)
ax.plot(
    n_samples,
    times_pick_points_weighted,
    marker=".",
    lw=2,
    color="tab:orange",
    label="pick_points (weighted)",
)
ax.plot(
    n_samples,
    times_nearest_points,
    marker=".",
    lw=1,
    color="k",
    ls="--",
    label="nearest_points (deprecated)",
)

ax.set_xscale("log")
ax.set_title("Performance: Herbie methods to get data at a point")
ax.set_xlabel("Number of points picked")
ax.set_ylabel("Query Time (seconds)")
ax.set_ylim(ymin=0)
ax.set_xlim(xmin=1)
ax.grid(True, alpha=0.3)
ax.legend(loc="upper center", bbox_to_anchor=[0.5, -0.2])

✓ Completed 1 points
✓ Completed 3 points
✓ Completed 5 points
✓ Completed 7 points
✓ Completed 9 points
✓ Completed 10 points
✓ Completed 30 points
✓ Completed 50 points
✓ Completed 70 points
✓ Completed 90 points
✓ Completed 100 points
✓ Completed 300 points
✓ Completed 500 points
✓ Completed 700 points
✓ Completed 900 points
CPU times: user 1min, sys: 20.9 s, total: 1min 21s
Wall time: 1min 23s

[26]:

<matplotlib.legend.Legend at 0x7d6e2157c2d0>

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_43_2.png

Nice! This BallTree method used by pick_points() scales well for harvesting many points 😎

Lets see how well this works for the GFS model…

[27]:

gfs = Herbie("2024-01-01", model="gfs").xarray(":TMP:2 m ")

✅ Found ┊ model=gfs ┊ product=pgrb2.0p25 ┊ 2024-Jan-01 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ local

[28]:

%%time
gfs.herbie.pick_points(points_self)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
File <timed eval>:1

NameError: name 'points_self' is not defined