Comment créer un profil radial à partir d'images astronomiques (Chandra, XMM etc.) en utilisant python
Contexte
En créant un profil radial d'une image astronomique, une méthode générale d'analyse d'image astronomique, par exemple si la source ponctuelle a une diffusion cohérente ou une diffusion plus large que la diffusion finie d'un télescope ou d'un détecteur. Il peut être inclus dans des outils d'analyse à usage général, mais si vous souhaitez effectuer une analyse détaillée, vous pouvez l'écrire vous-même.
Comment faire
Déterminez simplement le centre, préparez un anneau à partir de celui-ci, comptez le nombre de pixels contenus dans l'anneau et divisez par la zone de l'anneau.
python
def calc_radialprofile(data, xg, yg, ds = 10, ndiv = 10, debug = False):
"""
calc simple peak (just for consistendety check) and baricentric peak.
"""
tmp = data.T
nr = np.linspace(0, ds, ndiv)
rc = (nr[:-1] + nr[1:]) * 0.5
rp = np.zeros(len(nr)-1)
nrp = np.zeros(len(nr)-1)
for i in range(ds*2):
for j in range(ds*2):
itx = int(xg - ds + i)
ity = int(yg - ds + j)
val = tmp[itx][ity]
dist = dist2d(itx, ity, xg, yg)
for k, (rin,rout) in enumerate(zip(nr[:-1], nr[1:])):
if dist >= rin and dist < rout:
rp[k] = rp[k] + val
# normalize rp
for m, (rin,rout) in enumerate(zip(nr[:-1], nr[1:])):
darea = np.pi * ( np.power(rout,2) - np.power(rin,2) )
nrp[m] = rp[m] / darea
if debug:
print(m, rp[m], nrp[m], rin, rout, darea)
return rc, nrp
À proprement parler, si l'élément d'image est un CCD, il peut avoir un GAP, et selon l'emplacement, il peut ne pas s'agir d'un anneau mais d'une partie sans pixels, donc c'est un anneau et efficace. Il est nécessaire de calculer avec une certaine zone sensible. Outil Pileup pour détecteur XIS (CCD) du satellite "Suzaku" [Recette de l'effet Pileup sur le Suzaku XIS](http://www-x.phys.se.tmu.ac.jp/~syamada/ana/suzaku /XISPileupDoc_20120221/XIS_PileupDoc_20190118_ver2.html) se divise par la surface effective en fonction du mode fenêtre 1/4 ou 1/8, du mode de synchronisation, etc.
Exemple de code pour créer un profil radial en python
Aperçu du code
Saisissez le fichier image des ajustements. Pour créer un profil radial d'une "structure ponctuelle" dans l'image
- Coordonnées du centre (x, y) --Taille (ds) --Nombre de divisions de l'anneau (ndiv)
Est le paramètre minimum requis. L'image peut être un nombre de photons ou un taux de comptage, et la plage dynamique est différente, donc si vous spécifiez -m -a VMAX -i VMIN, la plage d'affichage et l'aspect seront dans la plage de VMIN à VMAX. , Si rien n'est spécifié, la proximité du minimum et du maximum des données est affichée de manière appropriée.
python
[syamada] $ python qiita_mk_radialprofile.py -h
Usage: qiita_mk_radialprofile.py FileList
Options:
--version show program's version number and exit
-h, --help show this help message and exit
-d, --debug The flag to show detailed information
-m, --manual The flag to use vmax, vmin
-x DETX, --detx=DETX det x
-y DETY, --dety=DETY det y
-a VMAX, --vmax=VMAX VMAX
-i VMIN, --vmin=VMIN VMIN
-s DS, --dataExtractSize=DS
data size to be extracted from image
-n NDIV, --numberOfDivision=NDIV
number of division of annulus
Comment utiliser
python
python qiita_mk_radialprofile.py ss433_659_broad_flux.img
Ensuite, l'image suivante est créée. En même temps, il vide le profil radial dans le texte et enregistre l'image pour une utilisation de débogage.
Télécharger un échantillon
Exemple de satellite à rayons X
Pour XMM Newton
J'ai essayé d'expliquer l'analyse des données du satellite XMM-Newton en japonais simple et [xmm_process_radialprofile.py](http: //www-x.phys) .se.tmu.ac.jp /% 7Esyamada / syamadatools / xmm-newton / py3 / xmm_process_radialprofile.py) est introduit dans le code source.
Pour Chandra
Il peut également être créé à l'aide d'un groupe d'outils appelé ciao. La méthode suivante est introduite dans le fil créer un profil radial.
python
$ download_chandra_obsid 17395
$ chandra_repro 17395 outdir=./
$ fluximage hrcf17395_repro_evt2.fits coma
$ dmextract "hrcf17395_repro_evt2.fits[bin sky=annulus(16635,15690,0:3200:100)]" coma_prof.rad \
exp=coma_wide_thresh.expmap clob+ op=generic
Pour des informations générales sur le satellite Chandra, reportez-vous à J'ai essayé d'expliquer l'analyse des données du satellite Chandra en japonais facile.
code
qiita_mk_radialprofile.py
import os
import sys
import math
import subprocess
import re
import sys
# for I/O, optparse
import optparse
# numpy
import numpy as np
from numpy import linalg as LA
# conver time into date
import datetime
# matplotlib http://matplotlib.org
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter, FixedLocator
import matplotlib.colors as cr
from matplotlib.font_manager import fontManager, FontProperties
import matplotlib.pylab as plab
## for FITS IMAGE
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import get_pkg_data_filename
params = {'xtick.labelsize': 10, # x ticks
'ytick.labelsize': 10, # y ticks
'legend.fontsize': 8
}
plt.rcParams['font.family'] = 'serif'
plt.rcParams.update(params)
# Fits I.O class
class Fits():
# member variables
debug = False
def __init__ (self, eventfile, debug):
self.eventfilename = eventfile
self.debug = debug
#extract header information
if os.path.isfile(eventfile):
self.filename = get_pkg_data_filename(eventfile)
self.hdu = fits.open(self.filename)[0]
# self.dateobs = self.hdu.header['DATE-OBS']
# self.object = self.hdu.header['OBJECT']
self.wcs = WCS(self.hdu.header)
self.cdelt = self.wcs.wcs.cdelt
self.p2deg = np.abs(self.cdelt[0]) # [deg/pixel]
self.p2arcsec = 3600. * np.abs(self.cdelt[0]) # [arcsec/pixel]
self.dirname = eventfile.replace(".fits","").replace(".gz","")
print("..... __init__ ")
# print " self.dateobs = ", self.dateobs
# print " self.object = ", self.object
print(" self.p2deg = ", self.p2deg, " [deg/pixel]")
print(" self.p2arcsec = ", self.p2arcsec, " [arcsec/pix]")
self.data = self.hdu.data
self.data[np.isnan(self.data)] = 0 # nan ==> 0, note that radio data include negative values
self.data = np.abs(self.data) # data is forced to be positive (Is it right?)
if self.debug:
print("... in __init__ Class Fits ()")
print("filename = ", self.filename)
print("hdu = ", self.hdu)
print("data.shape = ",self.data.shape)
self.wcs.printwcs()
else:
print("ERROR: cat't find the fits file : ", self.eventfilename)
quit()
def plotwcs(self, vmin = 1e-1, vmax = 20, manual = False):
"""
just plot the entire image
"""
print("\n..... plotwcs ...... ")
if not manual:
vmin = np.amin(self.hdu.data) + 1e-10
vmax = np.amax(self.hdu.data)
# plt.figtext(0.1,0.97, "OBJECT = " + self.object + " DATE = " + self.dateobs)
plt.figtext(0.1,0.95, "Unit = " + str(self.p2arcsec) + " [arcsec/pix]" )
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
plt.colorbar()
plt.xlabel('RA')
plt.ylabel('Dec')
outputfiguredir = "fig_image_" + self.dirname
subprocess.getoutput('mkdir -p ' + outputfiguredir)
plt.savefig(outputfiguredir + "/" + "plotwcs.png ")
def plotwcs_ps(self, detx, dety, ds = 40, vmin = 1e-1, vmax = 20, manual = False):
"""
just plot the enlarged image around (detx,dety)
"""
print("\n..... plotwcs_ps ...... ")
if not manual:
vmin = np.amin(self.hdu.data) + 1e-10
vmax = np.amax(self.hdu.data)
self.detx = detx
self.dety = dety
gpixcrd = np.array([[ self.detx, self.dety]], np.float_)
gwrdcrd = self.wcs.all_pix2world(gpixcrd,1)
ra = gwrdcrd[0][0]
dec = gwrdcrd[0][1]
self.ra = ra
self.dec = dec
print("detx, dety = ", detx, dety)
print("ra, dec = ", ra, dec)
F = plt.figure(figsize=(12,8))
# plt.figtext(0.1,0.97, "OBJECT = " + self.object + " DATE = " + self.dateobs)
plt.figtext(0.1,0.95, "Unit = " + str(self.p2arcsec) + " [arcsec/pix]" )
ax = plt.subplot(111, projection=self.wcs)
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
print(" self.hdu.data[detx][dety] = ", self.hdu.data[detx][dety] , " (detx,dety) = (", detx, ",", dety, ")")
ax.set_xlim(detx - ds, detx + ds)
ax.set_ylim(dety - ds, dety + ds)
plt.colorbar()
plt.xlabel('RA')
plt.ylabel('Dec')
outputfiguredir = "fig_image_ps_" + self.dirname
subprocess.getoutput('mkdir -p ' + outputfiguredir)
plt.savefig(outputfiguredir + "/" + "plotwcs_ps_detx" + str("%d" % detx) + "_dety" + str("%d" % dety) + ".png ")
def mk_radialprofile(self, detx, dety, ds = 40, ndiv = 20, vmin = 1e-1, vmax = 20, manual = False):
"""
create the radial profiles
"""
print("\n..... mk_radialprofile ...... ")
if not manual:
vmin = np.amin(self.hdu.data) + 1e-10
vmax = np.amax(self.hdu.data)
self.detx = detx
self.dety = dety
gpixcrd = np.array([[ self.detx, self.dety]], np.float_)
gwrdcrd = self.wcs.all_pix2world(gpixcrd,1)
ra = gwrdcrd[0][0]
dec = gwrdcrd[0][1]
self.ra = ra
self.dec = dec
print("detx, dety = ", detx, dety)
print("ra, dec = ", ra, dec)
# radial profiles around the input (ra, dec)
rc, rp = calc_radialprofile(self.data, detx, dety, ds = ds, ndiv = ndiv)
# plot images and radial profiles
F = plt.figure(figsize=(10,12))
F.subplots_adjust(left=0.1, bottom = 0.1, right = 0.9, top = 0.87, wspace = 0.3, hspace=0.3)
# plt.figtext(0.1,0.97, "OBJECT = " + self.object + " DATE = " + self.dateobs)
plt.figtext(0.1,0.95, "Unit = " + str(self.p2arcsec) + " [arcsec/pix]" )
plt.figtext(0.1,0.93, "input center [deg](x,c,ra,dec) = (" + str("%3.4f" % detx) + ", " + str("%3.4f" % dety) + ", "+ str("%3.4f" % ra) + ", " + str("%3.4f" % dec) + ")")
ax = plt.subplot(3,2,1)
plt.title("(1) SKY image")
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
# ax.set_xlim(detx - ds, detx + ds)
# ax.set_ylim(dety - ds, dety + ds)
# plt.imshow(self.hdu.data, origin='lower', cmap=plt.cm.viridis)
plt.colorbar()
plt.scatter(detx, dety, c="k", s=300, marker="x")
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.subplot(3,2,2, projection=self.wcs)
plt.title("(2) Ra Dec image (FK5)")
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
plt.colorbar()
plt.scatter(detx, dety, c="k", s=300, marker="x")
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.subplot(3,2,3)
plt.title("(3) SKY image")
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
ax.set_xlim(detx - ds, detx + ds)
ax.set_ylim(dety - ds, dety + ds)
plt.colorbar()
plt.scatter(detx, dety, c="k", s=300, marker="x")
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.subplot(3,2,4, projection=self.wcs)
plt.title("(4) Ra Dec image (FK5)")
plt.imshow(self.data, origin='lower', cmap=plt.cm.jet, norm=cr.LogNorm(vmin=vmin,vmax=vmax))
ax.set_xlim(detx - ds, detx + ds)
ax.set_ylim(dety - ds, dety + ds)
plt.colorbar()
plt.scatter(detx, dety, c="k", s=300, marker="x")
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.subplot(3,2,5)
plt.title("(5) radial profile (pix)")
plt.errorbar(rc, rp, fmt="ko-", label="input center")
plt.xlabel('radial distance (pixel)')
plt.ylabel(r'c s$^{-1}$ deg$^{-2}$')
plt.grid(True)
plt.legend(numpoints=1, frameon=False)
ax = plt.subplot(3,2,6)
plt.title("(6) radial profile (arcsec)")
plt.errorbar(rc * self.p2arcsec, rp, fmt="ko-", label="input center")
plt.xlabel('radial distance (arcsec)')
plt.ylabel(r'c s$^{-1}$ deg$^{-2}$')
plt.legend(numpoints=1, frameon=False)
plt.grid(True)
outputfiguredir = "fig_image_xy_" + self.dirname
subprocess.getoutput('mkdir -p ' + outputfiguredir)
plt.savefig(outputfiguredir + "/" + "mkrad_detx" + str("%d" % detx) + "_dety" + str("%d" % dety) + ".png ")
# dump the radial profiles into the text file
outputfiguredir = "txt_image_xy_" + self.dirname
subprocess.getoutput('mkdir -p ' + outputfiguredir)
fout = open(outputfiguredir + "/" + "mkrad_detx" + str("%d" % detx) + "_dety" + str("%d" % dety) + ".txt", "w")
for onex, oney1 in zip(rc * self.p2arcsec, rp):
outstr=str(onex) + " " + str(oney1) + " \n"
fout.write(outstr)
fout.close()
def calc_radialprofile(data, xg, yg, ds = 10, ndiv = 10, debug = False):
"""
calc simple peak (just for consistendety check) and baricentric peak.
"""
tmp = data.T
nr = np.linspace(0, ds, ndiv)
rc = (nr[:-1] + nr[1:]) * 0.5
rp = np.zeros(len(nr)-1)
nrp = np.zeros(len(nr)-1)
for i in range(ds*2):
for j in range(ds*2):
itx = int(xg - ds + i)
ity = int(yg - ds + j)
val = tmp[itx][ity]
dist = dist2d(itx, ity, xg, yg)
for k, (rin,rout) in enumerate(zip(nr[:-1], nr[1:])):
if dist >= rin and dist < rout:
rp[k] = rp[k] + val
# normalize rp
for m, (rin,rout) in enumerate(zip(nr[:-1], nr[1:])):
darea = np.pi * ( np.power(rout,2) - np.power(rin,2) )
nrp[m] = rp[m] / darea
if debug:
print(m, rp[m], nrp[m], rin, rout, darea)
return rc, nrp
def dist2d(x, y, xg, yg):
return np.sqrt( np.power( x - xg ,2) + np.power( y - yg ,2) )
def calc_center(data, idetx, idety, ds = 10):
"""
calc simple peak (just for consistendety check) and baricentric peak.
"""
tmp = data.T
xmax = -1.
ymax = -1.
zmax = -1.
xg = 0.
yg = 0.
tc = 0.
for i in range(ds*2):
for j in range(ds*2):
tx = idetx - ds + i
ty = idety - ds + j
val = tmp[tx][ty]
tc = tc + val
xg = xg + val * tx
yg = yg + val * ty
if val > zmax:
zmax = val
xmax = idetx - ds + i
ymax = idety - ds + i
if tc > 0:
xg = xg / tc
yg = yg / tc
else:
print("[ERROR] in calc_center : tc is negaive. Something wrong.")
sys.exit()
print("..... in calc_center : [simple peak] xmax, ymax, zmax = ", xmax, ymax, zmax)
print("..... in calc_center : [total counts in ds] tc = ", tc)
if zmax < 0:
print("[ERROR] in calc_center : zmax is not found. Something wrong.")
sys.exit()
return xg, yg, tc
def main():
""" start main loop """
curdir = os.getcwd()
""" Setting for options """
usage = '%prog FileList '
version = __version__
parser = optparse.OptionParser(usage=usage, version=version)
parser.add_option('-d', '--debug', action='store_true', help='The flag to show detailed information', metavar='DEBUG', default=False)
parser.add_option('-m', '--manual', action='store_true', help='The flag to use vmax, vmin', metavar='MANUAL', default=False)
parser.add_option('-x', '--detx', action='store', type='int', help='det x', metavar='DETX', default=61)
parser.add_option('-y', '--dety', action='store', type='int', help='det y', metavar='DETY', default=45)
parser.add_option('-a', '--vmax', action='store', type='float', help='VMAX', metavar='VMAX', default=4e-6)
parser.add_option('-i', '--vmin', action='store', type='float', help='VMIN', metavar='VMIN', default=1e-10)
parser.add_option('-s', '--dataExtractSize', action='store', type='int', help='data size to be extracted from image', metavar='DS', default=40)
parser.add_option('-n', '--numberOfDivision', action='store', type='int', help='number of division of annulus', metavar='NDIV', default=20)
options, args = parser.parse_args()
print("---------------------------------------------")
print("| START : " + __file__)
debug = options.debug
manual = options.manual
detx = options.detx
dety = options.dety
vmax = options.vmax
vmin = options.vmin
ds = options.dataExtractSize
ndiv = options.numberOfDivision
print("..... detx ", detx)
print("..... dety ", dety)
print("..... ds ", ds, " (bin of the input image)")
print("..... ndiv ", ndiv)
print("..... debug ", debug)
print("..... manual ", manual)
print("..... vmax ", vmax)
print("..... vmin ", vmin)
argc = len(args)
if (argc < 1):
print('| STATUS : ERROR ')
print(parser.print_help())
quit()
filename=args[0]
print("\n| Read each file and do process " + "\n")
print("START : ", filename)
eve = Fits(filename, debug)
eve.plotwcs(vmax = vmax, vmin = vmin, manual = manual) # plot the entire image
eve.plotwcs_ps(detx, dety, ds = ds, vmax = vmax, vmin = vmin, manual = manual) # plot the enlarged image around (detx,dety).
eve.mk_radialprofile(detx, dety, ds = ds, ndiv = ndiv, vmax = vmax, vmin = vmin, manual = manual) # create detailed plots and radial profiles
print("..... finish \n")
if __name__ == '__main__':
main()
Recommended Posts