I am trying to reproject an observation of the Sun in ICRS coordinates on to an observation in Heliographic coordinates. The header of the input map is:
Input map header
(‘SIMPLE’, True),
(‘BITPIX’, -32),
(‘NAXIS’, 4),
(‘NAXIS1’, 556),
(‘NAXIS2’, 516),
(‘EXTEND’, True),
(‘BLOCKED’, True),
(‘OBJECT’, ‘Apr3’),
(‘TELESCOP’, ‘EVLA’),
(‘INSTRUME’, ‘EVLA’),
(‘OBSERVER’, ‘Dr. Debo’),
(‘DATE-OBS’, ‘2020-04-03’),
(‘DATE-MAP’, ‘2021-10-07’),
(‘BSCALE’, 1.0),
(‘BZERO’, 0.0),
(‘BUNIT’, ‘KELVIN’),
(‘EQUINOX’, 2000.0),
(‘VELREF’, 257),
(‘ALTRPIX’, 1.0),
(‘OBSRA’, 12.9474625144),
(‘OBSDEC’, 5.54566945688),
(‘NUMDOCAL’, 0),
(‘NUMDOBAN’, 0),
(‘NUMDOPOL’, 0),
(‘DATAMAX’, 271617.4688),
(‘DATAMIN’, -8853.824219),
(‘CTYPE1’, ‘RA—SIN’),
(‘CRVAL1’, 12.9474625144),
(‘CDELT1’, -0.0004166666768),
(‘CRPIX1’, 537.0),
(‘CROTA1’, 0.0),
(‘CTYPE2’, ‘DEC–SIN’),
(‘CRVAL2’, 5.54566945688),
(‘CDELT2’, 0.0004166666768),
(‘CRPIX2’, -228.0),
(‘CROTA2’, 0.0)
Trying to do the reprojection on to a Heliographic system doesn’t work however, resulting in the error:
UnitsError: The input ICRS coordinates do not have length units. This probably means you created coordinates with lat/lon but no distance. Heliocentric<->ICRS transforms cannot function in this case because there is an origin shift.
I think I have reduced the issue to the following short, self contained example:
from sunpy.coordinates import Helioprojective
from astropy.coordinates import ICRS, SkyCoord
import astropy.units as u
# Create a helioprojective coordinate without an input distance; sunpy
# will implicitly put this coordinate on the solar disk when transforming
helio_frame = Helioprojective(observer='earth', obstime='2021-01-01')
c = SkyCoord(0*u.deg, 0*u.deg, frame=helio_frame)
# Transform to ICRS, and drop the distance
c_icrs = c.transform_to('icrs')
c_icrs = SkyCoord(c_icrs.ra, c_icrs.dec, frame='icrs')
# Try transforming back to the Helioprojective frame. This errors
c_icrs.transform_to(helio_frame)
I think what I want to happen is for the transformation from ICRS > Helioprojective to implicitly assume my coordinate is on the surface of the Sun, in a similar way to how the Helioprojective > ICRS transformation works. Is this a missing feature in sunpy’s coordinate transformations? Or is there a way I can manually force the solar surface condition when doing this transformation?