拿太阳做壁纸-用Python的astropy和matplotlib库处理FITS文件

如果你想要

请DIY吧

——苑长

• 机缘巧合🍀
• 数据的获取📑
  • 让我们先Bing一下⌨️
  • 拿下数据🗄
• 开始处理吧
  • 再Bing一下🔔
  • Astropy库🔭
  • 你好~太阳☀️
  • 修改色阶🌈
  • GIMP处理
• 后记

⚓️机缘巧合🍀

之前网上闲逛的时候，突然想起 Edge for Chrome 浏览器的 网页App功能 ，又捣鼓了一番。而我之前就装了微博（主要是看航天的微博），便打开测试一下。

于是就看到了一篇微博：新科学家杂志上刊登的高清太阳日冕图，一下子就被这篇文章的封面的照片吸引住了。

配色醒目，细节突出，层次感非常好，而且还是太阳（也许是一种对太阳的特殊感觉,赞美太阳!）。我便想着搞到高清的图片做壁纸，但是一番搜索，并没有找到（原论文也没有）高清的大图。

由于有这方面的视界，顺其自然地便想着自己着手做了。因为大部分的天文数据都是开放的，一些处理数据的图像软件也是有的。于是和 Bing一起<(￣︶￣)↗[GO!]

⚓️数据的获取📑

有了微博的那张图(就是文章开头的那张图啦)作参考，直接搜索引擎就好啦，基本就是先找到新闻源，然后找到这个项目的相关信息，最后下到数据，作图即可。

⚓️让我们先Bing一下⌨️

可以得知是新科学家杂志发表的一篇新闻，那么直接检索 新科学家杂志 太阳来看看有没有有用的消息。发现好像并没有太多相关的内容,那毕竟是外网的内容,我们直接英文走起。

检索 新科学家杂志就可以从百度百科里找到它的英文名 New Scientist(这里建议还是要搜一下的, 毕竟专有名词可能会比较特殊), 直接检索 New Scientist sun我们就可以看到结果啦。

点开后直接就可以看到这篇文章的封面图了，我这就只放个链接吧。

Weird magnetic threads in sun's corona seen for the first time Read more

文章带的图片清晰度不够，所以还是要尝试自己做。另外可惜的是要订阅才能看到全文

收费的不要

看来得另辟蹊径。不过仔细看了一下图片的来源写着 NASA High Resolution Coronal Imager，有了些眉目，这应该就是来源了。直接一个 Bing检索这个词，结果非常 amazing.

基本可以得知这个是缩写为 Hi-C项目的成果了,点开右排第一个链接一探究竟.

一通阅读，可以得知这是一个 探空火箭项目🚀 ，通过探空火箭 Black Brant IX发射 Hi-C太空望远镜，把 Hi-C送入亚轨道，进行短暂的 弹道飞行(考虑了空气阻力的抛物线轨迹)并在 滑行段(只受引力阶段)对太阳成像。相当于用火箭大力地把望远镜从地球表面 抛起，望远镜可以依靠惯性冲到百公里高度。 这个高度已经脱离 浓密大气层了(事实上绝大部分气象和主要影响可见光观测的大气成分都集中在地表以上几公里范围内的 对流层；100km高度为国际公认的太空界面——卡门线)

The experiment reached a maximum velocity of Mach 7 and max altitude of 264 km. The experiment collected 345 seconds of EUV science images.

实验达到的最大速度为七马赫(七倍音速)，最大高度264千米，收集了345s的极紫外科学图像

资料来源：Media Resources: NASA Telescope’s Sharpest-Ever Images of Corona Catch Glimpse of Energy Source in the Solar Atmosphere

下图来自：Hi-C 项目NASA页面的探空火箭介绍子菜单

我做了简单的注释，原链接里还有更多的项目实施相关的图片。

望远镜的全称叫 The High-Resolution Coronal Imager，也是该项目的名字，通过检索这个名字找到了一篇期刊文章：The High-Resolution Coronal Imager, Flight 2.1 里面介绍了望远镜的情况。下图来自这篇文章，图解了望远镜的结构和实验组件：

Hi-C火箭发射视频:

既然项目找到了，还有了大致的了解了，我们继续吧。

⚓️拿下数据🗄

