python を用いて天体画像(Chandra, XMM etc.)から radial profile を作成する方法

背景

天体画像の radial profile を作ることで，望遠鏡や検出器の有限の広がりと比較し，点源のコンシステントな広がりなのか，それよりも広がっているのかなど，天体の画像解析の一般的な方法であり，汎用的な解析ツールに含まれることもあるが，細かい解析をする場合は，自分で書いてみたい時もあるだろう．

作り方

実直に，中心を決めて，そこから円環を用意して，円環の中に含まれるピクセルのカウント数を数えて円環の面積で割るだけである．

def calc_radialprofile(data, xg, yg, ds = 10, ndiv = 10, debug = False):
    """
    calc simple peak (just for check) and barycentric 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

一応，図式で説明すると，このように円環を考えて，そこに入るカウント数を数えて，円環の面積で割るだけである．

これだけであるが，厳密には，イメージ素子がCCDの場合は，GAPがあることや，場所によっては，円環ではなくて，一部が画素のない部分のこともあるので，円環かつ有効な感度がある面積で計算する必要がある．「すざく」衛星のXIS(CCD)検出器のパイルアップツール Recipe for Pileup Effect on the Suzaku XIS では，1/4 or 1/8 Window mode, timing mode, などに応じて，実効的な面積で割るようにしている．

python で radial profile を作るサンプルコード

コード概要

fits のイメージファイルを入力とする．イメージの中のある"点構造"の radial profile を作るためには，

• 中心座標(x,y)
• サイズ(ds)
• 円環の分割数(ndiv)

が最低限必要なパラメータとなる．イメージは光子のカウントの時もあれば，カウントレートの時もあり，ダイナミックレンジが違うので，表示の範囲や見栄えは -m -a VMAX -i VMIN を指定すると，VMIN から VMAX の範囲になるが，何も指定しなければ，適当にデータの最小最大付近を表示する．

[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

使い方

python qiita_mk_radialprofile.py ss433_659_broad_flux.img

とすると，次のような絵が作成される．同時に radial profile をテキストにダンプしたり，debug 用途のイメージの保存してくれる．

サンプルのダウンロード

• qiita_mk_radialprofile.py or qiita_mk_radialprofile.py
• Chandra の SS433 のイメージ or Chandra の SS433 のイメージ

X線衛星の例

XMM Newton の場合

XMM-Newton 衛星のデータ解析について易しい日本語で解説してみた に，xmm_process_radialprofile.py というソースコードで紹介している．

Chandra の場合

ciao というツール群を使っても作成可能．make a radial profile というスレッドに下記の方法が紹介されている．

$ 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

Chandra衛星の全般的なことは，Chandra 衛星のデータ解析について易しい日本語で解説してみたを参照されたい．

コード

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