Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
28
Help us understand the problem. What is going on with this article?
@yamadasuzaku

宇宙の研究を始める人向けのgoogle Colabで行うイメージ解析

はじめに

fits ファイルを扱うのが初めての人向けです.銀河の画像は親しみやすいので,これを例に google Colab だけでも動くサンプルです.もう少し上級者向けに書いた pythonのastroqueryを用いて検索し、skyviewでfits画像を取得する方法 もあります.

基本編

サーバーとクライエントの基本

宇宙天文では,大量のデータをやりとりするので,サーバーとクライエント,について知識が必要になります.クライアントサーバ (client server system) とは,宇宙天文でも基本となるコンセプトなので,知っておこう.簡単には,サーバーとは,ネットワーク上のたくさんの機器に情報を提供できる機能が備わったOSやソフトウェア群のことで,クライエントはサーバーからデータを取得するだめのマシンのことです.

OSによってもかなり違いがあります.宇宙天文では,linuxベースのサーバがほぼ100%だと思いますので,基本的なサーバーとのやりとりや仕組みを事前に知っておくと,宇宙天文のデータを取得する際に少しは気が楽になると思います.サーバーOSとは? などを参考ください.

urllib の例

サーバーとクライエントの仕組みを使う場合は,今ではとても簡単なので,まずは気楽に使ってみよう.例えば,NASAのheasoftという有名なページのファイルを取得してみよう.

import urllib.request
from bs4 import BeautifulSoup
url = 'https://heasarc.gsfc.nasa.gov/docs/software/heasoft/'
req = urllib.request.Request(url)
html = urllib.request.urlopen(req)
soup = BeautifulSoup(html, "html.parser")
topicsindex = soup.find('div', attrs={'class': 'navigation'})

これで,soup に全htmlの情報が取得されている.ブラウザはこれを表示しているのである.この中から,navigation という部分だけ取得したのが,topicsindex であり,

print(topicsindex) # extract sentences begining with div navigation 
<div class="navigation">
<div id="nasaHeader">
<div>
....

のように,一部の情報だけ表示される.これは,クライエントからNASAのサーバーにアクセスし,データをもらうという作業で,数行のプログラムで書けることに注目してほしい.データを取得するというのは,ブラウザからマウスでボタンを色々と押してデータを取得する,というだけではないことを認識して欲しいのと,宇宙天文の場合はこのようなコマンドラインでデータを取得することが多い.

astroquery の設定

pip install astroquery 

astroquery をインストールしよう.これはskyviewやvizierなど,多々のカタログから,同じコーディングでデータを取得するためのツールです.

vizier の使い方

vizierは論文ごとにカタログを抽出したものである.例えば,blackhole で検索すると,大量のカタログを含む論文のリストが現れる.vizier上でも簡単なプロットなどがすぐにできる.論文とも紐付けられているので,もと情報にもすぐに辿れるのも便利.astroquery からは,

# download galaxy data via Vizier
from astroquery.vizier import Vizier 
v = Vizier(catalog="J/ApJS/80/531/gxfluxes",columns=['Name',"logLx","Bmag","dist","type"],row_limit=-1)
data = v.query_constraints()
sname = data[0]["Name"] 
namelist = []
olist=["HRI","DSS"]

J/ApJS/80/531/gxfluxes は銀河のカタログで,ここから,'Name',"logLx","Bmag","dist","type" のデータを抽出して,sname = data[0]["Name"] として, Name だけ取り出している.

skyview の使い方

skyview は,天体の画像を取得するのに便利なカタログです.今回はここから銀河の画像を取得してみましょう.

# get image from skyview
from astroquery.skyview import SkyView
name=namelist[0]
paths = SkyView.get_images(position=name, survey=olist)
save(paths,name,olist)

ここでは,position="天体名", survey=["DSS",..]などを指定して,
skyview から画像を取得し,save という関数で fits ファイルに保存する.

fitsファイルへの保存方法

def save(p,name,obs): # save file into fits 
    for onep,oneo in zip(p,obs):
        onep.writeto(name+"_"+oneo+".fits",overwrite=True)

obs の中に書かれたカテゴリリストの名前をつけて,fits ファイルに保存する.

matplotlib を用いた天空座標でのプロット例

画像をプロットするのに,必要なモジュールをインポートする.

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from astropy.wcs import WCS
from astropy.io import fits
import glob 
import numpy as np