简单地浏览主页的菜单，第二个一级菜单就是Hi-C DATA PRODUCTS——Hi-C 数据产品看这名字就知道是我们想要的了。

二级菜单介绍：

• DATA POLICY & PRODUCTS：

  • Hi-C项目的简要介绍,引用数据规范，数据库链接
  • Hi-C 1 ：数据库
  • Hi-C 2.1：数据库

• DATA ANALYSIS RESOURCES

  • 官方的数据使用介绍，图像处理手册

获取图像只要点击数据库链接就好啦(通过数据库的介绍图片，可知是 Hi-C 2.1是我们找的照片的数据源)。

有用户手册，也有各种不同的打包文件。倒数第二个是 AIA和 IRIS的观测数据，应该是用于和 Hi-C拍摄的图像做对比分析的。这里可以随便选择包进行下载，也无所谓哪个观测设备拍的，毕竟只是做壁纸的话，对数据的要求并不高。这里我直接下了 AIA的数据了(介绍了半天居然没用Hi-C的数据，哈哈哈哈)

相关项目的介绍：

1.AIA:太阳大气成像仪（Atmospheric Imaging Assembly）,为SDO卫星的一个成像部件。

2.IRIS:表面区域成像光谱仪（the Interface Region Imaging Spectrograph），SMEX卫星。

⚓️开始处理吧

下载到的格式为 FITS,是个没见过的格式，这个时候又到了 Bing的时候啦。

不过英文好的话，也可以直接阅读数据库里的英文用户手册，按照上面的做应该也是可以处理的。但是考虑到涉及一些专业的软件，这里我们还是通过搜索引擎来帮助我们，看看有没有更适合小白的方法 :happy: 。

⚓️再Bing一下🔔

检索词就 FITS、处理、格式、Python之类的，相信可以得到比较好的结果。

这是搜到的较详细的 FITS格式的介绍：在Python中FITS格式文件数据的读取

这篇文章详细的介绍了FITS格式的结构内容，以及一个叫做Astropy的处理天文数据的Python库。

一个崭新的大门出现了，一个新的世界出现了。

⚓️Astropy库🔭

由于是 Python的第三方库，需要安装 Python, 可以参考这篇文章 Anaconda之路

如果是直接安装的 Python并且加入了环境变量的话, 按 Win+R, 在弹出的框框里输入 cmd，回车,在弹出的 cmd界面输入：

pip install astropy

安装搞定！

⚓️你好~太阳☀️

通过查阅Astropy的文档，有专门的 FITS 格式的读取模块

Win+R,在弹出的框框里输入 cmd后,在弹出的命令提示符界面输入 ipython调出交互式python

