[Bug]: Regridding error from differing time dimensions

[acordonez]

What happened?

In this case, I am using the Dataset.regridder.horizontal(tool='regrid2') method to regrid a reference dataset to the grid of a test dataset. The time dimensions are different lengths, for example the old grid data has 20 years and the target grid data has 50 years.
I do not have this issue when the time dimensions are the same length (the years can differ, but length must match), or when the target grid does not have a time dimension.

Running this command (target has 50 years, old has 20 years):
>> regridded = old.regridder.horizontal("ANN", target, tool='regrid2')

Throws:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[50], line 1
----> 1 regridded = old.regridder.horizontal("ANN", target, tool='regrid2')

File ~/miniconda3/envs/pcmdi_metrics_dev/lib/python3.9/site-packages/xcdat/regridder/accessor.py:316, in RegridderAccessor.horizontal(self, data_var, output_grid, tool, **options)
    311     raise ValueError(
    312         f"Tool {e!s} does not exist, valid choices {list(REGRID_TOOLS)}"
    313     )
    315 regridder = regrid_tool(self._ds, output_grid, **options)
--> 316 output_ds = regridder.horizontal(data_var, self._ds)
    318 return output_ds

File ~/miniconda3/envs/pcmdi_metrics_dev/lib/python3.9/site-packages/xcdat/regridder/regrid2.py:141, in Regrid2Regridder.horizontal(self, data_var, ds)
    137 ordered_axis_names = list(output_axis_sizes)
    139 output_data = self._regrid(input_data, output_axis_sizes, ordered_axis_names)
--> 141 output_ds = self._create_output_dataset(
    142     ds, data_var, output_data, axis_variable_name_map, ordered_axis_names
    143 )
    145 dst_mask = self._output_grid.get("mask", None)
    147 if dst_mask is not None:

File ~/miniconda3/envs/pcmdi_metrics_dev/lib/python3.9/site-packages/xcdat/regridder/regrid2.py:316, in Regrid2Regridder._create_output_dataset(self, input_ds, data_var, output_data, axis_variable_name_map, ordered_axis_names)
    313     else:
    314         coords[variable_name] = input_ds[variable_name].copy()
--> 316 output_da = xr.DataArray(
    317     output_data,
    318     dims=[axis_variable_name_map[x] for x in ordered_axis_names],
    319     coords=coords,
    320 )
    322 data_vars = {data_var: output_da}
    324 return xr.Dataset(data_vars, attrs=input_ds.attrs.copy())

File ~/miniconda3/envs/pcmdi_metrics_dev/lib/python3.9/site-packages/xarray/core/dataarray.py:417, in __init__(self, data, coords, dims, name, attrs, indexes, fastpath)
    415 if name is None:
    416     name = getattr(data, "name", None)
--> 417 if attrs is None and not isinstance(data, PANDAS_TYPES):
    418     attrs = getattr(data, "attrs", None)
    420 data = _check_data_shape(data, coords, dims)

File ~/miniconda3/envs/pcmdi_metrics_dev/lib/python3.9/site-packages/xarray/core/dataarray.py:163, in _infer_coords_and_dims(shape, coords, dims)
    158 for k, v in new_coords.items():
    159     if any(d not in dims for d in v.dims):
    160         raise ValueError(
    161             f"coordinate {k} has dimensions {v.dims}, but these "
    162             "are not a subset of the DataArray "
--> 163             f"dimensions {dims}"
    164         )
    166     for d, s in zip(v.dims, v.shape):
    167         if s != sizes[d]:

ValueError: conflicting sizes for dimension 'time': length 50 on the data but length 20 on coordinate 'time'

This does work when I slice "target" to be 20 years:
regridded = old.regridder.horizontal("ANN", target.sel(time=slice("1950-01-01","1969-12-31")), tool='regrid2')

This also works when I strip out the time coordinate:
regridded = old.regridder.horizontal("ANN", target.isel({"time":0}).drop("time"), tool='regrid2')

What did you expect to happen? Are there are possible answers you came across?

I would expect the time dimension to not affect horizontal regridding, or to get a clearer error message if the target grid should not have any time coordinates.

Minimal Complete Verifiable Example (MVCE)

import xcdat as xc

# I've been trying to create dummy arrays in-place that 
# reproduce this issue but haven't been successful.

# These files are actually on the same grid but serve as
# an example.
path1="/home/ordonez4/git/pcmdi_metrics/pcmdi_metrics/extremes/test_nc_files/netcdf/MRI-ESM2-0_r1i1p1f1_land_TXn.nc"
path2="/home/ordonez4/git/pcmdi_metrics/pcmdi_metrics/extremes/test_nc_files/netcdf/Reference_fake_tasmax_land_TXn.nc"

