use xESMF for interpolation on a section (2D to 1D interpolation)

[szwang1990]

Hi, I've been using xESMF for planar interpolation (2D to 2D). I just wonder is it possible for xESMF to interpolate a planar field onto a section (2D to 1D)? I'm asking because when I want to do the 2D to 1D interpolation, I got errors.

ds_out = xr.Dataset(
 {
 'lat': (['lat'], np.arange(0, 90, 1), {'units': 'degrees_north'}),
 'lon': (['lon'], np.arange(0, 180, 2), {'units': 'degrees_east'})
 },
)

s_an = xr.open_dataset(dname+fname, decode_times=False).squeeze().s_an

regridder = xe.Regridder(s_an, ds_out, 'conservative')

Then I got the following errors.

/home/swang/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/cf_xarray/accessor.py:1772: UserWarning: Variables {'lon_bnds'} not found in object but are referred to in the CF attributes.
warnings.warn(
--------------------------------------------------------------------------- KeyError Traceback (most recent call last)
File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xesmf/frontend.py:69, in _get_lon_lat_bounds(ds) 68 try: ---> 69 lon_bnds = ds.cf.get_bounds('longitude') 70 lat_bnds = ds.cf.get_bounds('latitude')

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/cf_xarray/accessor.py:2287, in CFDatasetAccessor.get_bounds(self, key) 2286 if not results: -> 2287 raise KeyError(f"No results found for {key!r}.") 2289 return self._obj[results[0] if len(results) == 1 else results] KeyError: "No results found for 'longitude'."

During handling of the above exception, another exception occurred: ValueError Traceback (most recent call last)
Cell In[126], line 1 ----> 1 regridder = xe.Regridder(s_an, ds_out, 'conservative') File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xesmf/frontend.py:785, in Regridder.init(self, ds_in, ds_out, method, locstream_in, locstream_out, periodic, **kwargs) 783 grid_out, shape_out, _ = ds_to_ESMFlocstream(ds_out) 784 else: --> 785 grid_out, shape_out, _ = ds_to_ESMFgrid(ds_out, need_bounds=need_bounds) 787 # Create the BaseRegridder 788 super().init(grid_in, grid_out, method, input_dims=input_dims, **kwargs)

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xesmf/frontend.py:136, in ds_to_ESMFgrid(ds, need_bounds, periodic, append) 133 grid = Grid.from_xarray(lon.T, lat.T, periodic=periodic, mask=None) 135 if need_bounds: --> 136 lon_b, lat_b = _get_lon_lat_bounds(ds) 137 lon_b, lat_b = as_2d_mesh(np.asarray(lon_b), np.asarray(lat_b)) 138 add_corner(grid, lon_b.T, lat_b.T)

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xesmf/frontend.py:77, in _get_lon_lat_bounds(ds) 75 lon_name = ds.cf['longitude'].name 76 lat_name = ds.cf['latitude'].name ---> 77 ds = ds.cf.add_bounds([lon_name, lat_name]) 78 lon_bnds = ds.cf.get_bounds('longitude') 79 lat_bnds = ds.cf.get_bounds('latitude')

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/cf_xarray/accessor.py:2387, in CFDatasetAccessor.add_bounds(self, keys, dim, output_dim) 2385 if bname in obj.variables: 2386 raise ValueError(f"Bounds variable name {bname!r} will conflict!") -> 2387 out = _guess_bounds( 2388 obj[var].reset_coords(drop=True), dim=dim, out_dim=output_dim 2389 ) 2390 if output_dim in obj.dims and (new := out[output_dim].size) != ( 2391 old := obj[output_dim].size 2392 ): 2393 raise ValueError( 2394 f"The {output_dim} dimension already exists but has a different length than the new one " 2395 f"({old} vs {new}). Please provide another bound dimension name with output_dim." 2396 )

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/cf_xarray/accessor.py:574, in _guess_bounds(da, dim, out_dim) 571 else: 572 dim = dim[0] --> 574 return _guess_bounds_1d(da, dim).rename(bounds=out_dim)

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/cf_xarray/helpers.py:26, in _guess_bounds_1d(da, dim) 23 ADDED_INDEX = True 25 diff = da.diff(dim) ---> 26 lower = da - diff / 2 27 upper = da + diff / 2 28 bounds = xr.concat([lower, upper], dim="bounds")

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/_typed_ops.py:225, in DataArrayOpsMixin.sub(self, other) 224 def sub(self, other): --> 225 return self._binary_op(other, operator.sub) File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/dataarray.py:4615, in DataArray._binary_op(self, other, f, reflexive) 4613 if isinstance(other, DataArray): 4614 align_type = OPTIONS["arithmetic_join"] -> 4615 self, other = align(self, other, join=align_type, copy=False) # type: ignore 4616 other_variable = getattr(other, "variable", other) 4617 other_coords = getattr(other, "coords", None)

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:787, in align(join, copy, indexes, exclude, fill_value, *objects) 591 """ 592 Given any number of Dataset and/or DataArray objects, returns new 593 objects with aligned indexes and dimension sizes. (...) 777 778 """ 779 aligner = Aligner( 780 objects, 781 join=join, (...) 785 fill_value=fill_value, 786 ) --> 787 aligner.align() 788 return aligner.results

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:580, in Aligner.align(self) 578 self.results = self.objects 579 else: --> 580 self.reindex_all() File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:555, in Aligner.reindex_all(self) 554 def reindex_all(self) -> None: --> 555 self.results = tuple( 556 self._reindex_one(obj, matching_indexes) 557 for obj, matching_indexes in zip( 558 self.objects, self.objects_matching_indexes 559 ) 560 ) File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:556, in (.0) 554 def reindex_all(self) -> None: 555 self.results = tuple( --> 556 self._reindex_one(obj, matching_indexes) 557 for obj, matching_indexes in zip( 558 self.objects, self.objects_matching_indexes 559 ) 560 )

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:540, in Aligner._reindex_one(self, obj, matching_indexes) 534 def _reindex_one( 535 self, 536 obj: DataAlignable, 537 matching_indexes: dict[MatchingIndexKey, Index], 538 ) -> DataAlignable: 539 new_indexes, new_variables = self._get_indexes_and_vars(obj, matching_indexes) --> 540 dim_pos_indexers = self._get_dim_pos_indexers(matching_indexes) 542 new_obj = obj._reindex_callback( 543 self, 544 dim_pos_indexers, (...) 549 self.exclude_vars, 550 ) 551 new_obj.encoding = obj.encoding

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/alignment.py:506, in Aligner._get_dim_pos_indexers(self, matching_indexes) 504 if obj_idx is not None: 505 if self.reindex[key]: --> 506 indexers = obj_idx.reindex_like(aligned_idx, **self.reindex_kwargs) 507 dim_pos_indexers.update(indexers) 509 return dim_pos_indexers

File ~/Software/anaconda3/envs/xesmf_env/lib/python3.10/site-packages/xarray/core/indexes.py:527, in PandasIndex.reindex_like(self, other, method, tolerance) 523 def reindex_like( 524 self, other: PandasIndex, method=None, tolerance=None 525 ) -> dict[Hashable, Any]: 526 if not self.index.is_unique: --> 527 raise ValueError( 528 f"cannot reindex or align along dimension {self.dim!r} because the " 529 "(pandas) index has duplicate values" 530 ) 532 return {self.dim: get_indexer_nd(self.index, other.index, method, tolerance)} ValueError: cannot reindex or align along dimension 'lon' because the (pandas) index has duplicate values

[sol1105]

There seems to be an issue when cf_xarray tries to calculate the missing bounds for lat and lon. Bounds are required when generating conservative weights.

As a workaround you could have xesmf create the output grid for you, incl. the bounds:

import xesmf as xe
# With cells centered at 0.5, 1.5 ... (lat) and 1, 3, ... (lon) 
ds_out=xe.util.cf_grid_2d(0, 180, 2, 0, 90, 1)
# With cells centered at 0, 1, 2 (lat) and 0, 2, 4, ... (lon) - so same as you defined it using np.arange
ds_out=xe.util.cf_grid_2d(-1, 179, 2, -1, 89, 1)

Then you can again try generating the weights.

[raphaeldussin]

@szwang1990 have you tried using the locstream output option?
https://xesmf.readthedocs.io/en/stable/notebooks/Using_LocStream.html