(下面的代码框里的 #后面的字符为注释, CTRL+C时无需复制)

参照着文档的描述，开始读取 FITS文件。以下内容来自我写的Astropy处理FITS.py，源文件可直接在左边蓝色的字符上右键另存为即可。主要是我个人的理解和对官方文档的翻译，所以为了方便大家对比，这里推荐一位大佬编写的Astropy的中文文档 ，右键另存为即可。

# astropy库fits模块读取FITS文件
# 参考翻译自：https://docs.astropy.org/en/stable/io/fits/index.html
# 源码下载：https://github.com/astropy/astropy
# 用 # 标注的主要是我对代码的注释
# 用 ''' ''' 包围的主要是官方的说明和代码对应的输出结果
# 这两种标记都是用来说明的
#-------------导入库---------------

# 从astropy库中导入‘io.fits’模块
from astropy.io import fits  

#-------------打开文件---------------

# 这个方法可以直接在线调用官方的示例图片,这里用不上，注释掉
# fits_image_filename = fits.util.get_testdata_filepath('test0.fits')  

# 我的本地FITS文件的文件名字符串赋给一个变量，方便后面调用
fits_image_filename = 'aia.lev1.2018-05-29T185610Z.172929830.image_lev1.fits'   
# 调用astropy.io.open()打开文件，这里终端打开的文件夹就是图片所在位置，所以直接调用文件名就可以打开
# 也可以使用如下的open()参数，mode='update'，实时更新
'''代码示例
with fits.open('original.fits', mode='update') as hdul:
    # Change something in hdul.
    # 对hdul做些修改
    hdul.flush()  # 把修改保存到原文件 changes are written back to original.fits
# 关闭文件将会保存修改，并阻止后面的写入
# closing the file will also flush any changes and prevent further writing
'''
# 默认是'只读'模式读取文件
hdul = fits.open(fits_image_filename) 

# 这里提供了另一种方式来读取文件
# with层级下的代码执行完后，会自动回收，也就可以直接关闭文件了
'''代码示例
with fits.open(fits_image_filename) as hdul:
     hdul.info()
     ···
'''

print(type(hdul))
'''输出结果
<class 'astropy.io.fits.hdu.hdulist.HDUList'>
'''

#-------------读取文件--------------
''' [HDUList的说明]
hdul是一个HDUList对象，即一个看起来像列表的 HDU 对象的集合. 
an HDUList which is a list-like collection of HDU objects.

HDU (Header Data Unit，头部数据集合) 是FITS文件构造中最高层级部分,
is the highest level component of the FITS file structure, 

由头部数据和数据矩阵或者是数据表组成。
consisting of a header and (typically) a data array or table.
'''
# 调用HDUList.info()方法来查看文件概况
hdul.info()
'''输出结果
Filename: aia.lev1.2018-05-29T185610Z.172929830.image_lev1.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       6   ()
  1                1 CompImageHDU    190   (4096, 4096)   int16
'''
# 可见这个文件包含有两个HDU数据集，索引号为0的那个HDU，有六个说明用的Cards
# 维度为空，所以初步推断是文件的说明部分。
# 索引号为1的那个HDU，有190个Cards,说明含有很多参数
# 维度为(4096,4096),数据类型为int16，这应该就是我们想要的太阳的图片数据了
# 而且Type属性也说了第二个HDU是CompImageHDU，石锤。



# astropy.io.open()还可以注入其他参数实现一些其他的功能
# 比如lazy_load_hdus参数，若设为'True'，文件被close()后还可以读取，数据的头部数据
'''代码示例
hdul = fits.open(fits_image_filename,lazy_load_hdus=False) 
'''

# 导入大文件时，memmap参数，可以通过内存映射的方式而不一次性把文件载入内存
# 默认为'True'，不过有一个问题就是文件close()后，可能还有残留的读取句柄在读取文件
# 参考：https://docs.astropy.org/en/stable/io/fits/index.html#working-with-large-files
'''代码示例
hdul = fits.open(fits_image_filename,memmap=True) 
'''

# FITS采用的是Fortran的风格的类型系统，不支持无符号的整数数据
# astropy导入FITS文件会自动转换有符号的整数数据为无符号的整数数据，uint=False参数来控制
# 参考：https://docs.astropy.org/en/stable/io/fits/index.html#unsigned-integers
''' 代码示例
hdul = fits.open(fits_image_filename,uint=False) 
'''

# astropy的open()可以无障碍地打开gzip，bzip2，pkzip等压缩的FITS文件
# 不过不同格式，有一些不同的限制。对于zip，open方法只会打开压缩包里的第一个FITS文件
# 对于bzip2则不支持使用append和update模式来open。

# 调用第一个头部数据集(HDU)的头部数据中的'COMMENT'这个Card的数据
# 也可以使用索引号，hdul[0].header[4] ，不建议使用，不同文件的特定关键字位置不一样
hdul[0].header['COMMENT']
'''输出结果
FITS (Flexible Image Transport System) format is defined in 'Astronomy
and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
'''
hdul[0].header[4]
'''输出结果
FITS (Flexible Image Transport System) format is defined in 'Astronomy
'''

# 修改FITS文件，内置了一个评论和历史记录系统，也有些意思
# 参考：https://docs.astropy.org/en/stable/io/fits/index.html#working-with-fits-headers
# 将header传给hdr，方便后面修改属性，添加了一条历史记录和一条评论
hdr = hdul[0].header
'''代码示例，这里注释掉，我已经添加过了，就不再添加了
hdr['history'] = 'YuanZ1949 changed this file on 2020.6.12'
hdr['comment'] = 'Thanks STEM, I am very happy to know the power of sun.'
'''
hdr  
# print(repr(hdr)) # 这个也可以
'''输出结果
SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                   16 / number of bits per data pixel
NAXIS   =                    0 / number of data axes
EXTEND  =                    T / FITS dataset may contain extensions
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
COMMENT Thanks STEM, I am very happy to know the power of sun.
HISTORY YuanZ1949 changed this file on 2020.6.12
'''

'''代码示例，修改其他属性
hdr = hdul[0].header

hdr.set('observer', 'Edwin Hubble')

hdr['targname'] = ('NGC121-a', 'the observation target')
hdr.comments['targname']
'''

hdr[:2]
'''输出结果
SIMPLE  =                    T / file does conform to FITS standard   
BITPIX  =                   16 / number of bits per data pixel  
'''
list(hdr.keys())  
'''输出结果
['SIMPLE',
 'BITPIX',
 'NAXIS',
 'EXTEND',
 'COMMENT',
 'COMMENT',
 'COMMENT',
 'HISTORY']
'''

# 如果读取的FIST文件不是非常标准的FIST文件(这是常有的事)
# 那可以使用fix工具尝试修复(仅限比较无关的差异)
'''读取报错时的输出
VerifyError: Unparsable card (OSCNMEAN), fix it first with .verify('fix').
'''
# 修复代码
hdul[1].verify('fix')
'''输出结果
WARNING: VerifyWarning: Verification reported errors: [astropy.io.fits.verify]
WARNING: VerifyWarning: Card 68: [astropy.io.fits.verify]
WARNING: VerifyWarning:     Card 'OSCNMEAN' is not FITS standard (invalid value string: '-nan').  Fixed 'OSCNMEAN' card to meet the FITS standard. [astropy.io.fits.verify]
WARNING: VerifyWarning: Card 69: [astropy.io.fits.verify]
WARNING: VerifyWarning:     Card 'OSCNRMS' is not FITS standard (invalid value string: '-nan').  Fixed 'OSCNRMS' card to meet the FITS standard. [astropy.io.fits.verify]
WARNING: VerifyWarning: Note: astropy.io.fits uses zero-based indexing.
 [astropy.io.fits.verify]
'''
# 读取带有图片数据的第二个HDU(hdul[1])的数据部分(.data)
data = hdul[1].data
# 还可以使用name来索引hdu，但是我们的文件里第二个HDU没有这个属性，就不用了
'''代码示例
data = hdul['某个name值'].data
data = hdul['某个name值',2].data  # 同名name的第二个HDU
'''

# 获取数据矩阵大小
data.shape
'''输出结果
(4096, 4096)
'''
# 数据的类型
data.dtype.name
'''输出结果
'int16'
'''

# 可以看出是numpy类型，那么就可以随意分割，查看，数学运算
type(data)
'''输出结果
numpy.ndarray
'''

# 取值,取指定行列的数据的值
print(data[1, 4])
'''输出结果
0
'''

# 切片,截取数据的一部分输出出来
data[30:40, 10:20]
'''输出结果
array([[ 0,  3,  0,  0, -2,  0,  1,  0,  0, -1],
       [ 0,  0, -1, -1,  1, -1,  0,  1, -1,  1],
       [ 0, -1, -1, -1,  0,  0, -1,  1, -1, -1],
       [-1,  0,  0,  0, -2,  2,  0,  0, -2,  0],
       [ 1,  0, -2,  1,  0,  0,  0,  1, -1, -1],
       [-2, -1,  0, -1, -1, -2, -1,  0, -1, -1],
       [ 0, -1, -2,  0,  0, -1,  0,  1, -1,  1],
       [-3,  1,  0,  1,  1,  1, -1, -1, -2,  1],
       [ 0, -2, -2, -1,  0, -2,  1,  1,  1,  1],
       [ 0, -1, -1, -2,  0,  1,  1,  0,  2,  0]], dtype=int16)
'''
# The next example of array manipulation is to convert the image data from counts to flux:
'''不太懂的一个操作，可能是特殊图片的处理
photflam = hdul[1].header['photflam']
exptime = hdr['exptime']
data = data * photflam / exptime
hdul.close()
'''

# Table数据处理(那种每个数据点是个元组的FITS文件)
# 参考：https://docs.astropy.org/en/stable/io/fits/index.html#working-with-table-data
# 暂且略过

# 保存文件更改，调用HDUList.writeto() 方法保存
'''代码示例:#这里注释掉,我下载的hdul不规范,无法保存
hdul.writeto('newtable.fits')
'''

# 读完数据自然要关掉文件，调用close()即可
hdul.close()

上面把图像数据赋值给了 data这个变量，并演示了一些处理矩阵的函数。

---

下面就开始处理 data, 图像处理参考：参考文章1;参考文章2(前面文章的汇总和翻译)

# 导入matplotlib绘图库
from matplotlib import image
 # 调用图片保存函数直接将数组储存为图片
image.imsave('FITS-image-out.png', data)

效果：(分辨率为原图片,想要可自行保存)

其实这里的上色采用的是默认的色阶上色的,而且这样直接保存，自定义的空间就太少了。

不过这个保存图片的 imsave函数还是可以导入几个参数的，也算是弥补了一些不足。

def imsave(fname, arr, vmin=None, vmax=None, cmap=None, format=None,
           origin=None, dpi=100):
    #Save an array as an image file.

比如里面的 camp就可以调节上色的色阶。

---

# 导入pyplot库,取plt做别名
from matplotlib import pyplot as plt
#plt绘制data
plt.imshow(data)
#plt显示绘制结果
plt.show()

可以看到下方的函数说明里，imshow明显比 imsave要多好多的参数，这就提供了相当多的自定义空间，不过这里还是用的初始化的设置，所以两者的导出效果是一样的。

imshow的使用可以参考这篇文章：Matplotlib imshow()函数

def imshow(
        X, cmap=None, norm=None, aspect=None, interpolation=None,
        alpha=None, vmin=None, vmax=None, origin=None, extent=None,
        shape=cbook.deprecation._deprecated_parameter, filternorm=1,
        filterrad=4.0, imlim=cbook.deprecation._deprecated_parameter,
        resample=None, url=None, *, data=None, **kwargs):

效果：

---

# 导入opencv库
import cv2
# 保存data
cv2.imwrite("FITS-cv2-out.png", data)

效果：(此为原图，可直接右键保存)

⚓️修改色阶🌈

cmap其实是 colormap的缩写，指代的意思是色彩映射图。

有关于不同的 cmap可以参考：matplotlib.pyplot——cmap直观理解

通过不同的色阶设置，可以得到各种各样的配色的图片。

⚓️GIMP处理

GIMP:the GNU Image Manipulation Program
GNU不是Unix(GNU’s Not Unix) 图像 处理 程序

a cross-platform image editor available for GNU/Linux, OS X, Windows and more operating systems.
一个跨平台的可用于 GNU/Linux, OS X, Window和其他操作系统的图像编辑器

Whether you are a graphic designer, photographer, illustrator, or scientist, GIMP provides you with sophisticated tools to get your job done. You can further enhance your productivity with GIMP thanks to many customization options and 3rd party plugins.
无论您是平面设计师，摄影师，插图画家还是科学家，GIMP都能为您提供完善的工具来完成您的工作。 借助许多自定义选项和第三方插件，您可以使用GIMP进一步提高生产力。

总的来说就是 GIMP是一个专业的开源免费的图像处理软件。

---

我们取上文的 黑白着色和 默认着色的两张导出图片，来尝试自己做壁纸。

黑白着色可见保留了相当多的中间亮度的细节。

默认着色则保留了暗部和亮部的一些细节。

通过将 黑白着色叠在 默认着色的上面，然后 黑白着色采用 叠加作为图层的模式，就可以较好的保留两者的细节。之后再调节 明度曲线、色相、曝光，优化得到的最后结果啦。

⚓️后记

最后的效果就是这样啦！

前前后后拖了蛮久，时值暑假之际，静下心来把这篇补的差不多了。

还有一些日后再补，主要是自定义 camp，这样可以尝试复原文章的封面图，也可以尝试修改色阶来改善图像的动态范围，把暗部和高光的细节都体现出来。

林林总总的，我把我整个瞎捣鼓的过程都记录了下来，希望给诸君一些工作流上的启发吧。

希望疫情☣️早点结束。

同时也祝愿所有人一帆风顺吧~⛵️