`find_optimal_celestial_wcs` fails with a stack of solar WCSes, when specifying a frame which is not the first element of the sequence

[Cadair]

This is a slightly weird bug, because it's not really reproject's fault, I'm also tracking this here: sunpy/sunpy#8494

This is the code which triggers the bug:

import sunpy.coordinates
from sunpy.data.sample import *
from astropy.io import fits
from astropy.wcs import WCS
from astropy.wcs.utils import wcs_to_celestial_frame
from reproject.mosaicking import find_optimal_celestial_wcs

file_seq = [AIA_193_CUTOUT01_IMAGE, AIA_193_CUTOUT02_IMAGE, AIA_193_CUTOUT03_IMAGE, AIA_193_CUTOUT04_IMAGE, AIA_193_CUTOUT05_IMAGE]

hdus = [fits.open(f)[1] for f in file_seq]
wcses = [WCS(h) for h in hdus]
frames = [wcs_to_celestial_frame(w) for w in wcses]

# This works
key_frame = frames[0]
# This does not
key_frame = frames[3]

with sunpy.coordinates.screens.SphericalScreen(key_frame.observer):
    wcs, shape = find_optimal_celestial_wcs(hdus, frame=key_frame, negative_lon_cdelt="auto")

The error raised is:

Show Error

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity.py:964, in Quantity.to_value(self, unit, equivalencies)
    963 try:
--> 964     scale = self.unit._to(unit)
    965 except Exception:
    966     # Short-cut failed; try default (maybe equivalencies help).

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/core.py:624, in UnitBase._to(self, other)
    622         return self_decomposed.scale / other_decomposed.scale
--> 624 raise UnitConversionError(f"'{self!r}' is not a scaled version of '{other!r}'")

UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("m")'

During handling of the above exception, another exception occurred:

UnitConversionError                       Traceback (most recent call last)
File <string>:20

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/reproject/mosaicking/wcs_helpers.py:206, in find_optimal_celestial_wcs(input_data, hdu_in, frame, auto_rotate, projection, resolution, reference, negative_lon_cdelt)
    202     resolutions.append(pixel_scale(wcs, shape))
    204 # We now stack the coordinates - however the frame classes can't do this
    205 # so we have to use the high-level SkyCoord class.
--> 206 corners = SkyCoord(corners)
    207 references = SkyCoord(references)
    209 # If no reference coordinate has been passed in for the final header, we
    210 # determine the reference coordinate as the mean of all the reference
    211 # positions. This choice is as good as any and if the user really cares,
    212 # they can set  it manually.

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/coordinates/sky_coordinate.py:229, in SkyCoord.__init__(self, copy, *args, **kwargs)
    225 # Parse the args and kwargs to assemble a sanitized and validated
    226 # kwargs dict for initializing attributes for this object and for
    227 # creating the internal self._sky_coord_frame object
    228 args = list(args)  # Make it mutable
--> 229 skycoord_kwargs, components, info = _parse_coordinate_data(
    230     frame_cls(**frame_kwargs), args, kwargs
    231 )
    233 # In the above two parsing functions, these kwargs were identified
    234 # as valid frame attributes for *some* frame, but not the frame that
    235 # this SkyCoord will have. We keep these attributes as special
    236 # skycoord frame attributes:
    237 for attr in skycoord_kwargs:
    238     # Setting it will also validate it.

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/coordinates/sky_coordinate_parsers.py:278, in _parse_coordinate_data(frame, args, kwargs)
    271         _components[frame_attr_name] = attr_class(arg, unit=unit)
    273 elif len(args) == 1:
    274     # One arg which must be a coordinate.  In this case coord_kwargs
    275     # will contain keys like 'ra', 'dec', 'distance' along with any
    276     # frame attributes like equinox or obstime which were explicitly
    277     # specified in the coordinate object (i.e. non-default).
--> 278     _skycoord_kwargs, _components = _parse_coordinate_arg(args[0], frame, units)
    280     # Copy other 'info' attr only if it has actually been defined.
    281     if "info" in getattr(args[0], "__dict__", ()):

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/coordinates/sky_coordinate_parsers.py:473, in _parse_coordinate_arg(coords, frame, units)
    469         skycoord_kwargs[fattrnm] = getattr(scs[0], fattrnm)
    471     # Now combine the values, to be used below
    472     values = [
--> 473         np.concatenate([np.atleast_1d(getattr(sc, data_attr)) for sc in scs])
    474         for data_attr, repr_attr in zip(frame_attr_names, repr_attr_names)
    475         if not_unit_sphere or repr_attr != "distance"
    476     ]
    477 else:
    478     is_radec = (
    479         "ra" in frame.representation_component_names
    480         and "dec" in frame.representation_component_names
    481     )

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity.py:1884, in Quantity.__array_function__(self, function, types, args, kwargs)
   1882 function_helper = FUNCTION_HELPERS[function]
   1883 try:
