WCS座標でイメージをマージしpythonでプロットする方法

# 背景

WCS(ra, dec)などの天球座標のイメージの四則演算や結合から、ds9 ではなく python でイメージの解析を行いたい場合にどうするのがよいかまとめる(2018/08/21時点)。動作は、OSは mac のバージョン 10.11 でのみ検証、原理的には linux でも動くはずである。

# 簡単なまとめ

## WCS座標でのイメージの足しあわせ方

1. IRAF 
IRAF の [imcombine](http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?imcombine) を使うのが一番、安定性とメンテの良さの点ではよい。mac にも簡単にインストールできる(下記に詳細)。

2. swarp (http://www.astromatic.net/software/swarp)
ピクセルサイズ程度の解析をしたい場合は、swarp で足したほうがよい。普通の ximage や iraf のsum は光子数 or イベント数が保存する足し方(おそらく)で、その場合、軸の異なる座標系に射影する際に、元来の輝度中心がずれる。swarp はその辺りも考慮して、resampling と interpolation も実装していて、フラックスの絶対値だけでなくて、空間分布も保存するように足すこともできる。

3. その他
ximage の sum 関数 (https://heasarc.gsfc.nasa.gov/ftools/fhelp/ximage.sum )のユーザーが多い。ただし情報が細い。


## python でのイメージの処理

1. astropy (http://www.astropy.org)
matplotlib ユーザーには一番わかりやすい。ただし、電波のコントアとX線に重ねる場合に計算がものすごい遅い。ax.get_transform が遅いと思われる。自分自身のコントアであれば座標変換が不要なので計算が早い。

2. aplpy (https://aplpy.github.io)
一気に論文クオリティの絵ができるのが売り。ちょっと癖があるが、pythonに慣れてれば問題ないと思う。コントアを重ねるのもものすごく早いので、この点でastropyではなくaplpyを使うメリットは高い。

3. ３色図の作り方
montage_wrapper を APLpy (http://aplpy.github.io) から使うと、
RGBの図であれば、裏で Montage (https://github.com/Caltech-IPAC/Montage) が走って繋げてくれる。http://python4esac.github.io/plotting/aplpy.html の "Color Images" の箇所を参考。ds9 の -rgb オプション (http://ds9.si.edu/doc/user/rgb/index.html )でも可能なので、あまり使う機会はないと思われる。


# 具体例

## IRAF のインストール方法 on mac 

http://iraf.noao.edu
から、Mac OSX 64-bit (183MB)
ftp://iraf.noao.edu/iraf/v216/PCIX/iraf.macx.x86_64.tar.gz
をダウンロードする。

``` bash:
# IRAF 2.16
iraf=/iraf/iraf
```

にダウンロードしたファイルを置いた場合、

``` bash:

[syamada] $ cd $iraf                                                        [~]
[syamada] $ ls                                                     [/iraf/iraf]
HS.PCIX.GEN    README.install doc            math          util
IRAF.NET       RELEASE.txt    extern        mkpkg          vendor
IS.PORT.GEN    bin            include        noao          vo
LICENSE        bin.generic    install        pkg
Makefile       bin.macintel  lib            sys
README         dev            local          unix
[syamada] $ $iraf/install 

Enter default terminal type (xgterm): xterm

```

で、xterm に変更して、== Installation Completed With No Errors　== が出たらインストールに成功。zshユーザーは手動で、

``` zsh:
# Add iraf setup commands
if [ -e /Users/syamada/.iraf/setup.sh ]; then
    source /Users/syamada/.iraf/setup.sh
fi
```

を.zshrcに加筆すると、

``` zsh:
[syamada] $ cl
setting terminal type to 'xterm' ...
NOAO/IRAF PC-IRAF Revision 2.16.1 EXPORT Mon Oct 14 21:40:13 MST 2013
This is the EXPORT version of IRAF V2.16 supporting PC systems.
  The following commands or packages are currently defined:
      dataio.     language.   obsolete.   softools.   vo.         
      dbms.       lists.      plot.       system.     
      images.     noao.       proto.      utilities.  
vocl> 
```
と出力されたらOK。最近は、ecl の方が新しいので、ecl を使うのが標準らしい。

## WCS上で 座標の異なる XMM Newton のイメージを足して、exposure map で割る

``` zsh:

ls xmm_750_1300_eV_*.fits > xmm1p3keV.list # 輝度分布画像fitsのファイル名リスト
ls exp_750_1300_eV_*.fits > exp1p3keV.list # exposure画像fitsのファイル名リスト

ecl <<EOF

imcombine @xmm1p3keV.list xmm1p3keV.fits combine="sum" offsets="wcs"
imcombine @exp1p3keV.list exp1p3keV.fits combine="sum" offsets="wcs"
imarith xmm1p3keV.fits / exp1p3keV.fits xmm1p3keV_cor.fits

EOF

```

このページ(https://qiita.com/yamadasuzaku/items/0c5f28438f7342b63975) のスクリプトを使った場合は、adapt-*fits が adaptively smooth された画像で、それに対応する exposure map が size-*fits である。その場合は、
http://www-x.phys.se.tmu.ac.jp/~syamada/syamadatools/xmm-newton/link_image.sh
を参考。

## astropy の使い方

特に癖はないので、image を読みこんで、matplotlib で通常の使い方をすればよい。

``` python:
        ax = plt.subplot(111,projection=self.wcs)
        ax.grid(color='black', ls='dotted')
```

のように、projection で wcs を指定すると自動で wcs で表示してくれる。grid は、ax の属性から呼ばないとだめなので、ax.grid を使う。

ある画像に、その画像のコントアを重ねたいのであれば、ax.contour を使うだけ。もし、座標系が異なる場合は、ax.contour で、transform=ax.get_transform オプションを指定すればよい。ただし、画像によって時間が著しくかかる。

- 参考 http://www-x.phys.se.tmu.ac.jp/~syamada/syamadatools/xmm-newton/xmm_process_plotimage.py

## aplpy の使い方

[aplpy.FITSFigure](https://aplpy.readthedocs.io/en/stable/api/aplpy.FITSFigure.html) で イメージの fits を開けばあとは用意されている関数を羅列するだけ。figure 属性も呼び出すのが astropy との違い。コントアは show_contour 関数を呼び出すだけでよく、astropy のような計算の遅さは今の所ない。matplotlib の rapper なので、オプションは matplotlib とほど同じものが使える。

- 参考 http://www-x.phys.se.tmu.ac.jp/~syamada/syamadatools/xmm-newton/xmm_process_plotimage_aplpy.py

上記スクリプトで生成された画像。

<img width="407" alt="ss433.png" src="https://qiita-image-store.s3.amazonaws.com/0/179652/76be36cc-3ea8-973d-32d1-0a8419f1d6e1.png">


## swarp のインストール方法と使い方

http://www.astromatic.net/download/swarp/
から swarp-2.16.4.tar をダウンロードし、mac の EL Capitan にインストールできた。

``` bash:
./configure CC=gcc
make 
sudo make install
```

これで、/usr/local/bin/swarp が生成される。

使い方は、マニュアルの example を参考にし、swarp -d で設定を確認し必要なら変更し、swarp *fits で実行するだけ。