Nearest points error with ECMWF forecasts

[svarlamov]

I'm running into an error when I try to run the nearest_points accessor on ECMWF forecasts. It appears that the error is related to the coordinate system, but I don't have enough background knowledge to fully understand why this is failing even after spending a few hours trying to debug this myself. This code reliably reproduces the error:

from herbie import Herbie
H = Herbie("2022-01-26", model="ecmwf", product="oper", fxx=12)
ds = H.xarray(":2t:")
dsi = ds.herbie.nearest_points(points=(-100, 40))

Error:

~/.pyenv/versions/my-env/lib/python3.11/site-packages/herbie/accessors.py:260: UserWarning: More than one time coordinate present for variable  "gribfile_projection".
  ds = ds.metpy.parse_cf()
~/.pyenv/versions/my-env/lib/python3.11/site-packages/herbie/accessors.py:264: UserWarning: More than one time coordinate present for variable  "t2m".
  ds = ds.metpy.assign_y_x()
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[4], line 1
----> 1 dsi = ds.herbie.nearest_points(points=(-100, 40))
      2 dsi

File ~/.pyenv/versions/my-env/lib/python3.11/site-packages/herbie/accessors.py:264, in HerbieAccessor.nearest_points(self, points, names, verbose)
    260     ds = ds.metpy.parse_cf()
    262 # Apply the MetPy method `assign_y_x` to the dataset
    263 # https://unidata.github.io/MetPy/latest/api/generated/metpy.xarray.html?highlight=assign_y_x#metpy.xarray.MetPyDataArrayAccessor.assign_y_x
--> 264 ds = ds.metpy.assign_y_x()
    266 # Convert the requested [(lon,lat), (lon,lat)] points to map projection.
    267 # Accept a list of point tuples, or Shapely Points object.
    268 # We want to index the dataset at a single point.
    269 # We can do this by transforming a lat/lon point to the grid location
    270 crs = ds.metpy_crs.item().to_cartopy()

File ~/.pyenv/versions/my-env/lib/python3.11/site-packages/metpy/xarray.py:1032, in MetPyDatasetAccessor.assign_y_x(self, force, tolerance)
   1029             grid_prototype = data_var
   1030         if (not force and (hasattr(data_var.metpy, 'y')
   1031                            or hasattr(data_var.metpy, 'x'))):
-> 1032             raise RuntimeError('y/x coordinate(s) are present. If you wish to '
   1033                                'overwrite these, specify force=True.')
   1035 # Calculate y and x from grid_prototype, if it exists, and assign
   1036 if grid_prototype is None:

RuntimeError: y/x coordinate(s) are present. If you wish to overwrite these, specify force=True.

Versions:

• MetPy: 1.6.1
• herbie-data: 2023.12.4
• Python: 3.11.4

This is possibly related to #179

Thank you so much for any thoughts on how to fix or perhaps work around this error!

[blaylockbk]

Thanks for reporting this! I can't promise I will have time to dig into this soon, but I can reproduce the error. Looks like something in how Herbie is using MetPy that is different for this model than other models like HRRR.

FYI, There is no error for the HRRR model

from herbie import Herbie

H = Herbie("2022-01-26", model="HRRR", fxx=12)
ds = H.xarray("TMP:2 m")
dsi = ds.herbie.nearest_points(points=(-100, 40))

[blaylockbk]

I rewrote this xarray accessor and named it pick_points, which resolves this issue.

In Herbie version 2024.5.0, you can now do this to get value at two different points...

from herbie import Herbie
import pandas as pd

H = Herbie("2022-01-26", model="ifs", product="oper", fxx=12)
ds = H.xarray(":2t:")

points = pd.DataFrame({"latitude": [40, 41], "longitude": [-100, 101]})
dsi = ds.herbie.pick_points(points)
dsi