Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pixel-to-pixel conversion failing because of doppler convention #312

Open
keflavich opened this issue Sep 24, 2022 · 0 comments
Open

Pixel-to-pixel conversion failing because of doppler convention #312

keflavich opened this issue Sep 24, 2022 · 0 comments

Comments

@keflavich
Copy link
Contributor

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:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/scratch/local/47745774/ipykernel_15582/2515218747.py in <module>
      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:

> /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
<SpectralCoord 
   (observer: <LSRK Coordinate: (ra, dec, distance) in (deg, deg, m)
                  (293.81585318, -21.63372113, 1.50879416e+11)
               (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                  (0., 0., 0.)>
    target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                (254.57170833, -42.86875, 1000.)
             (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                (0., 0., -0.)>
    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
<LSRK Coordinate: (ra, dec, distance) in (deg, deg, m)
    (293.81585318, -21.63372113, 1.50879416e+11)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (0., 0., 0.)>

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant