xarray interp variable with two dimensional gird into one point

Question

I have a dataset with two-dimension lon/lat and I want to calculate the value of one certain point.

The ncfile can be download from ftp://ftp.awi.de/sea_ice/product/cryosat2_smos/v204/nh/LATEST/ and the variable can be retrieved by the following ways

SAT = xr.open_dataset('W_XX-ESA,SMOS_CS2,NH_25KM_EASE2_20181231_20190106_r_v204_01_l4sit.nc')
SAT_lon,SAT_lat = SAT.lon,SAT.lat
SAT_SIT = SAT.analysis_sea_ice_thickness

The SAT_SIT is presented as

<xarray.DataArray 'analysis_sea_ice_thickness' (time: 1, yc: 432, xc: 432)>
[186624 values with dtype=float64]
Coordinates:
  * time     (time) datetime64[ns] 2019-01-03T12:00:00
  * xc       (xc) float32 -5.388e+03 -5.362e+03 ... 5.362e+03 5.388e+03
  * yc       (yc) float32 5.388e+03 5.362e+03 ... -5.362e+03 -5.388e+03
    lon      (yc, xc) float32 -135.0 -135.1 -135.3 -135.4 ... 44.73 44.87 45.0
    lat      (yc, xc) float32 16.62 16.82 17.02 17.22 ... 17.02 16.82 16.62
Attributes:
    units:          m
    long_name:      CS2SMOS merged sea ice thickness
    standard_name:  sea_ice_thickness
    grid_mapping:   Lambert_Azimuthal_Grid

Now, I want to check the value of SAT_SIT in one certain point, for example, lon=100 lat=80. Is there any potential way to handle this?

Note: the xarray only support interpolation of 1d coordinate, "Currently, our interpolation only works for regular grids. Therefore, similarly to sel(), only 1D coordinates along a dimension can be used as the original coordinate to be interpolated."

There are some possible ways to solve this questions in this link. However, the first method is somehow complex because I need to interp the variable to many points. The second methods comes up with a error.

SAT = xr.open_dataset('W_XXESA,SMOS_CS2,NH_25KM_EASE2_20181231_20190106_r_v204_01_l4sit.nc')
SAT_lon,SAT_lat = SAT.lon,SAT.lat
lon_test,lat_test = -80,80
data_crs = ccrs.LambertAzimuthalEqualArea(central_latitude=90,central_longitude=0,false_easting=0,false_northing=0)
x, y = data_crs.transform_point(lon_test, lat_test, src_crs=ccrs.PlateCarree())
SAT.sel(xc=x,yc=y)

the error is

Traceback (most recent call last):
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3441, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-54-b210190142ea>", line 6, in <module>
    SAT.sel(xc=x,yc=y)
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/xarray/core/dataset.py", line 2475, in sel
    self, indexers=indexers, method=method, tolerance=tolerance
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/xarray/core/coordinates.py", line 422, in remap_label_indexers
    obj, v_indexers, method=method, tolerance=tolerance
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/xarray/core/indexing.py", line 117, in remap_label_indexers
    idxr, new_idx = index.query(labels, method=method, tolerance=tolerance)
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/xarray/core/indexes.py", line 225, in query
    label_value, method=method, tolerance=tolerance
  File "/Users/osamuyuubu/anaconda3/envs/xesmf_env/lib/python3.7/site-packages/pandas/core/indexes/base.py", line 3363, in get_loc
    raise KeyError(key) from err
KeyError: 2087687.75

Can anyone help me? Thanks in advance.

I cant acces that ftp. Do you have another place to share the nc file? — Bert Coerver
– Bert Coerver, Commented Jan 20, 2022 at 15:53
Hi, Bert! @BertCoerver The data can be downloaded here, any nc file is ok, they have the same lon/lat, respectively. Thanks in advance! — Yongwu Xiu
– Yongwu Xiu, Commented Jan 21, 2022 at 2:02
try adding method="nearest" to .sel like SAT.sel(xc=x,yc=y, method="nearest") which will give you the nearest point to your coordinates, which is most likely the one you are looking for as indexing with exact coordinates is near impossible in that setting — Val
– Val, Commented Jan 21, 2022 at 8:15
Hi, Val, the method you provide doesn't come up with an error but the values of x and y are somehow strange, with x=-1098463 and y=-193688. It should be noted that the xc and yc of the dataarray range from -5000 to 5000, which are much smaller than x and y. So, thanks your advice, but I think there is something wrong with the method, which we should be careful. — Yongwu Xiu
– Yongwu Xiu, Commented Jan 21, 2022 at 8:30
well yes, you need to correctly transform your lat/lon coordinates into the array's coordinate system and make sure that your test point is even within the array's extent — Val
– Val, Commented Jan 21, 2022 at 13:32

Bert Coerver · Accepted Answer · 2022-01-21 15:32:31Z

1

Don't think you can use ds.analysis_sea_ice_thickness.sel(lat=80.0, lon=100.0, method = "nearest") directly, because lat and lon are not dimensions (see ds.dims), but coordinates (ds.coords). So you could do this for example:

import xarray as xr
import numpy as np

# Open dataset
ds = xr.open_dataset(r"/Path/to/your/folder/W_XX-ESA,SMOS_CS2,NH_25KM_EASE2_20211218_20211224_o_v204_01_l4sit.nc")

# Define point-of-interest.
lat = 80.0
lon = 100.0

# Find indices where lon and lat are closest to point-of-interest.
idxs = (np.abs(ds.lon - lon) + np.abs(ds.lat - lat)).argmin(dim = ["xc", "yc"])

# Retrieve value of variable at indices
value = ds.analysis_sea_ice_thickness.isel(idxs).values

# Check the actual lat and lon
lat_in_ds = ds.lat.isel(idxs).values
lon_in_ds = ds.lon.isel(idxs).values

# Print some results.
print(f"Thickness at ({lat_in_ds:.3f}, {lon_in_ds:.3f}) = {value[0]} {ds.analysis_sea_ice_thickness.units}.")

Thickness at (80.107, 99.782) = 0.805 m.

edited Jan 21, 2022 at 15:32

answered Jan 21, 2022 at 15:26

Bert Coerver

2821 silver badge10 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Yongwu Xiu Over a year ago

Thanks, this method works fine.

Bert Coerver Over a year ago

ok, feel free to accept it as the correct answer...

Yongwu Xiu Over a year ago

Sorry for a little bit

Collectives™ on Stack Overflow

xarray interp variable with two dimensional gird into one point

1 Answer 1

3 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related