I need to plot some SDSS images with python. I downloaded the original frames from the SDSS DR7 and I'm trying to display (a zoom of) the image with the WCS coordinates projection. The original image has north on the right and east on top. According to DS9 the alignment is not perfect. I would like to have the final image with north on top and east left. My attempt is this one:
import numpy as np
import matplotlib.pyplot as plt
import astropy.visualization as vis
import astropy.units as un
from astropy.io import fits
from astropy.wcs import WCS
from astropy.nddata.utils import Cutout2D
def open_image(imagename):
hdu = fits.open(imagename)
image = hdu[0].data
header = hdu[0].header
wcs = WCS(header)
return image, header, wcs
def plot_image(imagename):
image, header, wcs = open_image(imagename) #loading image
print(wcs)
#selecting a region centered on the target
zoom = Cutout2D(image,(412,559),51*un.pixel, wcs = wcs, copy=True)
#setting axis format
fig,ax = plt.subplots(1,1, subplot_kw=dict(projection=zoom.wcs))
ra = ax.coords[0]
dec = ax.coords[1]
ra.set_major_formatter('hh:mm:ss.s')
dec.set_major_formatter('dd:mm:ss')
#setting image scale
interval = vis.PercentileInterval(99.9)
vmin,vmax = interval.get_limits(zoom.data)
norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))
ax.imshow(zoom.data, cmap =plt.cm.Reds, norm = norm, origin = 'lower')
ax.set_ylabel('Dec.')
ax.set_xlabel('RA')
plt.show()
if __name__ == '__main__':
imagename = '../SDSS_g_orig.fits'
plot_image(imagename)
You can find the results in attachment. The coordinates are not present on the axes. The image also kept its original orientation. I assumed that the WCS projection automatically rotates the image to align it, but evidently it is not true.
How can I align the image with north on top and east left?
P.s.: this is the output of the print(wcs) line:
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 195.77014030000001 16.47383778
CRPIX : 1024.5 744.5
CD1_1 CD1_2 : 6.1883041204266702e-06 0.000109834787851833
CD2_1 CD2_2 : 0.000109820625 -6.2125672043002999e-06
NAXIS : 2048 1489
EDIT: I tried to use the SkyView module but the result is not satisfactory. The intensity contours of the images (image + contours is my final goal) are different and the skyview image shows some kind of artifact right in the middle of the target.
With the help of the Python users in Astronomy Facebook group I found an answer. It requires the installation of Montage and montage-wrapper but it allows to rotate the image and to continue work with astropy and matplotlib.
I just had to modify the open_image() function in this way:
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
import montage_wrapper as montage
def open_image(imagename):
hdu = fits.open(imagename)
hdu = montage.reproject_hdu(hdu[0], north_aligned=True)
image = hdu.data
nans = np.isnan(image)
image[nans] = 0
header = hdu.header
wcs = WCS(header)
return image, header, wcs
The resulting image will be rotated in such a way that north is up. The remaining part of the script does not need to be modified. Only the Cutout2D object must be recentered on the target.
I'm attaching the final image with the contours.
