Plotting SDSS images with python

Question

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

output image

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.

Original with contours Skyview with contours

Cloud182 · Accepted Answer · 2018-08-22T13:46:02.717

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.

final+contours

FWIW to fill the `NaN`s you can also use [`numpy.nan_to_num(image, copy=False)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nan_to_num.html) to modify the array in-place. — Iguananaut, Aug 22 '18 at 13:13
Or use `astropy.convolution.interpolate_replace_nans`: http://docs.astropy.org/en/stable/convolution/index.html#using-astropy-s-convolution-to-replace-bad-data — Thriveth, Sep 17 '18 at 13:27

Brett Morris · Answer 2 · 2018-08-21T16:57:05.467

An easier way might be to use astroquery's SkyView module. For example:

import matplotlib.pyplot as plt
from astroquery.skyview import SkyView
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

# Query for SDSS g images centered on target name
hdu = SkyView.get_images("M13", survey='SDSSg')[0][0]

# Tell matplotlib how to plot WCS axes
wcs = WCS(hdu.header)
ax = plt.gca(projection=wcs)

# Plot the image
ax.imshow(hdu.data)
ax.set(xlabel="RA", ylabel="Dec")
plt.show()

which produces an image like this:

You can replace "M13" with your target name, or coordinates -- see the SkyView.get_images docs for all of the finer knobs you can tune.

Then, if you want to flip the directions of either axis, you can call

ax.invert_xaxis()

or

ax.invert_yaxis()

This is the first thing I tried. However I found that the image downloaded in this way shows an artifact I do not like. I'm editing the original message with the result of the skyview module. However, thanks for the answer. — Cloud182, Aug 22 '18 at 08:01

Plotting SDSS images with python

2 Answers2