A common usage of xarray is to retrieve climate model data from the grid cells closest to a weather station. That might look like this

``` import xarray as xr import numpy as np

ds = xr.tutorial.open_dataset("air_temperature")

Define taget latitude and longitude

tgt_lat = np.linspace(0, 100, num=10) tgt_lon = np.linspace(0, 100, num=10)

Retrieve data at target latitude and longitude

tgt_data = ds['air'].sel(lon=tgt_lon, lat=tgt_lat, method='nearest') ```

My problem is that I am trying subset ds to 10 points in space (which is the length of tgt_lat and tgt_lon), but in fact xarray retrieves 100 points (10 latitude by 10 longitude). I can get around this by calling tgt_data = tgt_data.values.diagonal(). But this results in a non-xarray object. Furthermore, if instead of querying for 10 points in space, I query for 10,000, I run out of memory because xarray retrieves 100,000,000 points in space (10,000^2).

Is there a way to only retrieve the diagonal elements? If not, is this something that should be added?