How to superimpose fits images using python well in WCS space

Environment

OK if matplotlib and astropy are installed in python. If anconda3 system is included, it is OK by default.

An example of overlapping two images with Galactic coordinates

Taking the visible light and X-ray NGC 4051 images as examples, the plot of the whole image and the plot of the enlarged image near the center are shown.

point

--Just generate a WCS class from the image fits file and put it in the matplotlib projection option. However, with this alone, it is displayed only in the initial image size.

--If you want to enlarge a certain place, cut the WCS object in the same way as the image file. The part of the sample code, hriwcs [hricx-hridx: hricx + hridx, hricy-hridx: hricy + hridx]. Give this to matplotlib's wcs option.

――If you want to display two images on top of each other, it becomes a difficult problem. In that case, it is briefly summarized in How to merge images with WCS coordinates and plot with python.

Code and file storage

qiita_plot_fits_2image.py --Visible and X-ray image of NGC 4051
- http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/NGC4051_DSS.fits
- http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/NGC4051_HRI.fits

For how to take an image, How to search using python's astroquery and get fits image with skyviewを参考にされたい。

Sample code

`qiita_plot_fits_2image.py`



import matplotlib 
matplotlib.use('TkAgg')
import glob 
import matplotlib.pyplot as plt
import sys
from matplotlib.colors import LogNorm
from astropy.wcs import WCS
from astropy.io import fits

name="NGC4051"

hri = glob.glob(name + '*HRI*.fits')
dss = glob.glob(name + '*DSS*.fits')	

if len(hri) == 1 and len(dss) == 1:
	pass
else:
	print("need to store HRI and DSS fits")

F = plt.figure(figsize=(10,8))

hrifilename = hri[0]
hriname = hrifilename.replace(".fits",".png ")
hrihdu = fits.open(hrifilename)[0]
hriwcs = WCS(hrihdu.header)
hridata = hrihdu.data

hrixlen, hriylen = hridata.shape
hricx = int(0.5 * hrixlen)
hricy = int(0.5 * hriylen)
hridx = int(hrixlen*0.1)
hriwcscut = hriwcs[hricx-hridx:hricx+hridx,hricy-hridx:hricy+hridx]

dssfilename = dss[0]
dssname = dssfilename.replace(".fits",".png ")
dsshdu = fits.open(dssfilename)[0]
dsswcs = WCS(dsshdu.header)
dssdata = dsshdu.data

dssxlen, dssylen = dssdata.shape
dsscx = int(0.5 * dssxlen)
dsscy = int(0.5 * dssylen)
dssdx = int(dssxlen*0.1)
dsswcscut = dsswcs[dsscx-dssdx:dsscx+dssdx,dsscy-dssdx:dsscy+dssdx]

plt.figtext(0.45,0.93, name, size="large")
plt.figtext(0.15,0.9, "X-ray, Rosat HRI")
plt.figtext(0.57,0.9, "Optical, UK Shimidt")				

plt.subplot(221, projection=dsswcs)
#		plt.subplot(221, projection=hriwcs)		
try:
	plt.imshow(hridata, origin='lower', norm=LogNorm())
	plt.colorbar()					
except:
	print("ERROR, couldn't plot log-z scale")			
	plt.close()

plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')
plt.ylabel('Galactic Latitude')

plt.subplot(223, projection=dsswcscut)
#		plt.subplot(223, projection=hriwcscut)		
plt.imshow(hridata[hricx-hridx:hricx+hridx,hricy-hridx:hricy+hridx], origin='lower')
plt.colorbar()		
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')
plt.ylabel('Galactic Latitude')


plt.subplot(222, projection=dsswcs)
try:
	plt.imshow(dssdata, origin='lower', norm=LogNorm())
	plt.colorbar()								
except: 
	print("ERROR, couldn't plot log-z scale")
	plt.imshow(dssdata, origin='lower')
	plt.colorbar()		
	
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')

plt.subplot(224, projection=dsswcscut)
plt.imshow(dssdata[dsscx-dssdx:dsscx+dssdx,dsscy-dssdx:dsscy+dssdx], origin='lower')
plt.colorbar()		
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')

plt.savefig(name  + ".png ")
plt.close()

Output result

The top two overviews show the z-axis in log, and the bottom two show in linear. When displaying the z-axis in log, the data can only be empty or zero, so try except is used to avoid it.

Reference example

[How to search using python's astroquery and get fits images with skyview]( An example of a gif animation created by acquiring X-ray (ROSAT) and visible light (DSS) images using https://qiita.com/yamadasuzaku/items/df185d151248578e8e1e). To make a gif animation, I generated it with convert -delay 100 -loop 1 * .png xrayopt.gif.

--X-ray and visible light (~ 200 celestial bodies): http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/xrayopt.gif --Visible light (~ 400 celestial bodies): http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/opt.gif

How to plot multiple fits images side by side in galactic coordinates using python