transform.rescale no longer handles xarray object directly if preserve_range enabled

[yichechang]

Description:

Starting in 0.21.0, passing an xarray object (xarray.DataArray) into transform.rescale no longer works if preserve_range is set to True. Like:

rescale(..., scale=..., preserve_range=True)

Above would work in both 0.19.3 and 0.20.0.

More detail can be found below, but before that, I just wanted to check if scikit-image would like to support passing say xarray.DataArray into most scikit-image functions? Because if that's not within the scope, it's perfectly fine that I can always pass in the underlying numpy array via DataArray's .data attribute. (In that case, I'll simply close this issue...)

---

Traceback

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[51], line 1
----> 1 rescale(image, scale=0.5, preserve_range=True)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/skimage/_shared/utils.py:316, in channel_as_last_axis.__call__.<locals>.fixed_func(*args, **kwargs)
    313 channel_axis = kwargs.get('channel_axis', None)
    315 if channel_axis is None:
--> 316     return func(*args, **kwargs)
    318 # TODO: convert scalars to a tuple in anticipation of eventually
    319 #       supporting a tuple of channel axes. Right now, only an
    320 #       integer or a single-element tuple is supported, though.
    321 if np.isscalar(channel_axis):

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/skimage/transform/_warps.py:289, in rescale(image, scale, order, mode, cval, clip, preserve_range, anti_aliasing, anti_aliasi
ng_sigma, channel_axis)
    286 if multichannel:  # don't scale channel dimension
    287     output_shape[-1] = orig_shape[-1]
--> 289 return resize(image, output_shape, order=order, mode=mode, cval=cval,
    290               clip=clip, preserve_range=preserve_range,
    291               anti_aliasing=anti_aliasing,
    292               anti_aliasing_sigma=anti_aliasing_sigma)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/skimage/transform/_warps.py:188, in resize(image, output_shape, order, mode, cval, clip, preserve_range, anti_aliasing, anti_
aliasing_sigma)
    184 zoom_factors = [1 / f for f in factors]
    185 out = ndi.zoom(filtered, zoom_factors, order=order, mode=ndi_mode,
    186                cval=cval, grid_mode=True)
--> 188 _clip_warp_output(image, out, mode, cval, clip)
    190 return out

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/skimage/transform/_warps.py:692, in _clip_warp_output(input_image, output_image, mode, cval, clip)
    689     min_val = min(min_val, cval)
    690     max_val = max(max_val, cval)
--> 692 np.clip(output_image, min_val, max_val, out=output_image)

File <__array_function__ internals>:200, in clip(*args, **kwargs)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/numpy/core/fromnumeric.py:2180, in clip(a, a_min, a_max, out, **kwargs)
   2111 @array_function_dispatch(_clip_dispatcher)
   2112 def clip(a, a_min, a_max, out=None, **kwargs):
   2113     """
   2114     Clip (limit) the values in an array.
   2115
   (...)
   2178
   2179     """
-> 2180     return _wrapfunc(a, 'clip', a_min, a_max, out=out, **kwargs)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/numpy/core/fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds)
     54     return _wrapit(obj, method, *args, **kwds)
     56 try:
---> 57     return bound(*args, **kwds)
     58 except TypeError:
     59     # A TypeError occurs if the object does have such a method in its
     60     # class, but its signature is not identical to that of NumPy's. This
   (...)
     64     # Call _wrapit from within the except clause to ensure a potential
     65     # exception has a traceback chain.
     66     return _wrapit(obj, method, *args, **kwds)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/numpy/core/_methods.py:161, in _clip(a, min, max, out, casting, **kwargs)
    158     return _clip_dep_invoke_with_casting(
    159         um.maximum, a, min, out=out, casting=casting, **kwargs)
    160 else:
--> 161     return _clip_dep_invoke_with_casting(
    162         um.clip, a, min, max, out=out, casting=casting, **kwargs)

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/numpy/core/_methods.py:115, in _clip_dep_invoke_with_casting(ufunc, out, casting, *args, **kwargs)
    113 # try to deal with broken casting rules
    114 try:
--> 115     return ufunc(*args, out=out, **kwargs)
    116 except _exceptions._UFuncOutputCastingError as e:
    117     # Numpy 1.17.0, 2019-02-24
    118     warnings.warn(
    119         "Converting the output of clip from {!r} to {!r} is deprecated. "
    120         "Pass `casting=\"unsafe\"` explicitly to silence this warning, or "
   (...)
    123         stacklevel=2
    124     )

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/xarray/core/arithmetic.py:86, in SupportsArithmetic.__array_ufunc__(self, ufunc, method, *inputs, **kwargs)
     77     raise NotImplementedError(
     79         "for ufuncs. As an alternative, consider explicitly "
     80         "converting xarray objects to NumPy arrays (e.g., with "
     81         "`.values`)."
     82     )
     84 join = dataset_join = OPTIONS["arithmetic_join"]
