안녕하세요, 거의 1년만에 새로운 글로 찾아뵙게 되었습니다.
학생회장과 학부연구생, 대외활동으로 블로그에 글을 쓰지 못하였습니다.
작년 또한 마찬가지로 파이썬으로 천문학 공부를 해보았습니다. 좋은 데이터와 분석 코드가 있어 공유해드리고자합니다.
작년 1학기에는 천문학 데이터를 찾아보고 구글에 검색해보아서 나온 그래프를 똑같이 구현하고 그 의미를 분석하는 것이 저의 소소한 취미였습니다. 그러던 중 천체물리 수업에서 초기 우주에 대한 수업을 듣게 되었고, WMAP 데이터에 흥미가 가게 되었습니다.
그래서 곧바로 데이터를 서칭하여 열어보았더니, WMAP 데이터를 통하여 천문학 전용 라이브러리를 사용하지 않고 기초적인 파이썬 코딩으로도 데이터 분석이 가능하고 과학적으로 유의미한 결과를 낼 수 있었습니다. 그래서 흥미로운 주제로 파이썬을 통해 텍스트 파일을 다루어보고 싶거나, 천문학 데이터 분석에 입문하시는 분들을 위해서 제가 거쳐온 과정을 공유해드리고자 합니다. 그럼 시작하겠습니다~!!
(파이썬에서의 Class 문의 활용과 텍스트 파일을 열고 읽는 과정이 궁금하신 분은 제가 추후에 제작할 class 구문에 대한 게시물을 참고하시거나 조금 밑으로 내리셔서 코드블럭을 참고해주시길 바랍니다.)
https://himbopsa.tistory.com/50
[Python Class 구문 입문] Class 문을 활용하여 텍스트 형식의 WMAP 데이터 다루기 (천문학 데이터 분석
안녕하세요, AstroPenguin입니다! 2탄인 텍스트(.txt) 형식으로 된 WMAP데이터를 통하여 CMB의 Angular Power Spectrum 그래프를 그려낸 게시물로 먼저 찾아뵙고, 순서를 바꾸어 python으로 class구문을 작성하는.
himbopsa.tistory.com
WMAP (Wilkinson Microwave Anisotropy Probe) mission은 우주 마이크로파 배경 복사의 온도 분포를 전체 우주의 지도로 살펴봄으로써 우주의 기하학적인 내용 및 진화에 대하여 알아보는 임무입니다. 저희는 텍스트 형식의 WMAP 데이터를 파이썬을 통하여 분석해보고, 심화과정으로는 저희가 분석한 내용이 우주에서 뜻하는 의미에 대하여 이야기해볼 예정입니다.
우선 데이터를 다운받기 위하여 데이터 아카이브에 들어가겠습니다.
https://lambda.gsfc.nasa.gov/product/map/current/
LAMBDA - Wilkinson Microwave Anisotropy Probe
lambda.gsfc.nasa.gov
위의 링크를 타고 해당 페이지와 같은 곳으로 들어오셨다면 성공입니다.
그러면 아래 사진과 같은 과정을 거쳐 데이터를 다운받으시면 됩니다.
WMPA 데이터는 1~9년 동안의 데이터가 있으니 원하시는 데이터로 다운받으시면 되겠습니다. (저는 7년동안의 데이터를 사용하였습니다.)
Binned라고 쓰여진 데이터를 클릭하셔야 저희가 다루기 쉬운 데이터를 다운받을 수 있으니 유의 바랍니다!
다음과 같은 과정을 거치시면 데이터는 쉽게 다운받을 수 있습니다. 그렇다면 이제 데이터를 열어서 데이터에 대한 정보를 알아볼까요?
저희가 직접적으로 사용할 데이터는 1열의 multipole moment, 4열의 온도 변화에 대한 수치인 temperautre power spectrum, 7열의 에러 보정 수치입니다. 천문학 데이터는 보통 열어보면 관측기, 레퍼런스, 데이터 명과 목적, 각 열에 대한 정보 등 다루기 쉽게 설명이 지정되어있습니다. 그렇다면 지금부터 이 데이터를 파이썬을 통하여 분석해보도록 하겠습니다.
이 때 텍스트파일을 조금 가공하여야 합니다.
데이터 정보와 각 열에 대하여 적혀있는 정보들을 모두 지워주어야 합니다.
그 후 숫자로 나오는 데이터들 중 처음 데이터 값인 multipole moment = 2일 때의 수치를 지워주면 데이터 분석이 더욱 쉬워집니다. (첫 행에서 error값을 뺐을 때 음수가 나오게 됩니다. 제가 우주론 전공자가 아니라 그 이유는 자세히 모릅니다. 그래프를 그리는 실습을 위하여 이 점은 감안하고 넘어간 점 양해 부탁드립니다.)
저희가 최종적으로 이 데이터를 분석하여 그려야 할 그래프는 바로 다음과 같은 그림입니다.
해당 그래프가 뜻하는 의미와 왜 그려야 하는지는 다음 게시물을 통하여 설명해드리도록 하겠습니다. :)
우선 제가 사용한 class구문의 코드를 보여드리겠습니다.
import numpy as np
import matplotlib.pyplot as plt
class WMAP():
def __init__(self):
self.multipole_moment = []
self.TT_power_spec = []
self.power_err = []우선은 def __init__(self): 구문을 통하여 저희가 class 내에서 변수로 사용할 multiploe moment, power, error에 대해서 정의해줍니다. 우선은 빈 list로 두어 다음 구문에서 텍스트 파일을 읽어 각각에 대한 정보를 집어넣을 예정입니다. 이러한 방식은 제가 천문학 데이터를 분석 및 가공할 때 즐겨쓰는 방법입니다.
다음으로는 오늘의 하이라이트인 텍스트 파일을 파이썬으로 읽는 과정입니다.
file = open("파일명.txt", 'r') # 파일을 불러온다
header = file.readline() # 파일을 한 줄씩 읽어들인다
first_name = []
last_name = []
gender = []
adress = []
for line in file:
line = line.strip() # 맨 앞과 뒤의 공백을 없애준다
columns = line.split() # 각각의 열을 쪼개준다.
first_name.append(columns[0])
last_name.append(float(columns[1]))
gender.append(float(columns[4]))
adress.append(float(columns[6]))
f.close() # 열었던 파일을 닫는다신상정보에 대한 텍스트 파일이 있다는 가정하에, 다음과 같이 파일을 열고 한 줄 씩 불러들입니다. 그리고 원하는 정보들에 대한 빈 list를 선언해줍니다.
그 후 .strip()구문으로 맨 앞과 뒤의 공백을 없애준 후 .split 구문으로 각각의 열을 다른 list로 쪼개줍니다. 마지막으로 각각의 정보를 나타내는 빈 list에 원하는 정보를 append해준 후 열었던 파일을 도로 닫아주면 됩니다. 다만 저희의 wamp 데이터를 활용하는 코드는 class 내부에 있기 때문에 readflie(self, filename)이라는 새로운 함수를 만들어주었습니다.
텍스트로 된 wmap데이터를 읽어오는 코드를 추가한 class 구문은 다음과 같습니다.
import numpy as np
import matplotlib.pyplot as plt
class WMAP():
def __init__(self):
self.multipole_moment = []
self.TT_power_spec = []
self.power_err = []
def readfile(self, filename): # 텍스트 형식의 파일을 여는 코드
try:
file = open(filename,'r')
del_2nd_row = file.readline()
for line in file:
line = line.strip()
columns = line.split()
self.multipole_moment.append(float(columns[0]))
self.TT_power_spec.append(float(columns[3]))
self.power_err.append(float(columns[6]))
file.close()
except:
print('FileNotFoundError : Please check your filename')이 곳에서는 try-except 구문을 사용하여 파일을 찾을 수 없어 에러가 뜨는 경우에는 파일 이름을 다시 체크해달라는 에러메세지가 뜨도록 설정하였습니다.
파일을 불러들여오고, 원하는 정보들인 multipole moment, temperature power spectrum, error값들에 대하여 선언 및 수치들에 대한 append가 끝났으니, 이제 해당 정보들을 가공하는 과정을 거쳐주겠습니다. temperature power spectrum에서 error값을 빼주어 오차에 대한 보정을 진행해준 후 multipole moment - temperature power spectrum의 관계를 보여주는 그래프를 그려내는 함수를 추가하였습니다. 그리고 이러한 과정을 multipole moment에 log scale을 적용하여 다시금 진행하였습니다. 그리하여 나오게 된 전체 코드는 다음과 같습니다.
import numpy as np
import matplotlib.pyplot as plt
class WMAP():
def __init__(self):
self.multipole_moment = []
self.TT_power_spec = []
self.power_err = []
def readfile(self, filename): # 텍스트 형식의 파일을 여는 코드
try:
file = open(filename,'r')
del_2nd_row = file.readline()
for line in file:
line = line.strip()
columns = line.split()
self.multipole_moment.append(float(columns[0]))
self.TT_power_spec.append(float(columns[3]))
self.power_err.append(float(columns[6]))
file.close()
except:
print('FileNotFoundError : Please check your filename')
def real_WMAP_power(self): # 오차를 빼서 정확한 값을 찾는 코드
observated_value = np.array(self.TT_power_spec)
error_value = np.array(self.power_err)
real_value = observated_value - error_value
return real_value
def Angular_Spectrum_plot(self): # 그래프를 그리는 코드
x = np.array(self.multipole_moment)
y = np.array(self.real_WMAP_power())
self.size = 20
plt.figure(figsize=(9,7))
plt.errorbar(x,y,yerr=self.power_err,c='red')
plt.scatter(x,y,c='black',s=30,marker='d')
plt.title('The CMB Power Spectrum', fontsize = self.size)
plt.xlabel('Multipole Moment (l)', fontsize = self.size)
plt.ylabel(r"$\Delta T$" +" " + r"[$\mu K^2$]", fontsize = self.size)
plt.xticks([10, 100, 500, 1000], fontsize = self.size)
plt.yticks([1000*(y+1) for y in range(6)], fontsize = self.size)
plt.show()
# r"{$l(l+1)\C_1/2\pi$}"
def log10_Angular_Spectrum_plot(self): # 로그 스케일의 그래프를 그리는 코드
x = np.log10(np.array(self.multipole_moment))
y = np.array(self.real_WMAP_power())
self.size = 20
plt.figure(figsize=(9,7))
plt.errorbar(x,y,yerr=self.power_err,c='blue')
plt.scatter(x,y,c='black',s=30,marker='d')
plt.title('The Log-Scale CMB Power Spectrum', fontsize = self.size)
plt.xlabel('Log Multipole Moment (log l)', fontsize = self.size)
plt.ylabel(r"$\Delta T$" +" " + r"[$\mu K^2$]", fontsize = self.size)
plt.xticks([0.5*(x+1) for x in range(6)], fontsize = self.size)
plt.yticks([1000*(y+1) for y in range(6)], fontsize = self.size)
plt.show()
# r"{$l(l+1)\C_1/2\pi$}"
WMAP = WMAP()
WMAP.readfile('WMAP_TT_power_spec.txt')
print(WMAP.real_WMAP_power())
WMAP.Angular_Spectrum_plot()
WMAP.log10_Angular_Spectrum_plot()위와 같은 과정을 거쳐 나오게 되는 그래프는 아래의 이미지와 같습니다.
그래프들을 보시면 제가 만들고자 하였던 이미지와 많이 다른 것을 알 수 있습니다. 이는 제가 그린 빨간색 그래프에서는 x축 값인 1, 10, 100, 1000의 간격이 실제 숫자 간격에 맞게 배치되어있으나 궁극적으로 만들고자 하였던 그래프에서는 1, 10, 100, 1000의 간격이 모두 똑같게 배치되어있기 때문입니다.
다음 게시물에서는 이를 극복하는 과정과 제가 그려낸 그래프가 우주론에서 어떠한 의미를 가지고 있는지에 대하여 간단히 서술해드리도록 하겠습니다.
긴 글 읽어주셔서 감사드립니다 :)