파이썬을 이용하여 허셜 우주망원경의 적외선 관측 데이터 분석하기

오늘 파이썬 포스팅은 그동안 써 왔던 내용 중 가장 내용이 어렵지 않을까 싶습니다. 

 

그럼에도 궂이 이 내용을 올리는 것은 실제 학생들과 해 보니 꾀 신기해 하며 열심히 따라오려는 모습에 생각보다 관심이 많구나 라는 생각이 들어서입니다. 이번 포스팅에서는 파이썬을 이용해 허셜 우주망원경이 관측한 데이터를 시각화 하는 과정에 대해 알아보겠습니다.

 

  1. 허셜우주망원경

미국에 나사(NASA)가 있다면 유럽에는 유럽항공우주국(ESA)이 있습니다. 나사에서 허블 우주망원경을 운영했다면, 유럽에서는 허셜 우주망원경을 운영했는데, 허블 망원경과 달리 주로 적외선 파장대역에서 관측을 수행하였습니다.

<허셜 우주망원경>

운영기간이 2009~2013년으로 벌써 10년이 지나 지금은 운영하지 않는 것으로 알고 있는데, 허셜 망원경이 관측한 데이터는 지금도 연구에 활용하고 있습니다. 

 

제가 알기로 허셜 망원경이 관측한 데이터는 주로 성간물질, 항성 탄생 초기과정 등을 연구하는데 활용하는 것으로 알고 있습니다. 실제로 별 탄생 초기의 성간물질의 구조 등을 밝혀내는데 큰 역할을 하였다고 합니다.

 

허셜 망원경이 관측한 데이터는 누구나 얼마든지 다운받을 수 있습니다. 아래 링크에서 다운받을 수 있으니 참고하세요.

** ESA Herschel Science Archive

ESA Herschel Science Archive

여기서는 파일을 다운로드 하는 과정은 따로 언급하지 않고, 이는 다른 글에서 언급하겠습니다. 내용이 다소 길어질것 같아 다음 포스팅 때 이에 대해 상세히 서술하겠습니다. 그래서, 실습용 FIS 파일을 따로 업로드 해 두겠습니다.

  2. FITS라는 파일 형식

대부분의 분들에게 다소 생소한 파일 형식일 것이라 생각합니다. 보통은 거의 사용하지 않는 파일 형식인데, 다른 분야는 잘 모르겠지만 천문학자들은 굉장히 자주 사용하는 파일 형식이라고 합니다. 천문학자들은 방식의 차이는 있지만 사진을 촬영하여 연구합니다. 그럼 이미지 센서는 하나하나의 픽셀에 촬영한 대상의 빛의 세기를 저장하는데, FITS파일은 이런 데이터를 처리하는데 매우 유리힙니다. 또한, 하나의 파일에 여러 유형의 관측 결과를 모두 저장할 수 있는 장점이 있습니다. 여기서 이것까지 모두 설명하기에는 글의 논지를 한참 벗어나니 일단 그런게 있다 정도만 논의하겠습니다. 간단하게라도 언급하는 이유는, 여기서 다룰 파일 형식도 FITS이기 때문입니다.

  3. 자료 다운로드

자료를 다운로드 하는 과정은 여기서는 생략하고, 제가 미리 다운로드한 자료를 사용토록 하겠습니다. 여기서는 말머리성운 IC434를 관측한 결과를 이용하겠습니다. 이 아이를 사용하는 이유는 간단합니다. 대중적으로 많이 알려져있는 대상이기 때문입니다.

IC434.fits

1.82MB

위에 첨부한 파일을 다운로드 하신뒤 적당한 폴더에 두고 실습에 사용하시면 됩니다. 허셜 아카이브에서 자료를 다운로드하는 과정은 아래 포스팅을 참고해 주세요

허셜 우주망원경 관측데이터 다운로드하기

반응형

  4. 파이썬 코드

저 데이터를 분석하는데 매우 필수적으로 사용되는 라이브러리는 astropy라는 아이입니다. 일반적으로는 거의 사용되지 않으나 천문학자들은 굉장히 많이 사용합니다. 아나콘다를 설치할 때 보통 함께 설치되기는 하여 따로 설치할 필요가 없지만, 만약 없다면 따로 설치하셔야겠죠?

4-1. 필수 라이브러리 추출

import matplotlib.pyplot as plt ## 컬러맵을 그리기 위해
from astropy.wcs import WCS  ##관측 좌표계 불러오기
from astropy.io import fits  ## fits 파일 읽기
from astropy.utils.data import get_pkg_data_filename  ## fits 파일 불러오기

위의 라이브러리들이 필수로 사용하는 아이들입니다. 

4-2. 파일 읽기

filename = get_pkg_data_filename('C:\\111\\IC434.fits')
hdu = fits.open(filename)[1]

FITS 파일을 읽기 위한 기본적인 코드입니다. 저기서 hdu=fits.open(filename)[1]은 조금 눈여겨 보셔야 하는데, FITS 파일은 여러 형태의 관측 데이터를 하나의 파일에 저장할 수 있는 장점이 있다 라고 말씀드렸습니다. 보통 [0]에는 헤더 데이터, 그러니까 관측 배경정보가, [1]에는 관측 데이터가 저장되있다 생각하시면 됩니다. 우리는 관측 결과를 뽑아올 것이기 때문에 1이라고 합니다.

4-3. 좌표와 그래프 그리기

wcs = WCS(hdu.header)
fig = plt.figure(figsize=(18, 12))
plt.subplot(projection=wcs)
plt.grid(color='black', ls='dotted')