次に,fits ファイルを astropy で読んで,header から WCS 情報を引き出す.
WCS は,World Coordinate System (WCS) の意味で,天体画像の座標変換に必要な情報で,具体的には,検出器座標から天空座標や,銀経銀緯座標への変換に必要な情報です.座標変換は,matplotlib の projection というキーワードに指定するだけで変換してくれます.

# set IC342 imaga of DSS
dssname = "IC0342_DSS.fits"
dsshdu = fits.open(dssname)[0]
dsswcs = WCS(dsshdu.header)
dssdata = dsshdu.data

matplotlib でプロットする.

# plot image 
F = plt.figure(figsize=(8,8))
plt.subplot(111, projection=dsswcs)
plt.imshow(dssdata, origin='lower', norm=LogNorm())
plt.colorbar() 

スクリーンショット 2020-05-15 13.10.22.png

このような絵が表示されれば成功です.

画像解析の例

画像の回転,転置,差分,など基本的な操作をして画像で遊んでみよう.
例えば,差分画像も簡単に表示できます.

スクリーンショット 2020-05-15 13.04.37.png

fits ファイルとは,テキストで書かれた基礎情報と,バイナリで書かれたデータが合わさったファイル,というだけで,本当にそれだけなので,まずは苦手意識を持たずに扱おう.画像データは,openCV など,python にはたくさんのモジュールがあるので,まずはいろいろと遊んでみましょう.画像も,銀河でも波長が違えば見え方が違うし,超新星残骸,銀河団など,他の画像で遊んでみるのもよいです.画像の検索方法は,vizier でもよいですが,意外と欲しいものが見つからないこともあるので,適当にgoogleで検索した方が早いときもあります..

モルフォロジー変換を用いた画像解析の例

専門的な天体解析では,望遠鏡の応答関数(像の広がり方,感度など)を用いて,画像中にある"点"が点源なのか,広がった構造なのかを調べます.例えば,Chandra衛星には,点源の解析のためのパターンマッチングを行う wavdetect というツールが用意されています.

ここではもう少し簡単に天体画像を遊べる方法を紹介します.モルフォロジー変換というのもので,画像を2値化(0 or 1)して,適当な小さい領域との and/or 演算を繰り返して画像の形状変換を行う方法です.例えば,銀河の面積が100ピクセルで繋がった構造をしていて,点源は一ピクセルだけだとすると,収縮(Erosion)処理(=適当な小さい領域に1ピクセルでも0があれば0とする演算)を行うと,ノイズは消えます.その次に,膨張(Dilation)処理(==適当な小さい領域に1ピクセルでも1があれば1とする演算)を行うと元に戻りますが,一旦0になったノイズは出てきません.このように,収縮(Erosion)と膨張(Dilation)を繰り返して,ノイズを処理する方法を,オープニング処理と呼びます.cv2.morphologyEx(binarydssdata,cv2.MORPH_OPEN, kernel) を用いて実行できます.モルフォロジー変換が参考になりました.大元はC言語で書かれていて,Morphology_2_8cpp-example の実行例があります.

使い方は簡単で,画像を,

# convert the data into binary image
binarydssdata = 1.0 * (dssdata > np.median(dssdata)) 

のように,median よりも大小で2値化します.int だと動かない関数があるようなので,float の 1.0 にします.

次に,getStructuringElement(cv2.MORPH_ELLIPSE,(8,8)) で,8x8の画像で,楕円形の 1, 0 の行列を生成して,これを kernel として,cv2.morphologyEx として変換します.

# 1.0 is used because morphologyEx may not work for int32
import cv2
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(8,8))
opening = cv2.morphologyEx(binarydssdata,cv2.MORPH_OPEN, kernel)

コードはこれだけです.実行例はこちらです.

スクリーンショット 2020-05-19 13.43.35.png

左から順に,「元画像」「2値化画像」「モルフォロジー変換後」で,一番右は,「元画像」 x 「モルフォロジー変換後」をしたものです.実際の天体解析でも,点源除去は同じような操作をします.

google Colab での実行例

google Colab 上で実行するには,astroquery がデフォルトではインストールされてないので,astroquery の install が必要です.google Colab 上での実行例 を参照ください.

また,fits ファイルを google drive 上に保存しないと,データは外に出せないので,必要ならgoogle drive をマウントして,そこに fits ファイルを書き出しましょう.

関連ページ

28
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
yamadasuzaku
宇宙物理、X線観測、超伝導検出器など.
RikkyoU
立教大学のAI/IT愛好家による教育研究コンテンツの発信.

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
28
Help us understand the problem. What is going on with this article?