찐따의 파이썬 독학 - Astropy를 활용한 천문학 데이터 FITS 영상 처리

Astropy를 활용한 FITS 데이터 영상 처리

안녕하세요, 고졸 찐따인 흔한 찐따입니다.

요즘 천문학계에서는 제임스 웹 망원경으로 굉장히 핫한데,

그전부터 천문학에 관심이 많이 있었기도 했고, 한번 천문학 데이터를 다뤄보고 싶었습니다.

그래서 혼자 공부를 하면서 이 문서를 작성하게 되었습니다.

 

개요

본 문서에서는 astropy.utils.data 를 사용하여 데이터 파일을 다운로드한 다음,

astropy.io.fits 를 사용하여 파일을 열고

마지막으로 matplotlib 라이브러리를 사용하여 다양한 색상 스케일 및 스트레치로 이미지를 생성하고

히스토그램을 만드는 방법을 보여줍니다.
이 문서에서는 간단한 이미지 스태킹의 데모도 포함했습니다.

 

위와 같은 과정을 거치기 위해서 아래의 절차대로 진행합니다.

1. FITS 파일 열기 및 이미지 데이터 로드
2. 이미지 데이터로 2D 히스토그램 만들기
3. 여러 이미지를 단일 이미지로 쌓기
4. FITS 파일에 이미지 데이터 쓰기

다음의 공식 문서를 통하여 자세한 FITS 데이터를 불러오는 API 정보를 확인할 수 있습니다.

• 공식문서: FITS File Handling (astropy.io.fits)

 

Astropy 설치하기

Astropy 라이브러리를 설치하기 위해서는 다음과 같이 명령 프롬프트에 명령어를 입력합니다.

pip install astropy

 

Astropy 라이브러리에 대한 자세한 설명을 참고하려면 아래의 논문을 참고하시면 됩니다.

• 논문: Astropy: A community Python package for astronomy

 

FITS 파일 불러오기

FITS 파일

FITS 파일 데이터셋은 NASA에서 제공하는 오픈 데이터를 사용합니다.

 

앞서 설명한 바와 같이 matplotlib 라이브러리를 통해 다양한 색상 스케일 및 스트레치로 이미지를 생성하고 히스토그램을 만들기 위해서 다음과 같이 import 합니다.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

먼저, FITS 파일(말머리성운에 대한 이미지 데이터 포함)을 다운로드합니다.

 

프로젝트를 진행하기 위해서는 FITS 파일 이미지를 다운로드해야 합니다.
사진판으로 찍은 말머리성운의 천체 사진입니다.
참고로, 이미지가 디지털화되어 있습니다.
즉, 다시 말해 이미지가 컴퓨터로 스캔되어 2차원 배열로 변환되어 있습니다.
배열의 각 위치는 하늘에 투영된 위치에 해당하며 이미지의 밝은 영역은 높은 값을 갖고 어두운 영역은 낮은 값을 갖습니다.

 

CCD 또는 "전하 결합 장치"라고 하는 천체 장비로 촬영한 이미지도 유사하게 구성됩니다.
빛에 의해 조명되면 CCD는 전자를 축적하여 밝기 값을 전자 수로 변환합니다.
CCD 이미지는 기본적으로 배열의 각 위치가 단일 CCD 픽셀을 나타내고 해당 배열의 값이 해당 픽셀에 등록된 개수를 나타내는 2차원 배열입니다.

 

CCD Data를 사용하는 방법에 대해서는 아래의 논문을 확인하시면 됩니다.

• “A User’s Guide to CCD Reductions with IRAF” (Massey 1997)
• “A User’s Guide to Stellar CCD Photometry with IRAF” (Massey and Davis 1992).

 

혹은 아래의 동영상을 통해 CCD Data를 사용하는 방법을 참고하면 좋을 것 같습니다.

 

AAVSO How-to Hour: How to Start with CCD Photometry

 

from astropy.io import fits
from astropy.utils.data import download_file

url_fits_data = "http://data.astropy.org/tutorials/FITS-images/HorseHead.fits"
image_file = download_file(url_fits_data, cache=True )

불러온 데이터를 사용하려면 astropy.io.fits.open() 을 사용하면 됩니다.

 

FITS 파일 열기 및 이미지 데이터 로드

FITS 파일을 열어 내용을 알아보겠습니다.

hdu_list = fits.open(image_file)
hdu_list.info()

결과

Filename: C:\Users\iamjjintta\.astropy\cache\download\url\ff6e0b93871033c68022ca026a956d87\contents
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     161   (891, 893)   int16   
  1  er.mask       1 TableHDU        25   1600R x 4C   [F6.2, F6.2, F6.2, F6.2]

일반적으로 이미지 정보는 PRIMARY 블록에 위치합니다.
블록에는 번호가 매겨져 있으며, hdu_list 변수를 인덱싱을 통해 액세스 할 수 있습니다.

image_data = hdu_list[0].data

위의 데이터는 이제 2D numpy 배열로 저장됩니다.

numpy 의 shape 프로퍼티를 통해서 배열의 크기를 확인할 수 있습니다.

print(type(image_data))
print(image_data.shape)

결과

<class 'numpy.ndarray'>
(893, 891)

이 시점에서 원하는 모든 것을 변수에 저장했기 때문에 FITS 파일을 close 할 수 있습니다.

FITS 파일을 close 함으로써 컴퓨터의 과도한 메모리나 파일 핸들링을 계속 사용하지 않도록 할 수 있습니다.

hdu_list.close()

 

더 빠르게 불러오기

FITS 헤더를 검사할 필요가 없는 경우 fits.getdata 를 호출하여 이전 단계를 건너뛸 수 있습니다.

