I’m looking for guidance on the best way to have GCRS coordinates interact with WCS objects, and whether I’m handling coordinate frames correctly. For the PUNCH mission (a set of four wide-field imager satellites observing the heliosphere), we’re using the stars to align our images. To run our starfield-subtraction step, we need alignment to within ~0.1 px (0.0025 deg or 9 arcsec). I’m calculating that, across our 90-degree composite field-of-view, relativistic aberration due to Earth’s orbital motion will have an effect that varies spatially from ~0.15 to 0.2 px, and aberration due to our spacecraft’s orbital motion around Earth will be up to 0.05 px, so we’re looking to handle these effects, and their variability across our field-of-view, correctly.
When we run our alignment process, we identify stars in our images and compare them to the Gaia catalog, and we iterate our WCS until the star locations agree. As I understand it, the Gaia coordinates are in ICRS, and if I use Astropy to convert those coordinates to GCRS, that corrects the stellar coordinates for aberration, giving the coordinates at which we can expect to see the stars at a given time.
I’ve been struggling to find the best way to transform those coordinates to GCRS and then project them into our image plane, as it seems WCS objects want to accept and produce SkyCoords in ICRS (or FK*, etc.), but I haven’t been able to find a way to tell a WCS it should accept/produce GCRS coordinates. If I have a FITS header with these WCS keys:
WCSAXESA= 2 /Number of coordinate axes
CRPIX1A = 1024.5 /Pixel coordinate of reference point
CRPIX2A = 1024.5 /Pixel coordinate of reference point
PC1_1A = 0.99250327546037 /Coordinate transformation matrix element
PC1_2A = -0.12221803549573 /Coordinate transformation matrix element
PC2_1A = 0.12221803549573 /Coordinate transformation matrix element
PC2_2A = 0.99250327546037 /Coordinate transformation matrix element
CDELT1A = -0.02444444444 /[deg] Coordinate increment at reference point
CDELT2A = 0.02444444444 /[deg] Coordinate increment at reference point
CUNIT1A = 'deg' /Units of coordinate increment and value
CUNIT2A = 'deg' /Units of coordinate increment and value
CTYPE1A = 'RA---AZP' / Right ascension, zenithal/azimuthal perspective
CTYPE2A = 'DEC--AZP' / Declination, zenithal/azimuthal perspective pro
CRVAL1A = 192.29073737789 /[deg] Coordinate value at reference point
CRVAL2A = 19.754133770555 /[deg] Coordinate value at reference point
PV2_1A = 0.0 /AZP projection parameter
LONPOLEA= 180.0 /[deg] Native longitude of celestial pole
LATPOLEA= 19.754133770555 /[deg] Native latitude of celestial pole
TIMESYS = 'UTC' /Time scale
MJDREF = 0.0 /[d] MJD of fiducial time
DATE-OBS= '2025-10-03T00:00:29.244' /ISO-8601 time of observation
MJD-OBS = 60951.000338472 /[d] MJD of observation
DATE-BEG= '2025-10-03T00:00:17.244' /ISO-8601 time at start of observation
MJD-BEG = 60951.000199583 /[d] MJD at start of observation
DATE-AVG= '2025-10-03T00:00:29.244' /ISO-8601 time at midpoint of observation
MJD-AVG = 60951.000338472 /[d] MJD at midpoint of observation
DATE-END= '2025-10-03T00:00:41.244' /ISO-8601 time at end of observation
MJD-END = 60951.000477361 /[d] MJD at end of observation
RADESYSA= 'ICRS' /Equatorial coordinate system
RSUN_REF= 695700000.0 /[m] Solar radius
DSUN_OBS= 149699314230.95 /[m] Distance from centre of Sun to observer
CRLN_OBS= 45.530540253641 /[deg] Carrington heliographic lng of observer
CRLT_OBS= 6.6274351073017 /[deg] Heliographic latitude of observer
HGLN_OBS= 0.0007774395069309 /[deg] Stonyhurst heliographic lng of observer
HGLT_OBS= 6.6274351073017 / [deg] Heliographic latitude of observer
and if I have a SkyCoord in GRCS (containing the CRVAL coordinate):
c = SkyCoord(192.29073738*u.deg, 19.75413377*u.deg, frame='gcrs', obstime='2025-10-03T00:00:29.244')
then it appears that world_to_pixel converts this GCRS coordinate to ICRS and then to pixel coordinates:
wcs.world_to_pixel(c)
>>> (array(1023.27694139), array(1023.44011703))
wcs.world_to_pixel(c.transform_to('icrs'))
>>> (array(1023.27694139), array(1023.44011703))
Whereas if I put the same RA and dec into an ICRS SkyCoord, I get a different result (it comes out to CRPIX):
c = SkyCoord(192.29073738*u.deg, 19.75413377*u.deg, frame='icrs', obstime='2025-10-03T00:00:29.244')
wcs.world_to_pixel(c)
>>> (array(1023.49999992), array(1023.49999999))
And pixel_to_world produces ICRS SkyCoords. So the WCS seems to accept and produce ICRS coordinates. I haven’t been able to find a way to tell the WCS it should accept/produce GCRS coordinates instead, so that I can give it an ICRS SkyCoord containing a star’s Gaia coordinates, and get the pixel location at which I can expect to see that star given aberration. The best approach I’ve been able to put together is to take my ICRS catalog locations, transform them to GCRS, and then re-package them as a SkyCoord that says it’s ICRS, so that the WCS will project the coordinates without further transformation:
def prep_star_coords(stars_in_image: pd.DataFrame, image_header: NormalizedMetadata) -> SkyCoord:
"""
Convert ICRS coordinates to GCRS and put in a SkyCoord that says it's ICRS.
"""
# Convert stellar coordinates to GCRS centered on the spacecraft location
sc_location = EarthLocation.from_geodetic(lon=image_header["GEOD_LON"].value * u.deg,
lat=image_header["GEOD_LAT"].value * u.deg,
height=image_header["GEOD_ALT"].value * u.m)
geoloc, geovel = sc_location.get_gcrs_posvel(image_header.astropy_time)
catalog_stars = SkyCoord(
np.array(stars_in_image["RAdeg"]) * u.degree,
np.array(stars_in_image["DEdeg"]) * u.degree,
np.array(stars_in_image["Dist_ly"]) * u.lyr,
frame="icrs", obsgeoloc=geoloc, obsgeovel=geovel, obstime=image_header.astropy_time,
).transform_to("gcrs")
return SkyCoord(catalog_stars.ra, catalog_stars.dec, frame="icrs")
Is there a better way to do this? I’ve also thought about making a class that wraps my WCS and re-packages coordinates to/from GCRS in the overridden pixel_to_world and world_to_pixel.
I also only think I’m understanding these coordinate frames and their differences correctly—I’d appreciate comment there. Thank you!