1

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.

Yongwu Xiu
  • 125
  • 1
  • 9
  • I cant acces that ftp. Do you have another place to share the nc file? – Bert Coerver Jan 20 '22 at 15:53
  • Hi, Bert! @BertCoerver The data can be downloaded [here](ftp://ftp.awi.de/sea_ice/product/cryosat2_smos/v204/nh/LATEST/), any nc file is ok, they have the same lon/lat, respectively. Thanks in advance! – Yongwu Xiu Jan 21 '22 at 02: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 Jan 21 '22 at 08: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 Jan 21 '22 at 08: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 Jan 21 '22 at 13:32
  • Hi, Val, that is where the problem exists. I'm new to the projection, and the attempt follows the [answer](https://stackoverflow.com/questions/58758480/xarray-select-nearest-lat-lon-with-multi-dimension-coordinates/64431809#64431809) comes along with an error. And I'm tarpped here. – Yongwu Xiu Jan 21 '22 at 13:39

1 Answers1

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.

Bert Coerver
  • 252
  • 1
  • 9