astropy: WCSaxes show data flipped (labels are attached to wrong axes)
Minimalist WE:
from astropy.wcs import WCS
from astropy.io import fits
from astropy.visualization import simple_norm
import pylab as pl
from astropy import units as u
hdul = fits.open('https://stpubdata-jwst.stsci.edu/ero/jw02731/L3/t/jw02731-o002_t017_miri_f1800w_i2d.fits')
pl.figure(figsize=(10,10))
ww = WCS(hdul[1].header)
ax = pl.subplot(projection=ww)
ax.imshow(hdul[1].data, norm=simple_norm(hdul[1].data, max_percent=99.95))
lon = ax.coords[0]
lat = ax.coords[1]
lon.set_major_formatter('dd:mm:ss.ss')
lat.set_major_formatter('dd:mm:ss.ss')
lon.set_ticks(spacing=30. * u.arcsec)
lat.set_ticks(spacing=30. * u.arcsec)
lon.set_ticklabel(rotation=45, pad=30)
Result:
Problem: the X-axis is declination!
The header is this:
XTENSION= 'IMAGE ' / Image extension
BITPIX = -32 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 3470
NAXIS2 = 1174
PCOUNT = 0 / number of parameters
GCOUNT = 1 / number of groups
EXTNAME = 'SCI ' / extension name
MJD-BEG = 59741.12639981297 / [d] exposure start time in MJD
MJD-AVG = 59741.34160363826 / [d] exposure mid-point in MJD
MJD-END = 59741.55679824074 / [d] exposure end time in MJD
TDB-BEG = 59741.13163161634 / [d] TDB time of exposure start in MJD
TDB-MID = 59741.34683091205 / [d] TDB time of exposure mid-point in MJD
TDB-END = 59741.56202095382 / [d] TDB time of exposure end in MJD
XPOSURE = 5994.08 / [s] Effective exposure time
TELAPSE = 5994.08 / [s] Total elapsed exposure time
JWST ephemeris information
REFFRAME= 'EME2000 ' / Ephemeris reference frame
EPH_TYPE= 'Definitive' / Definitive or Predicted
EPH_TIME= 59741.14305555556 / UTC time of position and velocity vectors in ep
JWST_X = -27666471.75619852 / [km] barycentric JWST X coordinate at MJD_AVG
JWST_Y = -138329023.1183411 / [km] barycentric JWST Y coordinate at MJD_AVG
JWST_Z = -60289095.90038356 / [km] barycentric JWST Z coordinate at MJD_AVG
OBSGEO-X= 270741107.153853 / [m] geocentric JWST X coordinate at MJD_AVG
OBSGEO-Y= -1382637909.38873 / [m] geocentric JWST Y coordinate at MJD_AVG
OBSGEO-Z= -955370885.732361 / [m] geocentric JWST Z coordinate at MJD_AVG
JWST_DX = 28.95847034488103 / [km/s] barycentric JWST X velocity at MJD_AVG
JWST_DY = -4.858839776968531 / [km/s] barycentric JWST Y velocity at MJD_AVG
JWST_DZ = -2.234909903693477 / [km/s] barycentric JWST Z velocity at MJD_AVG
OBSGEODX= 125.020158763399 / [m/s] geocentric JWST X velocity at MJD_AVG
OBSGEODY= 35.8482128539545 / [m/s] geocentric JWST Y velocity at MJD_AVG
OBSGEODZ= -114.373851163562 / [m/s] geocentric JWST Z velocity at MJD_AVG
PA_APER = 116.7940399545248 / [deg] Position angle of aperture used
VA_SCALE= 0.9999562813367563 / Velocity aberration scale factor
BUNIT = 'MJy/sr ' / physical units of the array values
Photometry information
PHOTMJSR= 0.4943844079971313 / Flux density (MJy/steradian) producing 1 cps
PHOTUJA2= 11.6202239773816 / Flux density (uJy/arcsec2) producing 1 cps
PIXAR_SR= 2.84403609523084E-13 / Nominal pixel area in steradians
PIXAR_A2= 0.0121 / Nominal pixel area in arcsec^2
Information about the coordinates in the file
RADESYS = 'ICRS ' / Name of the coordinate reference frame
Spacecraft pointing information
RA_V1 = 159.3540771893557 / [deg] RA of telescope V1 axis
DEC_V1 = -58.73506412600148 / [deg] Dec of telescope V1 axis
PA_V3 = 111.8777236761156 / [deg] Position angle of telescope V3 axis
WCS parameters
WCSAXES = 2 / number of World Coordinate System axes
CRPIX1 = 1731.103892447788 / axis 1 coordinate of the reference pixel
CRPIX2 = 581.0376680959058 / axis 2 coordinate of the reference pixel
CRVAL1 = 159.221073504818 / first axis value at the reference pixel
CRVAL2 = -58.61772542066391 / second axis value at the reference pixel
CTYPE1 = 'RA---TAN' / first axis coordinate type
CTYPE2 = 'DEC--TAN' / second axis coordinate type
CUNIT1 = 'deg ' / first axis units
CUNIT2 = 'deg ' / second axis units
CDELT1 = 3.0808700908093E-05 / first axis increment per pixel
CDELT2 = 3.0808700908093E-05 / second axis increment per pixel
PC1_1 = 0.4507846893158868 / linear transformation matrix element
PC1_2 = 0.8926327149944592 / linear transformation matrix element
PC2_1 = 0.8926327149944592 / linear transformation matrix element
PC2_2 = -0.4507846893158868 / linear transformation matrix element
S_REGION= 'POLYGON ICRS 159.144213298 -58.657226232 159.206276111 &'
CONTINUE '-58.673552892 159.298755295 -58.578102372 159.236818600 &'
CONTINUE '-58.561820231&'
CONTINUE '' / spatial extent of the observation
VELOSYS = 13087.39 / [m/s] Barycentric correction to radial velocity
EXTVER = 1 / extension value
Note that the PC matrix is bigger in the 1_2 and 2_1 terms - so, the image is rotated. But, the axis labels are attached as if the image is not rotated.
Workaround: Reproject the image to a north-facing WCS. There’s nothing wrong with the underlying axes, just the labels, afaict.
About this issue
- Original URL
- State: open
- Created 2 years ago
- Reactions: 1
- Comments: 23 (22 by maintainers)
That works:
Yes looking at this more I think this is what we might want to use to figure out which coordinate to show preferentially on which axis. We’ll have to think about whether we are ok with this potentially changing user figures - but I think it might be ok because for cases where we would actually swap the coordinates, the default output was probably not optimal and users might have anyway manually specified where to show the labels (in which case changing the default won’t change anything). If users didn’t fix it, they might welcome the improvement.
We already sort of do in that if there’s more than one world coord for an axis we use it? I see no reason why we could make the default behaviour auto detect based on the derivatives.
I agree this isn’t a bug, but it is definitely unexpected behavior and I think we’ll be seeing this often with JW data
I don’t think simply swapping the axes fixes it. The WCS itself is correct, it’s just that wcsaxes is labeling the wrong one. The wcsaxes default made sense for virtually every telescope before JW started scanning sideways, so I think this is a relatively simple edge case to catch.