Plucking nearest lat/lon point on a curvilinear grid

[nicolesplaining]

I am working with the HRRR data, which I loaded into xarray through a grib2 file. The file can be downloaded here.

<xarray.Dataset>
Dimensions:     (hybrid: 50, y: 1059, x: 1799)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
  * hybrid      (hybrid) float64 1.0 2.0 3.0 4.0 5.0 ... 47.0 48.0 49.0 50.0
    latitude    (y, x) float64 ...
    longitude   (y, x) float64 ...
    valid_time  datetime64[ns] ...
Dimensions without coordinates: y, x
Data variables: (12/17)
    pres        (hybrid, y, x) float32 ...
    clwmr       (hybrid, y, x) float32 ...
    unknown     (hybrid, y, x) float32 ...
    rwmr        (hybrid, y, x) float32 ...
    snmr        (hybrid, y, x) float32 ...
    grle        (hybrid, y, x) float32 ...
    ...          ...
    q           (hybrid, y, x) float32 ...
    u           (hybrid, y, x) float32 ...
    v           (hybrid, y, x) float32 ...
    w           (hybrid, y, x) float32 ...
    tke         (hybrid, y, x) float32 ...
    mdens       (hybrid, y, x) float32 ...

The latitude and longitude are in terms of y and x and is in the lambert conformal conic projection with the specifications:

projection = ccrs.LambertConformal(central_longitude=262.5, 
                                   central_latitude=38.5, 
                                   standard_parallels=(38.5, 38.5),
                                    globe=ccrs.Globe(semimajor_axis=6371229,
                                                     semiminor_axis=6371229))

I would like to get the data at the point nearest to a specific lat/lon point, say (40.76, -111.876183) - which is in the regular cartesian coordinates. I can't use xarray advanced selection because the data is in terms of y and x, not lat/lon. How do I approach this?

[kmuehlbauer]

Here in SO is a topic discussing several solutions. HTH.
https://stackoverflow.com/questions/58758480/xarray-select-nearest-lat-lon-with-multi-dimension-coordinates

[benbovy]

A bit late but you may look at https://xoak.readthedocs.io/en/latest/, which does specifically what you want.

[marceloandrioni]

Hello, I like to use pyinterp to do a RTree search, e.g.:

import numpy as np
import xarray as xr
import pyinterp


def get_indexes(*,
                longitude: float,
                latitude: float,
                lon2d: np.ndarray,
                lat2d: np.ndarray,
                ) -> tuple[int, int, float]:
    """Given a longitude and latitude, return the row and column indexes and
    distance (in meters) for the nearest grid point in the
    longitude[row,column] and latitude[row, column] matrices."""

    lon1d = lon2d.flatten()
    lat1d = lat2d.flatten()

    idxs1d = np.arange(lon1d.size)
    idxs = idxs1d.reshape(lon2d.shape)

    # rtree search
    mesh = pyinterp.RTree()
    mesh.packing(np.column_stack((lon1d, lat1d)), idxs1d)

    coordinates = np.atleast_2d([longitude, latitude])

    distances, values = mesh.query(coordinates)

    distance, value = distances[0][0], int(values[0][0])

    if distance == -1:
        raise ValueError("Could not find a nearest neighbor")

    # get row and columns indexes
    irow, icolumn = np.r_[np.where(idxs == value)]

    return irow, icolumn, distance


infile = "hrrr.t23z.wrfnatf00.grib2"

# requested coordinates
req_lat = 30.25
req_lon = 270.83

ds = xr.open_dataset(infile, filter_by_keys={'typeOfLevel': 'hybrid'})

iy, ix, distance = get_indexes(longitude=req_lon,
                               latitude=req_lat,
                               lon2d=ds["longitude"].values,
                               lat2d=ds["latitude"].values)

# nearest grid point
near_lat = ds["latitude"].sel(y=iy, x=ix).values
near_lon = ds["longitude"].sel(y=iy, x=ix).values

print(f"Requested location: lat={req_lat}, lon={req_lon}")
print(f"Nearest grid point: lat={near_lat}, lon={near_lon}")
print(f"Distance: {distance} meters")

print(ds["pres"].sel(y=iy, x=ix).values)

Requested location: lat=30.25, lon=270.83
Nearest grid point: lat=30.23842555308348, lon=270.82210024377054
Distance: 1491.4686944201699 meters
[102069.55   101765.95   101208.37   100245.04    98823.94    96995.484
  94808.195   92415.79    89814.87    86904.46    83685.24    80159.61
  76334.25    72204.95    67768.375   63026.43    58037.113   53212.656
  48881.477   45015.67    41547.617   38415.348   35606.957   33085.613
  30819.87    28788.402   26965.203   25274.863   23751.11    22358.748
  20968.7     19579.754   18190.89    16802.025   15413.166   14024.307
  12635.452   11246.598    9857.746    8567.393    7533.1406   6690.9634
   5907.8867   5174.0596   4494.408    3864.0068   3273.0054   2721.4038
   2209.2021   1731.4756]