Search code examples
pythoninterpolationpython-xarrayspatial-interpolation

xarray interp variable with two dimensional gird into one point

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.

Solution

  • 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.