[AstroPy] how to display an image on an different WCS

[Hessman, Frederic]

There are tutorials and plenty of examples showing how to put coordinates onto a plot of an image given the image's WCS, but what about the other way around:  define the coordinate system of the plot and then place the image in that frame?

Let's say one has a 101x101 image with a flat blob on it and an outer frame (so one can easily see rotations) and whose pixels are in an unrotated FK5 frame, centered at (00:11:22,-33:44:55) : the accompanying script will demo all of this.

I've got a WCS for the image, called "wcs"

CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 2.8416666666666663 33.74861111111111
CRPIX : 51.0 51.0
PC1_1 PC1_2  : 1.0 0.0
PC2_1 PC2_2  : 0.0 1.0
CDELT : 0.000277777777777777 0.000277777777777777
NAXIS : 101  101

I want to display the image where the galactic coordinates are the main axes, so I've created a galactic WCS from scratch at the same central position, called "gal_wcs".

CTYPE : 'GLON-TAN' 'GLAT-TAN'
CRVAL : 113.46911966371806 -28.38577360534605
CRPIX : 51.0 51.0
PC1_1 PC1_2  : 1.0 0.0
PC2_1 PC2_2  : 0.0 1.0
CDELT : 0.000277777777777777 0.000277777777777777
NAXIS : 101  101

Naively, I'd expect to be able to do it this way :

fig = plt.figure(figsize=(4,4))

ax = plt.subplot (projection=gal_wcs, label='galactic')
ax.coords.grid (True, color='gray',ls='solid')
ax.coords[0].set_axislabel (r'Galactic Longitude')
ax.coords[1].set_axislabel (r'Galactic Latitude')

overlay = ax.get_coords_overlay ('fk5')
overlay.grid (True, color='gray',ls='dotted')
overlay[0].set_axislabel (r'Right Ascension (J2000)')
overlay[1].set_axislabel (r'Declination (J2000)')

fk5_2_gal = ax.get_transform(wcs) # TRANSFORMATION FROM wcs=FK5 TO ax=galactic
ax.imshow (data,origin='lower',interpolation='nearest',transform=fk5_to_gal)

but this doesn't work: the image contents are rotated but matplotlib creates data outside of the frame to maintain an un-rotated area!  😳

In order to see what's happening, one can increase the size of the coordinate space, e.g.

plt.xlim (right=0-0.2*nx,left=nx-1+0.2*nx)
plt.ylim (top=0-0.2*ny,bottom=ny-1+0.2*ny)

but this doesn't really change anything: one sees the displayed image is not rotated, just the contents, and one sees the area of data created by matplotlib outside of the rotated image.  The area used is not the same in the two cases, so it's not simply matplotlib filling in the area defined by gal_wcs.

Where am I confused?  Can someone suggest an example I've missed?  Running python 3.9.11, astropy 6.0.0, matplotlib 3.5.1 on Darwin 12.7.4 / Monterey.

Rick



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20240423/a8c38ab7/attachment-0001.html>
-------------- next part --------------
# test.py

"""
Plots a fake (RA,DEC) image in a galactic reference frame.
"""

import astropy.units      as u
import matplotlib.pyplot  as plt
import numpy              as np

from astropy.coordinates import FK5, Galactic, SkyCoord
from astropy.wcs         import WCS, utils

# ---- DEFINE CENTER OF IMAGE

center = SkyCoord ('00:11:22','+33:44:55',frame='icrs',unit=(u.hourangle,u.deg))
print ('Center:',center)

# ---- CREATE FAKE IMAGE

data = np.zeros((101,101))
ny,nx = data.shape
x,y = np.meshgrid (np.arange(nx),np.arange(ny))
x0,y0 = nx//2+1,ny//2+1
data = np.exp(-((x-x0)**2/20**2+(y-y0)**2/3**2))
mask = (x == 0)+(x == nx-1)+(y == 0//3)+(y == ny-1)
data[np.where(mask)] = 1.0
plt.imshow (data)
plt.show ()

hdr = utils.celestial_frame_to_wcs (FK5()).to_header()
hdr['CRPIX1'] = nx//2+1
hdr['CRPIX2'] = ny//2+1
hdr['CDELT1'] = 1./3600	# ARCSEC
hdr['CDELT2'] = 1./3600 # ARCSEC
hdr['CRVAL1'] = center.ra.deg
hdr['CRVAL2'] = center.dec.deg
hdr['NAXIS']  = 2
hdr['NAXIS1'] = nx
hdr['NAXIS2'] = ny

"""
hdr['PC1_1'] = np.cos(45*np.pi/180)
hdr['PC1_2'] = np.sin(45*np.pi/180)
hdr['PC2_1'] = hdr['PC1_2']
hdr['PC2_2'] = hdr['PC1_1']
"""
wcs  = WCS(hdr)
print ('\nWCS of image:\n',wcs)

# ---- CREATE A CORRESPONDING GALACTIC WCS

adjacent = utils.pixel_to_skycoord (nx//2,ny//2,wcs)
cdelt = np.abs(center.separation (adjacent).to(u.deg))

gal_hdr = utils.celestial_frame_to_wcs (Galactic()).to_header()
gal_hdr['CRPIX1'] = nx//2+1
gal_hdr['CRPIX2'] = ny//2+1
gal_hdr['CDELT1'] = hdr['CDELT1']
gal_hdr['CDELT2'] = hdr['CDELT2']
gal_hdr['CRVAL1'] = center.galactic.l.value
gal_hdr['CRVAL2'] = center.galactic.b.value
gal_hdr['NAXIS']  = 2
gal_hdr['NAXIS1'] = nx
gal_hdr['NAXIS2'] = ny
gal_wcs = WCS(gal_hdr)
print ('\nWCS of equivalent galactic frame:\n',gal_wcs)

# ---- PLOT

fig = plt.figure(figsize=(4,4))

ax = plt.subplot (projection=gal_wcs, label='galactic')
ax.coords.grid (True, color='gray',ls='solid')
ax.coords[0].set_axislabel (r'Galactic Longitude')
ax.coords[1].set_axislabel (r'Galactic Latitude')

overlay = ax.get_coords_overlay ('fk5')
overlay.grid (True, color='gray',ls='dotted')
overlay[0].set_axislabel (r'Right Ascension (J2000)')
overlay[1].set_axislabel (r'Declination (J2000)')

fk5gal = ax.get_transform(wcs)
ax.imshow (data,origin='lower',interpolation='nearest',transform=fk5gal)

# ---- SHOW RESULT

plt.xlim (right=0-0.2*nx,left=nx-1+0.2*nx)
plt.ylim (top=0-0.2*ny,bottom=ny-1+0.2*ny)
plt.show ()