image_data = fits.getdata(image_file)

위의 방식대로 불러오면 이미지 데이터 역시 2D numpy 배열에 저장됩니다.

2차원 배열 모양을 print 하여 이미지의 픽셀 수를 볼 수도 있습니다.
확인해보면 이미지의 크기가 893 x 891 픽셀임을 알 수 있습니다.

print(type(image_data)) # Show the Python type for image_data
print(image_data.shape) # Show the number of pixels per side in the 2-D image

결과

<class 'numpy.ndarray'>
(893, 891)

 

이미지 데이터 보기 및 기본 통계정보 가져오기

plt.imshow(image_data, cmap='gray')
plt.colorbar()

결과

<matplotlib.colorbar.Colorbar at 0x1b502c37400>

이미지에 대한 몇 가지 기본적인 통계 정보를 살펴보겠습니다.

print('Min:', np.min(image_data))
print('Max:', np.max(image_data))
print('Mean:', np.mean(image_data))
print('Stdev:', np.std(image_data))

결과

Min: 3759
Max: 22918
Mean: 9831.481676287574
Stdev: 3032.3927542049046

 

히스토그램 플로팅

matplotlib.pyplot.hist() 로 히스토그램을 만들려면 2D 배열의 데이터를 1차원으로 변환해야 합니다.

즉, 다시 말해 1차원 벡터로 치환하는 과정이 필요합니다.
이 경우, ndarray.flatten() 를 사용하여 1D numpy 배열을 반환해 보겠습니다.

여기서 print 함수를 사용하여 "평평한" 배열이 여전히 numpy 배열인지,

그리고 평면화 된 배열의 항목 수가 이미지의 총 픽셀 수와 같은지 확인합니다.

(즉, 픽셀 총계; 893 x 891 = 795663 )

print(type(image_data.flatten()))
print(image_data.flatten().shape)

결과

<class 'numpy.ndarray'>
(795663,)

이제 matplotlib.pyplot.hist 을 사용하여 플롯을 만들면 됩니다.

histogram = plt.hist(image_data.flatten(), bins='auto')

결과

 

로그 스케일로 이미지 표시

로그 색상 스케일을 사용하려면 matplotlib 에서 LogNorm 객체를 생성하면 됩니다.

from matplotlib.colors import LogNorm

plt.imshow(image_data, cmap='gray', norm=LogNorm())

# I chose the tick marks based on the histogram above
cbar = plt.colorbar(ticks=[5.e3,1.e4,2.e4])
cbar.ax.set_yticklabels(['5,000','10,000','20,000'])

결과

[Text(1, 5000.0, '5,000'),
 Text(1, 10000.0, '10,000'),
 Text(1, 20000.0, '20,000')]

 

기본적인 영상 처리: 이미지 스태킹

다른 numpy 배열과 마찬가지로 영상 데이터를 사용하여 계산을 수행할 수도 있습니다.
본 문서에서는 최대 10인치 망원경으로 촬영한 M13의 여러 이미지를 쌓을 것입니다.

먼저, 일련의 FITS 파일을 열고 image_concat 이라는 이름의 리스트에 데이터를 저장해 보겠습니다.

base_url = 'http://data.astropy.org/tutorials/FITS-images/M13_blue_{0:04d}.fits'

image_list = [
    download_file(base_url.format(n), cache=True)
    for n in range(1, 5+1)
]
image_concat = [fits.getdata(image) for image in image_list]

이제 연결된 리스트를 합산하여 이미지를 쌓습니다.

# The long way
final_image = np.zeros(shape=image_concat[0].shape)

for image in image_concat:
    final_image += image

# The short way
# final_image = np.sum(image_concat, axis=0)

바로 이미지를 확인할 수 있지만, 그전에 가장 좋은 스트레치를 결정해야 합니다.
이렇게 하려면 먼저 데이터의 히스토그램을 표시합니다.

image_hist = plt.hist(final_image.flatten(), bins='auto')

결과

vmin 및 vmax 파라미터를 사용하여 imshow 의 색상 스케일링 제한을 설정합니다.
위의 히스토그램에 따라 한계를 설정함으로써 보다 뚜렷한 색상 분포를 가진 이미지를 생성할 수 있습니다.

plt.imshow(final_image, cmap='gray', vmin=2E3, vmax=3E3)
plt.colorbar()

결과

<matplotlib.colorbar.Colorbar at 0x1b506927370>

 

FITS 파일에 이미지 데이터 쓰기 및 불러오기

writeto() 메서드를 사용하여 이 작업을 수행할 수 있습니다.

 

주의: 쓰려는 파일이 이미 있으면 오류가 발생합니다.
그래서 overwrite=True 로 설정했습니다.

 

outfile = 'stacked_M13_blue.fits'

hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, overwrite=True)

저장된 데이터를 불러온 후 확인해봅니다.

stacked_M13_blue = fits.open(outfile)
M13_data = stacked_M13_blue[0].data

plt.imshow(M13_data, cmap='gray', vmin=2E3, vmax=3E3)
plt.colorbar()

stacked_M13_blue.close()

결과

 

참고

• Astropy: A community Python package for astronomy
• Astropy - FITS File Handling (astropy.io.fits)
• Definition of the Flexible Image Transport System (FITS)
• NASA Open Science Data - MISR Level 1A CCD Science data, all cameras (MIL1A_V1)
• The AAVSO Guide to CCD Photometry
• “A User’s Guide to CCD Reductions with IRAF” (Massey 1997)
• “A User’s Guide to Stellar CCD Photometry with IRAF” (Massey and Davis 1992).
• YouTube - AAVSO How-to Hour: How to Start with CCD Photometry