OK if matplotlib and astropy are installed in python. If anconda3 system is included, it is OK by default.
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.
--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.
For how to take an image, How to search using python's astroquery and get fits image with skyviewを参考にされたい。
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()
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.
[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
Recommended Posts