Pixel-to-pixel conversion failing because of doppler convention

I'm trying to reproject one cube to another. Both are in velocity coordinates and occupy very tiny regions on the sky.

I get this traceback:

```python
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/scratch/local/47745774/ipykernel_15582/2515218747.py in
2 # "ValueError: doppler_convention not set, cannot convert to/from velocities"
3 water_cube_spectral = water_cube.spectral_interpolate(stackAcube.spectral_axis)
----> 4 water_cube_regridA = water_cube_spectral.reproject(stackAcube.header)

/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/utils.py in wrapper(self, *args, **kwargs)
47 PossiblySlowWarning
48 )
---> 49 return function(self, *args, **kwargs)
50 return wrapper
51

/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/spectral_cube.py in reproject(self, header, order, use_memmap, filled, **kwargs)
2695 outarray = None
2696
-> 2697 newcube, newcube_valid = reproject_interp((data,
2698 self.header),
2699 newwcs,

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
544 warnings.warn(msg, warning_type, stacklevel=2)
545
--> 546 return function(*args, **kwargs)
547
548 return wrapper

/blue/adamginsburg/adamginsburg/repos/reproject/reproject/interpolation/high_level.py in reproject_interp(input_data, output_projection, shape_out, hdu_in, order, independent_celestial_slices, output_array, return_footprint, output_footprint, block_size, parallel, roundtrip_coords)
118 output_footprint=output_footprint)
119 else:
--> 120 return _reproject_full(array_in, wcs_in, wcs_out, shape_out=shape_out,
121 order=order, array_out=output_array,
122 return_footprint=return_footprint, roundtrip_coords=roundtrip_coords)

/blue/adamginsburg/adamginsburg/repos/reproject/reproject/interpolation/core.py in _reproject_full(array, wcs_in, wcs_out, shape_out, order, array_out, return_footprint, roundtrip_coords)
83 # For each pixel in the ouput array, get the pixel value in the input WCS
84 if roundtrip_coords:
---> 85 pixel_in = efficient_pixel_to_pixel_with_roundtrip(
86 wcs_out, wcs_in, *pixel_out[::-1])[::-1]
87 else:

/blue/adamginsburg/adamginsburg/repos/reproject/reproject/wcs_utils.py in efficient_pixel_to_pixel_with_roundtrip(wcs1, wcs2, *inputs)
230 def efficient_pixel_to_pixel_with_roundtrip(wcs1, wcs2, *inputs):
231
--> 232 outputs = efficient_pixel_to_pixel(wcs1, wcs2, *inputs)
233
234 # Now convert back to check that coordinates round-trip, if not then set to NaN

/blue/adamginsburg/adamginsburg/repos/reproject/reproject/wcs_utils.py in efficient_pixel_to_pixel(wcs1, wcs2, *inputs)
206 if not isinstance(world_outputs, (tuple, list)):
207 world_outputs = (world_outputs,)
--> 208 pixel_outputs = wcs2.world_to_pixel(*world_outputs)
209
210 for ipix in range(wcs2.pixel_n_dim):

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/wcs/wcsapi/high_level_api.py in world_to_pixel(self, *world_objects)
307 def world_to_pixel(self, *world_objects):
308
--> 309 world_values = high_level_objects_to_values(*world_objects, low_level_wcs=self.low_level_wcs)
310
311 # Finally we convert to pixel coordinates

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/wcs/wcsapi/high_level_api.py in high_level_objects_to_values(low_level_wcs, *world_objects)
231 for key, _, attr in components:
232 if callable(attr):
--> 233 world.append(attr(objects[key]))
234 else:
235 world.append(rec_getattr(objects[key], attr))

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/wcs/wcsapi/fitswcs.py in value_from_spectralcoord(spectralcoord)
600 return spectralcoord.to_value(**kwargs)
601 else:
--> 602 return spectralcoord.with_observer_stationary_relative_to(observer).to_value(**kwargs)
603
604 classes['spectral'] = (u.Quantity, (), {}, spectralcoord_from_value)

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/units/decorators.py in wrapper(*func_args, **func_kwargs)
300 # Call the original function with any equivalencies in force.
301 with add_enabled_equivalencies(self.equivalencies):
--> 302 return_ = wrapped_function(*func_args, **func_kwargs)
303
304 # Return

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/coordinates/spectral_coordinate.py in with_observer_stationary_relative_to(self, frame, velocity, preserve_observer_frame)
609
610 # Apply transformation to data
--> 611 new_data = _apply_relativistic_doppler_shift(self, fin_obs_vel - init_obs_vel)
612
613 new_coord = self.replicate(value=new_data, observer=observer)

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/coordinates/spectral_coordinate.py in _apply_relativistic_doppler_shift(scoord, velocity)
62 return squantity / doppler_factor
63 elif squantity.unit.is_equivalent(KMS): # velocity
---> 64 return (squantity.to(u.Hz) / doppler_factor).to(squantity.unit)
65 else: # pragma: no cover
66 raise RuntimeError(f"Unexpected units in velocity shift: {squantity.unit}. "

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/units/decorators.py in wrapper(*func_args, **func_kwargs)
300 # Call the original function with any equivalencies in force.
301 with add_enabled_equivalencies(self.equivalencies):
--> 302 return_ = wrapped_function(*func_args, **func_kwargs)
303
304 # Return

/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/coordinates/spectral_quantity.py in to(self, unit, equivalencies, doppler_rest, doppler_convention)
277
278 if doppler_convention is None:
--> 279 raise ValueError("doppler_convention not set, cannot convert to/from velocities")
280
281 if doppler_rest is None:

ValueError: doppler_convention not set, cannot convert to/from velocities
```

When I drill in with the debugger, I see:
```python
> /orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/wcs/wcsapi/fitswcs.py(602)value_from_spectralcoord()
600 return spectralcoord.to_value(**kwargs)
601 else:
--> 602 return spectralcoord.with_observer_stationary_relative_to(observer).to_value(**kwargs)
603
604 classes['spectral'] = (u.Quantity, (), {}, spectralcoord_from_value)

ipdb> spectralcoord

target:
observer to target (computed from above):
radial_velocity=5.37398150550586 km / s
redshift=1.792583345272547e-05)
[1970.36615651, 1970.36615651, 1970.36615651, ..., 1970.36615651,
1970.36615651, 1970.36615651] m / s>
ipdb> observer

```

I don't know where this coordinate (293.8 -21.6) is coming from; none of my data or headers are anywhere near that location. That's specifying the "observer" in the SpectralCoord. I think that's the source of this error?

Note that the final error message is extremely confusing, since both headers explicitly define the velocity convention (`CTYPE3 = 'VRAD'`).

Further, this is happening at the step in the loop when celestial coordinates are being transformed, not velocity coordinates, so it shouldn't even care about the spectral dimension.

This is new behavior in the last ~6 months.

Reproduce the traceback from the reported cube reprojection, then inspect reproject/interpolation/core.py and reproject/wcs_utils.py, especially _reproject_full and efficient_pixel_to_pixel. Trace how celestial-coordinate transformations pass spectral values into Astropy WCS, and verify that the reprojection no longer fails on the reported velocity-coordinate cubes.

Written by the indexing model from the issue text.