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