---> 86 return apply_ufunc(
     87     ufunc,
     88     *inputs,
     89     input_core_dims=((),) * ufunc.nin,
     90     output_core_dims=((),) * ufunc.nout,
     91     join=join,
     92     dataset_join=dataset_join,
     93     dataset_fill_value=np.nan,
     94     kwargs=kwargs,
     95     dask="allowed",
     96     keep_attrs=_get_keep_attrs(default=True),
     97 )

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/xarray/core/computation.py:1208, in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, datas
et_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args)
   1206 # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
   1207 elif any(isinstance(a, DataArray) for a in args):
-> 1208     return apply_dataarray_vfunc(
   1209         variables_vfunc,
   1210         *args,
   1211         signature=signature,
   1212         join=join,
   1213         exclude_dims=exclude_dims,
   1214         keep_attrs=keep_attrs,
   1215     )
   1216 # feed Variables directly through apply_variable_ufunc
   1217 elif any(isinstance(a, Variable) for a in args):

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/xarray/core/computation.py:315, in apply_dataarray_vfunc(func, signature, join, exclude_dims, keep_attrs, *args)
    310 result_coords, result_indexes = build_output_coords_and_indexes(
    311     args, signature, exclude_dims, combine_attrs=keep_attrs
    312 )
    314 data_vars = [getattr(a, "variable", a) for a in args]
--> 315 result_var = func(*data_vars)
    317 out: tuple[DataArray, ...] | DataArray
    318 if signature.num_outputs > 1:

File ~/mambaforge/envs/test_skimage21/lib/python3.8/site-packages/xarray/core/computation.py:796, in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_a
ttrs, dask_gufunc_kwargs, *args)
    794 data = as_compatible_data(data)
    795 if data.ndim != len(dims):
--> 796     raise ValueError(
    797         "applied function returned data with unexpected "
    798         f"number of dimensions. Received {data.ndim} dimension(s) but "
    799         f"expected {len(dims)} dimensions with names: {dims!r}"
    800     )
    802 var = Variable(dims, data, fastpath=True)
    803 for dim, new_size in var.sizes.items():

ValueError: applied function returned data with unexpected number of dimensions. Received 2 dimension(s) but expected 0 dimensions with names: ()

My understanding is that in skimage/transform/_warp.py, in function _clip_warp_output (called by resize/rescale) there's

np.clip(output_image, min_val, max_val, out=output_image)

that is now broken in 0.21.0.

(Specifically, I think this behavior might be introduced by PR #6852 while fixing other things?)

Way to reproduce:

import numpy as np
import xarray as xr
from skimage.transform import rescale

image = xr.DataArray(np.ones((8, 8)))

# this works with DataArray passed in directly
rescale(image, scale=0.5, preserve_range=False)

# but it doesn't work if turning preserve_range on
rescale(image, scale=0.5, preserve_range=True)

# again it works if we pass the underlying ndarray in
rescale(image.data, scale=0.5, preserve_range=True)

Version information:

3.8.17 | packaged by conda-forge | (default, Jun 16 2023, 07:11:32)
[Clang 14.0.6 ]
macOS-13.5.1-arm64-arm-64bit
scikit-image version: 0.21.0
numpy version: 1.24.4

[soupault]

Hi @yichechang, thank you for reporting and investigating the issue! I would say that native support for xr.DataArrays is something we may consider in the future, but not before they are widely supported by numpy and scipy (i.e. without issues similar to this one). We would also need to properly think through the assumptions we make about the users' xarrays (e.g., how we treat DataArray.coords and .dims), such that the library functions work in a consistent and expected way. That sounds like a big amount of work, and I would prefer to start it with specification (https://scientific-python.org/specs/).

Regarding the regression you reported, I have opened a PR with the fix - #7159. Please, feel free to jump in and make a review.

[yichechang]

Hi @soupault,

Thank you for your time and the PR! I'm currently traveling but plan to test it with my own data this weekend.

I also want to say I agree with your thoughts on supporting xarray and the first step being coming up with a spec. I was pleasantly surprised when things worked in skimage after accidentally passing an xr.DataArray instead of its .data attribute. I began using DataArray objects in some skimage functions without much verification, until this issue arose. So, I should make it clear that it’s not a bug and it didn't break anything, but rather an inconsistency for an unsupported array type.

[soupault]

So, I should make it clear that it’s not a bug and it didn't break anything, but rather an inconsistency for an unsupported array type

I am of the same opinion. However, since the fix is trivial, it is fine to have it.

[yichechang]

Apologies for my delayed reply. I think the reasoning and implementation are correct, and based on some tests you posted in PR #7159 these look promising!

Thank you again for fixing this, even if it doesn't get merged at the end, I still want to express my appreciation!

(I attempted to test with my real data locally but there was some build issues that I didn't know how to resolve when trying to install this branch locally.)