FITS 파일에 관측 정보는 각 픽셀별 빛의 세기와 좌표가 저장되어 있습니다. 먼저 좌표를 불러옵니다. WCS는 world coordinate system의 약자인데, 일반적으로 알고있는 적경(RA), 적위(Dec)의 픽셀별 좌표가 여기에 모두 있습니다. 그래서 먼저 좌표정보를 불러오기 위해 wcs=WCS(hdu.header)라고 입력합니다.

 

그 다음은 그래프 크기이니 넘어가고,

 

plt.subplot(projection=wcs)인데, 사실 이번에 그리는 그래프는 엄밀히 말하면 2개의 그래프를 그리는 것입니다. 하나는 좌표, 하나는 좌표별 빛의 세기입니다. 그래서 2개 그래프를 그릴때 많이들 사용하는 plt.subplot를 사용하였고, projection으로 좌표인 wcs를 지정해 주었습니다. projection으로 wcs를 지정함으로써, matplotlib의 기본적인 좌표 사용방식을 따르지 않고, FITS 파일의 좌표형식을 따르겠다는 의미로 생각하시면 됩니다.

 

그리고 마지막에 이것을 plt.grid로 바둑판 형태로 불러옵니다.

여기까지 실행하면 위와 같은 그림이 나옵니다. 따로 x축 라벨과 y축 라벨을 지정하지 않았음에도, 저렇게 x축과 y축이 달려 나왔습니다. 이는 기본 좌표정보로 wcs를 사용했기 때문입니다. 위에 projection=wcs라고 쓴 이유이기도 합니다.

4-4. 픽셀별 에너지 세기를 컬러 이미지로 불러오기

bar=plt.imshow(hdu.data,  origin='lower', cmap=plt.cm.gist_heat)
plt.colorbar(bar, label='Herschel 250$\mu$m (MJy/sr)')

이제 각 픽셀별 빛의 세기를 색상 지도로 표현합니다. 첫 번째 imshow는 픽셀별 색상을 지정해 주는 것입니다. cmap은 표현할 색상의 유형을, vmin과 vmax는 각각 해당 색상에서의 최소값과 최대값입니다. 관측 대상의 구조를 더 잘 볼 수 있는 형태로 숫자를 조정하시면 됩니다.

plt.colorbar는 색상의 스케일을 보여주는 것입니다. 이렇게 하면 아래와 같은 그림이 똭 그려집니다.

 

앗. 뭔가 이미지가 튀어나왔습니다. 컬러 스케일을 보니 0을 검정색, 5000정도를 흰색으로 표현했습니다. 그러니 1000이나 2000정도의 값은 너무 작아서 굉장히 어둡게 표현되었습니다. 구조를 좀 잘 보기 위해서는 컬러 스케일 범위를 다르게 할 필요가 있습니다.

4-5. 컬러바의 스케일 변경하기

bar=plt.imshow(hdu.data,  origin='lower', cmap=plt.cm.gist_heat, vmin=0, vmax=2000)
plt.colorbar(bar, label='Herschel 250$\mu$m (MJy/sr)')

컬러바의 스케일 변경은 간단합니다. 위의 코드 첫 줄에 vmin에 최소값, vmax에 최대값을 입력하면 됩니다. 이렇게 하면 광량 0은 가장 어두운 검정색, 2000은 가장 밝은 흰색으로 표현합니다. 분명 저 데이터에는 2000보다 밝은 5000도 있습니다. 하지만 max값을 2000으로 잡았기에, 5000도 흰색이 되어버립니다.

 

이건 연구자의 나름으로 정하면 됩니다. 본인이 보고 싶은 구조 영역이 잘 안보인다면 그 범위가 잘 보이게만 하면 됩니다. 정답은 없습니다. 여기서는 연구목적이 아니니 전체적인 구조가 잘 보이도록 하였습니다. 만약 5000까지의 데이터가 모두 필요하다면, 절대 저렇게 vmax값을 잡아서는 안됩니다.

이제 그림이 완성되었습니다. 조금 이상합니다. 좌표값의 그리드가 검정색이라, 위에서는 선이 잘 보이지 않습니다. 그래서 plt.grid(colors='white')로 바꾸어 줍시다. 그리고 그래프 제목도 달아주고, x축 y축 이름도 정해주고 마지막으로 저장을 하면 됩니다. 이를 포함한 전체 코드는.

import matplotlib.pyplot as plt
from astropy.wcs import WCS  ##wcs-> world coordinate system
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename

filename = get_pkg_data_filename('C:\\111\\IC434.fits')
hdu = fits.open(filename)[1]
wcs = WCS(hdu.header)
fig = plt.figure(figsize=(18, 12))
plt.subplot(projection=wcs)
plt.grid(color='white', ls='dotted')
bar=plt.imshow(hdu.data,  origin='lower', cmap=plt.cm.gist_heat, vmin=0, vmax=2000)
plt.colorbar(bar, label='Herschel 250$\mu$m (MJy/sr)')
plt.xlabel('RA')
plt.ylabel('Dec')
plt.title('IC 434 dust continuum map(250'r'$\mu$m)', size=15)
plt.savefig('C:\\Users\\kalch\\Desktop\\IC434.png',bbox_inches='tight', pad_inches=0.1, dpi=300)
plt.show()

저기서 plt.savefig라던가 filename=변수뒤에 있는 경로는 본인 컴퓨터에 맞는 경로를 쓰셔야 합니다. 안그러면 계속 오류가 떠요. 어쨌든 완성한 결과를 보면

픽셀별로 광세기가 같은 곳을 선으로 연결 할 수도 있습니다. 

이런 식으로요.

 

이건 다음번 포스팅때 말씀드리도록 하겠습니다.