-> 1884     args, kwargs, unit, out = function_helper(*args, **kwargs)
   1885 except NotImplementedError:
   1886     return self._not_implemented_or_raise(function, types)

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity_helper/function_helpers.py:498, in concatenate(arrays, axis, out, **kwargs)
    494 @function_helper
    495 def concatenate(arrays, /, axis=0, out=None, **kwargs):
    496     # TODO: make this smarter by creating an appropriately shaped
    497     # empty output array and just filling it.
--> 498     arrays, kwargs, unit, out = _iterable_helper(
    499         *arrays, out=out, axis=axis, **kwargs
    500     )
    501     return (arrays,), kwargs, unit, out

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity_helper/function_helpers.py:478, in _iterable_helper(out, *args, **kwargs)
    475 if out is not None:
    476     kwargs["out"] = _quantity_out_as_array(out)  # raises if not Quantity.
--> 478 arrays, unit = _quantities2arrays(*args)
    479 return arrays, kwargs, unit, out

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity_helper/function_helpers.py:466, in _quantities2arrays(unit_from_first, *args)
    462 # We use the private _to_own_unit method here instead of just
    463 # converting everything to quantity and then do .to_value(qs0.unit)
    464 # as we want to allow arbitrary unit for 0, inf, and nan.
    465 try:
--> 466     arrays = tuple((q._to_own_unit(arg)) for arg in args)
    467 except TypeError:
    468     raise NotImplementedError

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity_helper/function_helpers.py:466, in <genexpr>(.0)
    462 # We use the private _to_own_unit method here instead of just
    463 # converting everything to quantity and then do .to_value(qs0.unit)
    464 # as we want to allow arbitrary unit for 0, inf, and nan.
    465 try:
--> 466     arrays = tuple((q._to_own_unit(arg)) for arg in args)
    467 except TypeError:
    468     raise NotImplementedError

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity.py:1685, in Quantity._to_own_unit(self, value, check_precision, unit)
   1683     unit = self.unit
   1684 try:
-> 1685     _value = value.to_value(unit)
   1686 except AttributeError:
   1687     # We're not a Quantity.
   1688     # First remove two special cases (with a fast test):
   1689     # 1) Maybe masked printing? MaskedArray with quantities does not
   1690     # work very well, but no reason to break even repr and str.
   1691     # 2) np.ma.masked? useful if we're a MaskedQuantity.
   1692     if value is np.ma.masked or (
   1693         value is np.ma.masked_print_option and self.dtype.kind == "O"
   1694     ):

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity.py:967, in Quantity.to_value(self, unit, equivalencies)
    964     scale = self.unit._to(unit)
    965 except Exception:
    966     # Short-cut failed; try default (maybe equivalencies help).
--> 967     value = self._to_value(unit, equivalencies)
    968 else:
    969     value = self.view(np.ndarray)

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/quantity.py:873, in Quantity._to_value(self, unit, equivalencies)
    870     equivalencies = self._equivalencies
    871 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    872     # Standard path, let unit to do work.
--> 873     return self.unit.to(
    874         unit, self.view(np.ndarray), equivalencies=equivalencies
    875     )
    877 else:
    878     # The .to() method of a simple unit cannot convert a structured
    879     # dtype, so we work around it, by recursing.
    880     # TODO: deprecate this?
    881     # Convert simple to Structured on initialization?
    882     result = np.empty_like(self.view(np.ndarray))

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/core.py:660, in UnitBase.to(self, other, value, equivalencies)
    658     return UNITY
    659 else:
--> 660     return self.get_converter(Unit(other), equivalencies)(value)

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/core.py:589, in UnitBase.get_converter(self, other, equivalencies)
    586             else:
    587                 return lambda v: b(converter(v))
--> 589 raise exc

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/core.py:572, in UnitBase.get_converter(self, other, equivalencies)
    570 # if that doesn't work, maybe we can do it with equivalencies?
    571 try:
--> 572     return self._apply_equivalencies(
    573         self, other, self._normalize_equivalencies(equivalencies)
    574     )
    575 except UnitsError as exc:
    576     # Last hope: maybe other knows how to do it?
    577     # We assume the equivalencies have the unit itself as first item.
    578     # TODO: maybe better for other to have a `_back_converter` method?
    579     if hasattr(other, "equivalencies"):

File ~/.local/share/mamba/envs/sunpy-dev/lib/python3.14/site-packages/astropy/units/core.py:523, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
    520 unit_str = get_err_str(unit)
    521 other_str = get_err_str(other)
--> 523 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: '' (dimensionless) and 'm' (length) are not convertible

but the failure is actually described better in the sunpy issue linked above.

If we decide to work around this in reproject, my thought is that the best thing to do would be to cast all the frames to unitspherical representation before stacking, as we are only dealing with celestial WCSes in this function anyway?