📍 Pick Points#

Herbie’s pick_points xarray accessor picks data from a model grid nearest locations of interest. Think of it like picking apples from a tree. Herbie’s implementation uses the scikit-learn BallTree algorithm with the haversine formula to get nearest-neighbor values. The haversine formula is important when working with latitude and longitude coordinates.

For regular latitude-longitude grids (e.g., GFS, IFS), you could use xarray’s advanced indexing instead to easily select your points of interest, but that solution doesn’t work in the following cases:

Models with curvilinear grids (i.e. not a regular latitude-longitude grid);
Longitude units, [0, 360) vs [-180, 180), differ between model and requested point.
The distance to the nearest neighbor needs to be known.
Picking more than one nearest neighbor points.

Picking values at nearest-neighbor points in a curvilinear grid has been one of my longstanding questions. 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 in four years. I suggested a rudimentary solution that determined the nearest neighbor point by finding the absolute minimum value between the point and grid difference, but this solution didn’t scale well for picking many points.

I iterated over different solutions since then.

In the deprecated Herbie accessor ds.herbie.nearest_points, I used MetPy’s assign_y_x, transformed the latitude and longitude points to the coordinate reference system coordinates (e.g. lambert conformal for HRRR), and then picked the nearest neighbor points. That worked fine, but it required a coordinate transformation by Cartopy, and not all xarray Datasets have sufficient information to determine the appropriate coordinate transformation.

This new approach using BallTree, we don’t need to do a coordinate transformation, we can pick the \(k\)-th nearest neighbors, and the distance to the nearest neighbors (haversine formula) is returned. I have also implemented the capability to compute the inverse-distance weighted mean of the four (or more) grid points nearest your points of interest.

As I write this, I learned about the xoak package, which also uses BallTree to query nearest neighbors.

Let’s get started…

To use Herbie’s xarray accessors, you just need to import anything from Herbie.

[1]:

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

import warnings

import herbie
from herbie import Herbie, FastHerbie

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: Herbie won't cache the BallTree because it
         doesn't know what to name it. Please specify
         `tree_name` to cache the tree for use later.
INFO: 🌱 Growing new BallTree...🌳 BallTree grew in 0.007 seconds.

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: Herbie won't cache the BallTree because it
         doesn't know what to name it. Please specify
         `tree_name` to cache the tree for use later.
INFO: 🌱 Growing new BallTree...🌳 BallTree grew in 0.0065 seconds.

That’s a lot of info to digest. I’ll help you 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 @ aws

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

INFO: 🌱 Growing new BallTree...🌳 BallTree grew in 3.3 seconds.
INFO: Saved BallTree to /home/blaylock/data/BallTree/hrrr_1799-1059.pkl
CPU times: user 3.36 s, sys: 210 ms, total: 3.57 s
Wall time: 3.55 s

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 132 ms, sys: 20.7 ms, total: 153 ms
Wall time: 150 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 toolbox import EasyMap, pc, ccrs

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

Pick Points: Timeseries#

[16]:

# 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")
    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 @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 03:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 06:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 09:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 12:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 15:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 18:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-01 21:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
✅ Found ┊ model=hrrr ┊ product=sfc ┊ 2024-Jan-02 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
👨🏻‍🏭 Created directory: [/home/blaylock/data/hrrr/20240102]

[16]:

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

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_29_2.png

Benchmark#

Let’s see how fast ds.herbie.pick_points is for harvesting many points…

[17]:

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 @ aws

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

[18]:

def generate_random_string(len=8):
    """Generate a random string."""
    import random
    import string

    return "".join(random.choices(string.ascii_letters + string.digits, k=len))


n = 100
points_self = (
    ds[["latitude", "longitude"]]
    .to_dataframe()[["latitude", "longitude"]]
    .sample(n)
    .reset_index(drop=True)
)
points_self["stid"] = [generate_random_string() for _ in range(n)]
points_self

[18]:

	latitude	longitude	stid
