I often have xarray datasets with latitude and longitude coordinates (not as xarray dimensions). When I need to pluck values from a point nearest to a specific lat/lon coordinate (like a weather station), this is what I use...

nearest_latlon_index.py

import numpy as np
import xarray as xr

def to_180(lon):
    """
    Wrap longitude from degrees [0, 360] to degrees [-180, 180].

    An alternative method is
        lon[lon>180] -= 360

    Parameters
    ----------
    lon : array_like
        Longitude values
    """
    lon = np.array(lon)
    lon = (lon + 180) % 360 - 180
    return lon

def nearest_latlon_index(ds, points, return_value=True, verbose=False):
    """
    Find the index for a point nearest a given latitude and longitude. 
    
    See GitHub notebook for a demo
        https://github.com/blaylockbk/pyBKB_v3/blob/master/demo/Nearest_lat-lon_Grid.ipynb
    
    Parameters
    ----------
    ds : xarray.Dataset
        The Dataset should include coordinates for both 'latitude' and 
        'longitude'.
    points : tuple or list of tuples
        The latitude and longitude (lat, lon) tuple for the point you 
        want to pluck from the Grid. Multple tuples may be given as a
        list to return the values from multiple points.
    return_value : bool
        If True, return the value of the Dataset at the indexes.
        If False, return just the index values.
    
    Returns
    -------
    The x and y index for the grid point nearest the point's latitude
    and longitude. The great-circle distance between each point pair is
    given as a Dataset attribute in kilometers.
    """
    if 'lat' in ds:
        ds = ds.rename(dict(lat='latitude', lon='longitude'))
    
    if isinstance(points, tuple):
        # If a tuple is give, turn into a one-item list.
        points = [points]
        
    xs = []         # x index values
    ys = []         # y index values
    distances = []  # distance between the pair of points
    
    for point in points:
        assert len(point) == 2, "``points`` should be a tuple or list of tuples (lat, lon)"
        
        p_lat, p_lon = point

        # Force longitude values to range from -180 to 180 degrees.
        p_lon = to_180(p_lon)
        ds['longitude'][:] = to_180(ds.longitude)

        # Find absolute difference between requested point and the grid coordinates.
        abslat = np.abs(ds.latitude - p_lat)
        abslon = np.abs(ds.longitude - p_lon)

        # Create grid of the maximum values of the two absolute grids
        c = np.maximum(abslon, abslat)

        # Find location where lat/lon minimum absolute value intersects
        y, x = np.where(c == np.min(c))
        xs.append(x[0])
        ys.append(y[0])

        # Matched Grid lat/lon
        g_lat = ds.latitude.isel(x=x, y=y).data[0][0]
        g_lon = ds.longitude.isel(x=x, y=y).data[0][0]
        
        
        # Calculate Great Circle distance between requeted point and 
        # matched grid point. 
        # Based on https://andrew.hedges.name/experiments/haversine/
        
        # approximate radius of earth in km
        R = 6373.0

        lat1 = np.deg2rad(p_lat)
        lon1 = np.deg2rad(p_lon)
        lat2 = np.deg2rad(g_lat)
        lon2 = np.deg2rad(g_lon)

        dlon = lon2 - lon1
        dlat = lat2 - lat1

        a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
        c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))

        distance = R * c
        distances.append(distance)
        
        if verbose: 
            print(f"🔎 Matched requested point ({p_lat:.3f}, {p_lon:.3f}) to grid point ({g_lat:.3f}, {g_lon:.3f})...distance of {distance:.2f} km.")
       
    if return_value:
        ds = ds.isel(x=xs, y=ys)
        ds.attrs['point x,y index'] = list(zip(xs, ys))
        ds.attrs['distances (in km)'] = distances
        return ds
    else:
        return list(zip(xs, ys))