Comment superposer des images en utilisant bien Python dans l'espace WCS

Environnement

OK si matplotlib et astropy sont installés en python. Si le système anconda3 est inclus, il est OK par défaut.

Un exemple de chevauchement de deux images avec une trame d'argent

En prenant comme exemples les images de lumière visible et de rayons X NGC4051, le tracé d'image entier et le tracé agrandi près du centre sont montrés à titre d'exemples.

point

Il suffit de générer une classe WCS à partir du fichier ajustements d'image et de la placer dans l'option de projection de matplotlib. Cependant, avec cela seul, il n'est affiché que dans la taille d'image initiale.

--Si vous souhaitez agrandir un certain endroit, coupez l'objet WCS de la même manière que le fichier image. La partie hriwcs [hricx-hridx: hricx + hridx, hricy-hridx: hricy + hridx] de l'exemple de code. Vous pouvez donner ceci à l'option wcs de matplotlib.

――Si vous souhaitez afficher deux images l'une sur l'autre, cela devient un problème difficile. Dans ce cas, il est brièvement résumé dans Comment fusionner des images avec des coordonnées WCS et tracer avec python.

Stockage de code et de fichiers

qiita_plot_fits_2image.py -Image visible et radiographique de NGC4051
- 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

Pour savoir comment prendre une image, Comment rechercher à l'aide de la requête astro de python et obtenir une image ajustée avec skyviewを参考にされたい。

Exemple de 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()

Résultat de sortie

Les deux aperçus du haut montrent l'axe z sous forme de journal et les deux derniers sous forme linéaire. Lors de l'affichage de l'axe z sous forme de journal, on peut dire qu'il n'y a que des données vides ou nulles, alors essayez sauf est utilisé pour l'éviter.

Exemple de référence

[Comment rechercher à l'aide de la requête astro de python et obtenir des images ajustées avec skyview]( Un exemple d'animation gif créée en acquérant des images de rayons X (ROSAT) et de lumière visible (DSS) à l'aide de https://qiita.com/yamadasuzaku/items/df185d151248578e8e1e). Pour créer une animation gif, je l'ai générée avec convert -delay 100 -loop 1 * .png xrayopt.gif.

Rayons X et lumière visible (~ 200 corps célestes): http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/xrayopt.gif
Lumière visible (~ 400 corps céleste): http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/opt.gif

Comment tracer plusieurs images ajustées côte à côte en utilisant Python