# 50 year Dataset
target = xc.open_dataset(path1)
# 20 year Dataset
old = xc.open_dataset(path2)

regridded = old.regridder.horizontal("ANN", target, tool='regrid2')

Relevant log output

No response

Anything else we need to know?

Preview of one of my data sets (they both have the same format, other than different dimensions):

Environment

xcdat = 0.5.0

INSTALLED VERSIONS

commit: None
python: 3.9.16 | packaged by conda-forge | (main, Feb 1 2023, 21:39:03)
[GCC 11.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-1160.71.1.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.12.2
libnetcdf: 4.8.1

xarray: 2023.4.2
pandas: 1.5.3
numpy: 1.21.3
scipy: 1.10.1
netCDF4: 1.6.2
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.6.2
nc_time_axis: 1.4.1
PseudoNetCDF: None
iris: None
bottleneck: 1.3.6
dask: 2023.2.0
distributed: 2023.2.0
matplotlib: 3.4.3
cartopy: 0.21.1
seaborn: 0.11.1
numbagg: None
fsspec: 2023.4.0
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 67.7.2
pip: 21.2.4
conda: None
pytest: 6.2.5
mypy: None
IPython: 8.13.2
sphinx: None

[pochedls]

I usually pass in a target_grid object with no time information. Can you do something like this:

ds_target = xc.open_dataset(path1)
target = xc.create_grid(ds_target.lat, ds_target.lon)
...
regridded = old.regridder.horizontal("ANN", target, tool='regrid2')

Also note #447 is on regrid2 (though I think the issue is different). Also...there is an active discussion about changing the API for xc.create_grids() in #388 / #496 – something to be aware of.

[lee1043]

Just curious, can you try below in your code?

regridded = old.regridder.horizontal("ANN", target.isel(time=0), tool='regrid2')

OR

regridded = old.regridder.horizontal("ANN", target.isel(time=0).squeeze(), tool='regrid2')

Thanks to @pochedls for the above second option, which I think it makes more sense because squeeze will drop time coordinate thus horizontal regridder will not influenced by it.

[acordonez]

@lee1043 It might have gotten lost in all my text, but pulling off the time coordinate like this does work:
regridded = old.regridder.horizontal("ANN", target.isel({"time":0}).drop("time"), tool='regrid2')

But @pochedls suggestion to make a new grid might be more robust, and it also worked in my case. I think #447 is different from my issue, but thanks for the heads-up on what else is happening in the regridding area.

[acordonez]

I mainly wanted to flag that I was seeing inconsistent results depending on the time axis, and if the target grid really shouldn't have any time dimension it would be helpful to see that stated more explicitly.

[lee1043]

@acordonez thank you for clarifying, I missed that part. Thank you for flagging this out!

[pochedls]

I think the main issue/takeaway here is to improve the documentation, which I think may be out-of-date. The documentation lists a data_var parameter and a ds parameter, though I think ds should be something like target_grid with a note saying that it is a dataset that contains the regridding coordinates (and if it includes time, the time dimension should be the same length as the source dataset). Is this right?

[acordonez]

@pochedls That summary sounds right to me, thanks for putting it together.

[lee1043]

@pochedls That could be one way, and I can also think of benefit of giving ds as an input for regridder, when a user has loaded just 2 datasets and regrid one to another, without writing additional row of a code for target_grid. I think another alternative way could be revising the spatial regridder to drop the time coordinate from given ds if exist. What do you think?

[pochedls]

I wasn't intending to comment on the name ds versus target_grid. I think it would be fine to pass in ds (and call it ds), but the documentation for horizontal is a little confusing. It says

ds (xr.Dataset) – The dataset containing data_var.

which I think isn't right?

I'd be fine with either:

• updating the documentation (to say that ds should not have a time dimension or the time dimension should be the same size as the input dataset)
• automatically dropping the time dimension of the target grid

In either case I think the documentation should be updated to say that

ds (xr.Dataset) – The dataset containing the target grid information

(if I am right that the current documentation is incomplete)

[jasonb5]

@pochedls @acordonez Circling around to this issue, I'm updating the documentation in another PR.

There's another issue here, the regrid2 code assumes output_grid has no time dimension, I'll add a check and raise an appropriate error.

Both @pochedls and @lee1043 solutions work, the intended way to handle this is:

regridded = old.regridder.horizontal("ANN", target.regridder.grid, tool='regrid2')

The grid property is intended to extract the grid information from any dataset.