3.1 - Displaying a coordinate layer via the Header

In [6]:
%matplotlib inline
from astropy.wcs import WCS
import numpy as np
from matplotlib.colors import LogNorm
import warnings
warnings.filterwarnings("ignore", category=UserWarning)#delete after also
#warnings.filterwarnings("ignore", category=Warning)#delete after also

#Create WCS object
wcs_for_plot = WCS(image_header)

#prepare plot
fig5 = plt.figure(figsize=(20, 20))
ax = plt.subplot(projection=wcs_for_plot) #add wcs informations on plot
plt.imshow(image_data, origin='lower', cmap='gray', aspect='equal',norm=LogNorm())
plot_title = "Gamma Cassiopae "
plt.title(plot_title ,fontsize=16, loc='left')

#add and configure axis informations
ax = plt.gca()
ax.grid(True, color='white', ls='solid')
ra = ax.coords[0]
dec = ax.coords[1]

ra.set_axislabel('RA (J2000)',  minpad=1)
dec.set_axislabel('DEC (J2000)',  minpad=1)
ra.set_major_formatter('hh:mm:ss.s')
dec.set_major_formatter('dd:mm:ss.s')

#add other grid as overlay
overlay = ax.get_coords_overlay('icrs') # or galactic, alt-az, ...
overlay.grid(True, color='lime', ls='solid')
overlay[0].set_axislabel('RA (deg)')
overlay[1].set_axislabel('DEC (deg)')

#uncomment the line after for show contour on objects
#ax.contour(image_data, transform=ax.get_transform(wcs_for_plot), levels=np.logspace(-3, 0, 3), colors='orange', alpha=1)

#show
plt.show()

3.2 - Retrieval of coordinates by the astrometry.net service

To be able to retrieve the coordinates corresponding to the image, a request to the astrometry service (astrometry.net) is made through the AstroQuery library.

Converting the FITS image to JPEG

As the initial image is 80 MB, it is possible to convert this FITS image to jpeg/png, before using the astrometry service, in order to have a lighter image of about 100kb and speed up the process.

In [7]:
#image name (FITS) to convert : image_path
#image name convert in jpeg/png : gamcas_jpeg.jpeg
from matplotlib.colors import LogNorm
%matplotlib inline
fig6 = plt.figure(figsize=(16, 9))
plt.imshow(image_data, cmap='gray', norm=LogNorm())
plt.grid()
plt.axis('off')
plt.savefig("dataset/gamcas_jpeg.jpeg", bbox_inches='tight')

Astrometry.net API Query

It is necessary to specify the API key linked to your account. This can be found on your account page at this address: http://nova.astrometry.net/api_help. The present warning is normal and requires you to specify our configuration.

In [9]:
import warnings
warnings.filterwarnings("ignore", category=Warning)

#Import astroquery/astrometry library and World Coordinate System module from astropy
from astroquery.astrometry_net import AstrometryNet
from astropy.wcs import WCS


ast = AstrometryNet()
ast.api_key = '##############' # PUT YOUR API KEY HERE

try_again = True
submission_id = None

while try_again:
    try:
        if not submission_id:
            wcs_header = ast.solve_from_image('dataset/gamcas_jpeg.jpeg',
                                              submission_id=submission_id)
        else:
            wcs_header = ast.monitor_submission(submission_id,
                                                solve_timeout=120)
    except TimeoutError as e:
        submission_id = e.args[1]
    else:
        # got a result, so terminate
        try_again = False


#when solve is ok, record data in wcs_gamcas object
if wcs_header:
    print('Solve OK - Print wcs result ...')
    wcs_gamcas = WCS(wcs_header)
    print(wcs_gamcas)
    
else:
    print('nok')
    

WARNING: Astrometry.net API key not found in configuration file [astroquery.astrometry_net.core]
WARNING: You need to manually edit the configuration file and add it [astroquery.astrometry_net.core]
WARNING: You may also register it for this session with AstrometryNet.key = 'XXXXXXXX' [astroquery.astrometry_net.core]

Solving.............................................................Solve OK - Print wcs result ...
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  
CRVAL : 13.8500869723  60.6227667979  
CRPIX : 272.770778656  185.717145443  
CD1_1 CD1_2  : 0.00069678872339  0.00126336322432  
CD2_1 CD2_2  : 0.0012606991432  -0.000696346507525  
NAXIS : 0  0

Display of celestial coordinates after return from the Astrometry.net API

The coordinates can then be displayed over the image in the desired coordinate system. More information about the Astropy WCS management system at these addresses :

In [10]:
import matplotlib.image as mpimg
    
#Prepare plot
img = mpimg.imread('dataset/gamcas_jpeg.jpeg') #Image name_file
fig7 = plt.figure(figsize=(20, 20))

ax = plt.subplot(projection=wcs_gamcas) #Add wcs projection with wcs_gamcas file
plt.imshow(img, origin='lower', cmap='gray', aspect='equal',norm=LogNorm())

plt.title('gam Cas - Astrometry.net API',fontsize=16, loc='center')

ax.coords['ra'].set_axislabel('RA')
ax.coords['dec'].set_axislabel('DEC')

overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')

plt.show()

Retrieving and displaying coordinates without going through the API in line of code

In the event that astrometric tracking is not performed by the Python API via the previous piece of code. Typical use via the FITS file upload on astrometry.net will retrieve the single WCS file, or the FITS image including the WCS header. The code below shows the display of coordinates on the image using this method. The file retrieved from astrometry.net is called 'wcs.fit'

In [12]:
#Get wcs file path
wcs_file = fits.open('dataset/wcs.fits')
header_wcs_file = wcs_file[0].header

#Create WCS object
wcs_for_plot2 = WCS(header_wcs_file)

#Image plot must be the same than upload image on astrometry.net
import matplotlib.image as mpimg
img = mpimg.imread('dataset/gamcas_jpeg.jpeg')

fig8 = plt.figure(figsize=(20, 20))
ax = plt.subplot(projection=wcs_for_plot2)
plt.imshow(img, origin='lower', cmap='gray', aspect='equal',norm=LogNorm())

plot_title = "Gamma Cassiopae"
plt.title(plot_title ,fontsize=16, loc='left')

ax.coords['ra'].set_axislabel('RA')
ax.coords['dec'].set_axislabel('DEC')

overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')

plt.show()