xarray interp 变量与二维网格成一个点

Question

我有一个二维 lon/lat 的数据集，我想计算某个点的值。

ncfile可以从ftp://ftp.awi.de/sea_ice/product/cryosat2_smos/v204/nh/LATEST/下载，变量可以通过以下方式获取

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

SAT_SIT 表示为

<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

现在，我想检查 SAT_SIT 在某个特定点的值，例如 lon=100 lat=80。有没有可能的方法来处理这个？

注意：xarray 只支持一维坐标的插值，"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."

在link 中有一些可能的方法来解决这个问题。但是，第一种方法有点复杂，因为我需要将变量插入到很多点。第二种方法出现错误。

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)

错误是

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

谁能帮帮我？提前致谢。

Answer 1

不要以为可以直接使用ds.analysis_sea_ice_thickness.sel(lat=80.0, lon=100.0, method = "nearest")，因为lat 和lon 不是维度（参见ds.dims），而是坐标（ds.coords）。例如，您可以这样做：

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.