0	28.118041	251.508694	A4Bt4ArI
1	45.433523	236.549102	M8TVLmnW
2	29.326794	290.301011	EFnRTuds
3	48.066147	293.088237	LGzLiHaU
4	32.020415	247.694797	tzABimIV
...	...	...	...
95	44.124212	279.341471	zA7MhOEm
96	44.208684	270.780126	OSrB5uVX
97	30.286338	254.814174	e5xYHMth
98	38.231154	256.410611	ZwVijXzk
99	42.481527	268.409357	PsiEP6f0

100 rows × 3 columns

[19]:

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

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

[20]:

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

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

[21]:

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

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

[22]:

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

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

/tmp/ipykernel_1317/3015965854.py:4: DeprecationWarning: The accessor `ds.herbie.nearest_points` is deprecated in favor of the `ds.herbie.pick_points` which uses the BallTree algorithm instead.
  y3 = ds.herbie.nearest_points(points_self)
/home/blaylock/GITHUB/Herbie/herbie/accessors.py:524: UserWarning: More than one time coordinate present for variable  "t".
  crs = ds.metpy_crs.item().to_cartopy()

[22]:

(True, True)

[23]:

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()

[23]:

<GeoAxes: >

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_38_1.png

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

Now let’s try even more…

[32]:

n_samples = (
    [1, 2, 3, 4, 5, 6, 7, 8, 9]
    + list(range(10, 100, 10))
    + list(range(100, 1_000, 100))
    + list(range(1_000, 10_000, 1_000))
    + [100_000, 1_000_000]
)

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

    # Timers for nearest_points (deprecated)
    times_np = []
    samples_np = []

    # Timers for pick_points method='nearest'
    times_pp = []
    samples_pp = []

    # Timers for pick_points method='weighted'
    times_w = []
    samples_w = []

    for s in n_samples:
        samples_np.append(s)
        timer = pd.Timestamp("now")
        y1 = ds.herbie.nearest_points(points_self)
        times_np.append((pd.Timestamp("now") - timer).total_seconds())

        samples_pp.append(s)
        timer = pd.Timestamp("now")
        y1 = ds.herbie.pick_points(points_self, method="nearest")
        times_pp.append((pd.Timestamp("now") - timer).total_seconds())

        samples_w.append(s)
        timer = pd.Timestamp("now")
        y1 = ds.herbie.pick_points(points_self, method="weighted")
        times_w.append((pd.Timestamp("now") - timer).total_seconds())


plt.figure(figsize=[8, 4])
plt.plot(
    samples_pp,
    times_pp,
    marker=".",
    lw=2,
    color="tab:blue",
    label="pick_points (nearest)",
)
plt.plot(
    samples_w,
    times_w,
    marker=".",
    lw=2,
    color="tab:orange",
    label="pick_points (weighted)",
)
plt.plot(
    samples_np,
    times_np,
    marker=".",
    lw=1,
    color="k",
    ls="--",
    label="nearest_points (deprecated)",
)

plt.title("Herbie methods to get data at a point")
plt.xlabel("Number of points picked")
plt.ylabel("Query Time (s)")
plt.ylim(ymin=0)
plt.xlim(xmin=1)
plt.legend(loc="upper center", bbox_to_anchor=[0.5, -0.2])
plt.gca().set_xscale("log")

../../../_images/user_guide_tutorial_accessor_notebooks_pick_points_40_0.png

Nice! This method scales well for harvesting many points 😎

Lets see how well this works for the GFS model…

[25]:

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 @ aws

[26]:

%%time
gfs.herbie.pick_points(points_self)

INFO: 🌱 Growing new BallTree...🌳 BallTree grew in 1.7 seconds.
INFO: Saved BallTree to /home/blaylock/data/BallTree/gfs_1440-721.pkl
CPU times: user 1.8 s, sys: 100 ms, total: 1.9 s
Wall time: 1.89 s

📍 Pick Points#

A Very Simple Demonstration#

Pick points from HRRR data#

Pick Points: Sounding#

Pick Points: Timeseries#

Benchmark#

Advanced Options#