Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

.sel(...., method='nearest') fails for large requests. #4630

Closed
EricKeenan opened this issue Nov 30, 2020 · 8 comments · Fixed by #4711
Closed

.sel(...., method='nearest') fails for large requests. #4630

EricKeenan opened this issue Nov 30, 2020 · 8 comments · Fixed by #4711

Comments

@EricKeenan
Copy link
Contributor

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?

@dcherian
Copy link
Contributor

You should be able to do this with "vectorized indexing": https://xarray.pydata.org/en/stable/indexing.html#vectorized-indexing

@EricKeenan
Copy link
Contributor Author

@dcherian Thanks for pointing me in the right direction. I'm trying to implement this with vectorized indexing, but it seems that my queries need to exactly match the xarray object lat/lon, which is why I tried method='nearest'. Am I missing something?

@dcherian
Copy link
Contributor

dcherian commented Dec 1, 2020

TO use "vectorized indexing", tgt_lat and tgt_lon need to be DataArrays with a common dimension name that is not a dimenion in ds

import xarray as xr
import numpy as np

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

# Define taget latitude and longitude
tgt_lat = xr.DataArray(np.linspace(0, 100, num=10), dims="points") # <---
tgt_lon = xr.DataArray(np.linspace(0, 100, num=10), dims="points") # <---

# Retrieve data at target latitude and longitude
tgt_data = ds['air'].sel(lon=tgt_lon, lat=tgt_lat, method='nearest')
tgt_data

image

@EricKeenan
Copy link
Contributor Author

👏 👍 I didn't realize I needed to do that. Thanks for letting me know. Problem solved - marking this as closed.

@keewis
Copy link
Collaborator

keewis commented Dec 1, 2020

this trick is not mentioned in the narrative documentation (or rather: I can't find it), and the docstrings of isel and sel don't contain any examples at all.

Since I believe it should be documented somewhere I'm reopening this to make sure we don't forget. Also, we would definitely welcome a PR adding this, if you're up for it.

@keewis keewis reopened this Dec 1, 2020
@EricKeenan
Copy link
Contributor Author

I'd be happy to give this a shot. But I'm not sure where to start... Can you point me to an example PR that has done something similar?

@keewis
Copy link
Collaborator

keewis commented Dec 1, 2020

sure. #4621 added examples for interp / interpolate_na. For the narrative documentation I don't have a example PR but I think it's fine to just extend the Vectorized indexing section (thoughts, @dcherian?).

@EricKeenan
Copy link
Contributor Author

Thanks for sharing! I'll give this a first shot before the end of the year.

EricKeenan added a commit to EricKeenan/xarray that referenced this issue Dec 18, 2020
EricKeenan added a commit to EricKeenan/xarray that referenced this issue Dec 18, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants