ASDF 示例

ASDF 示例#

用例： 从 FITS 文件创建 ASDF（高级科学数据格式）文件。

数据： COSMOS 区域的 CANDELS 图像。

工具： asdf, gwcs, astrocut。

跨仪器： 所有仪器。

文档： 本笔记本是 STScI 更大 后处理数据分析工具生态系统 的一部分。

介绍#

JWST 数据文件使用 高级科学数据格式。ASDF 元数据存储在 FITS 扩展中。JWST 管道软件从内存中的 数据模型 中读取和写入这些数据。

然而，读取和写入纯 ASDF 文件相对简单，可以完全跳过 FITS 和数据模型。此笔记本展示了使用 FITS 文件作为起点的一些 ASDF 方面。

导入#

  • 使用 astrocut 服务从 MAST 获取数据的 astrocut

  • 用于读取 FITS 文件的 astropy 的 fits 库

  • 用于处理天体坐标的 astropy 坐标 SkyCoord 对象

  • 用于绘制图形的 matplotlib

  • asdf 及 AsdfFile 对象

  • 用于以笔记本友好的方式查看头信息的 astropy Table 对象

  • 来自 modeling、coordinates 和 wcs 库的项目，用于将 FITS 关键字中的世界坐标系统信息转换为 gwcs 数据结构的示例。

from astrocut import fits_cut  # 导入astrocut库中的fits_cut模块，用于处理FITS文件

from astropy.io import fits  # 导入astropy.io中的fits模块，用于读取和写入FITS文件

from astropy.coordinates import SkyCoord  # 导入SkyCoord类，用于处理天文坐标

import matplotlib.pyplot as plt  # 导入matplotlib.pyplot模块，用于绘图

import asdf  # 导入asdf库，用于处理ASDF文件格式

from asdf import AsdfFile  # 从asdf库中导入AsdfFile类，用于操作ASDF文件

# For example 4

from astropy.table import Table  # 导入astropy.table中的Table类，用于处理表格数据

# For example 6

from astropy.modeling import models  # 导入astropy.modeling中的models模块，用于模型构建和拟合

from astropy import coordinates as coord  # 导入astropy.coordinates模块并命名为coord，用于处理坐标

from astropy import units as u  # 导入astropy.units模块并命名为u，用于处理物理单位

from astropy.wcs import WCS  # 导入astropy.wcs中的WCS类，用于处理世界坐标系统

from gwcs.wcstools import wcs_from_fiducial  # 从gwcs.wcstools导入wcs_from_fiducial函数，用于从基准点生成WCS

%matplotlib inline  # 在Jupyter Notebook中内嵌显示Matplotlib图形

获取数据#

我们将使用 astroquery 从 COSMOS 的 CANDELS 观测中获取一个切片。

# 定义数据源的URL
url = "https://archive.stsci.edu/pub/hlsp/candels/cosmos/cos-tot/v1.0/"

# 创建一个包含输入文件路径的列表
input_files = [url + "hlsp_candels_hst_acs_cos-tot-sect23_f606w_v1.0_drz.fits"]

# 创建一个天球坐标对象，表示中心坐标，单位为度
center_coord = SkyCoord("150.0946 2.38681", unit='deg')

# 定义切割区域的大小，单位为像素
cutout_size = [100, 100]

# 使用fits_cut函数从输入文件中提取切片图像
cutout_file = fits_cut(input_files, center_coord, cutout_size, single_outfile=True)  # 从输入文件中提取以中心坐标和切片大小为参数的切片图像，并将结果保存为单个输出文件

# 打印切割文件的名称或路径
print(cutout_file)

读取FITS文件，查看其结构并显示数据#

cutout_hdulist = fits.open(cutout_file)  # 打开指定的FITS文件并返回一个HDU列表

cutout_hdulist.info()  # 输出HDU列表的信息，包括每个HDU的类型、名称和数据维度

将FITS组件拆分，以便后续使用。

data = cutout_hdulist[1].data  # 从cutout_hdulist中提取第二个HDUL的数据信息

header0 = cutout_hdulist[0].header  # 从cutout_hdulist中提取第一个HDUL的头信息

header1 = cutout_hdulist[1].header  # 从cutout_hdulist中提取第二个HDUL的头信息

import matplotlib.pyplot as plt  # 导入绘图库

# 显示数据
plt.imshow(data)  # 使用imshow函数显示数据
plt.colorbar()  # 添加颜色条以显示数据值的范围
plt.title("JWST Data")  # 设置图像标题
plt.xlabel("X-axis")  # 设置X轴标签
plt.ylabel("Y-axis")  # 设置Y轴标签
plt.show()  # 显示图像

看起来您提供的内容不完整，似乎是一个标题或文件名。如果您有具体的JWST望远镜数据处理的Python代码需要我添加中文注释，请将代码粘贴在这里，我将很乐意为您添加行级中文注释。

示例 1：仅存储 ASDF 中的元数据键值对#

基本的 asdf 数据结构是一个字典（dictionary）。Astropy FITS 头对象（header object）像一个 Python 字典一样工作。我们可以将其复制到一个纯字典中，这在我们想要添加数据时会很有用。

tree1 = {**header1}  # 创建一个新的字典tree1，内容为header1的所有键值对

再加一行代码将其转换为 asdf。

myfile = AsdfFile(tree1)  # 创建一个AsdfFile对象，使用tree1作为输入数据

我们还不将其保存到文件中。首先，让我们检查一下树结构。

看起来您提到的“myfile.tree”似乎是一个文件名，但没有提供具体的代码或内容供我进行注释。如果您能提供相关的Python代码或数据处理的示例，我将能够为您添加中文注释并保持代码结构不变。

请您提供具体的代码或内容，我将很乐意帮助您！

示例 2：保存FITS头部注释#

FITS头部包含许多关键字的注释。虽然有点笨重，但我们可以通过在字典中存储一个元组而不仅仅是值来保存这些注释。同时，我们也可以丢弃一些不太有用的FITS关键字。

# 定义需要丢弃的关键字列表
toss_these = ['XTENSION', 'BITPIX', 'NAXIS', 'NAXIS1', 'NAXIS2', 'CHECKSUM',
              'DATASUM', 'EXTNAME', 'FILETYPE', 'PCOUNT', 'GCOUNT',
              'IRAF-BPX', 'IRAF-MIN', 'IRAF-MAX', 'IRAFNAME', 'IRAFTYPE']

# 初始化一个空字典，用于存储注释树
annotated_tree = {}

# 遍历header1中的每个卡片
for card in header1.cards:
    # 如果卡片的第一个元素不在丢弃列表中
    if card[0] not in toss_these:
        # 将卡片的第一个元素作为键，第二和第三个元素作为值存入字典
        annotated_tree[card[0]] = (card[1], card[2])

myfile = AsdfFile(annotated_tree)  # 创建一个ASDF文件对象，传入注释树

看起来您提到的“myfile.tree”可能是一个文件名，但没有提供具体的代码或内容。如果您能提供与JWST望远镜数据处理相关的Python代码，我将能够为您添加中文注释并保持代码结构不变。请将代码粘贴在这里，我会尽快帮助您！

值现在位于每个元组的第一个元素中。例如，要仅获取 CRVAL1 的值，我们可以这样做：

myfile['CRVAL1'][0]  # 访问myfile字典中'CRVAL1'键的第一个元素

# 更新 O_EXT_NM 键，因为它为空，导致在保存时崩溃

myfile.tree['O_EXT_NM'] = 'Original Name Filler'  # 将 O_EXT_NM 的值设置为 'Original Name Filler'

示例 3：将文件视为可搜索的表格#

对于 FITS 和 ASDF 文件，有时在长头文件中搜索可能会很痛苦。这个示例展示了一种将 ASDF 文件放入 Astropy 表格中的方法，然后使用 show_in_notebook 方法提供一个交互式的表格视图。您可以通过关键字搜索或通过点击列标题进行排序。

在这个示例中，我们将两个元素的元组拆分为值和注释。对于其他数据结构，我们只需将它们放入注释列中。

def tree_to_table(tree):
    # 获取树的所有键
    keys = list(tree.keys())

    # 初始化值和其他数据的列表
    values, other = ([] for i in range(2))

    # 遍历每个键
    for k in keys:
        try:
            # 尝试获取每个键的第一个值
            values += [tree[k][0]]
            # 尝试获取每个键的第二个值
            other += [tree[k][1]]
        except Exception as e:
            # 如果发生异常，添加None到values，并将整个tree[k]添加到other
            values += [None]
            other += [tree[k]]
            # 打印错误信息
            print("An error occured ", e)

    # 返回一个包含键、值和其他数据结构的表
    return Table([keys, values, other], names=['key', 'value', 'comment or data structure'])

t = tree_to_table(myfile.tree)  # 将树形结构转换为表格格式

t.show_in_notebook()  # 在笔记本中显示表格

示例 4：添加数据并写入文件#

这就是头部。现在我们只需将数据添加到字典中。我们可以使用任何我们喜欢的描述性键。也许我们应该称它为 data。

myfile['data'] = data  # 将数据赋值给myfile的'data'键，相当于myfile.tree['data'] = data

# 访问myfile字典中的'data'键，获取对应的值
myfile['data']

# 检查变量 myfile 的类型
type(myfile)  # 输出 myfile 的数据类型

# 显示 myfile 的信息，包括数据的基本统计和结构
myfile.info()  # 输出 myfile 的详细信息

myfile.write_to('myfile.asdf')  # 将数据写入名为'myfile.asdf'的文件

读取磁盘上的asdf文件并查看树形结构和数据

# 导入asdf模块以处理ASDF文件
import asdf

# 打开名为'myfile.asdf'的ASDF文件
ff = asdf.open('myfile.asdf')

# 打印文件树结构，显示文件内容的层次结构
ff.tree

plt.imshow(ff['data'])  # 使用imshow函数显示ff字典中'data'键对应的数据

示例 5：存储多个扩展#

敏锐的读者会注意到，前面的示例仅处理了 FITS 文件的扩展 1，未将主头信息包含在 ASDF 文件中。将 FITS 文件的多个扩展整理成 ASDF 文件没有规定的方式。一种选择是为每个扩展创建一个单独的字典，然后将这些字典组合成一个字典，例如：

ext1, ext2 = dict(**header0), dict(**header1)

tree = {'ext1':ext1, 'ext2':ext2}

在这种情况下，这样做有点愚蠢，因为 ASDF 文件中唯一可能有价值的信息可能是 ORIGIN, DATE, PROC_VER, RA_OBJ 和 DEC_OBJ。然而，查看扩展时，会发现那里有一个 ORIGIN，这将与主头中的 ORIGIN 冲突。（有趣的是，根据注释，它们的含义不同。）

看起来您提供的内容不完整，似乎是一个标题或文件名。如果您有具体的Python代码需要添加中文注释，请提供该代码，我将帮助您添加行级中文注释并保持代码结构不变。

解决命名冲突的一个方法可能是将这些额外信息放入它自己的命名空间，作为原始字典的一个子项。

keywords = ['ORIGIN', 'DATE', 'PROCVER', 'RA_OBJ', 'DEC_OBJ']  # 定义需要提取的关键字列表

primary_header = {}  # 初始化一个空字典，用于存储提取的主头信息

for card in header0.cards:  # 遍历头部中的每一张卡片

    if card[0] in keywords:  # 如果卡片的关键字在我们定义的关键字列表中

        primary_header[card[0]] = (card[1], card[2])  # 将关键字及其对应的值存入字典中

ff.tree['primary_header'] = primary_header  # 将提取的主头信息存入ff.tree中

ff.tree  # 返回ff.tree对象

示例 6：从 FITS WCS 关键字转换为 gwcs 对象#

这有点复杂，对于这个例子来说并不是特别有利。

通用世界坐标系统 gwcs 包旨在支持探测器坐标与世界坐标之间的复杂映射。在这种情况下，这有些过于复杂。但它确实展示了如何将复杂数据对象保存到 asdf 文件中。gwcs 包扩展了 asdf，以指定 wcs 对象。在此过程中，它广泛使用了在 asdf 标准 中定义的变换。

本示例遵循 gwcs 文档 中的说明。由于此示例中没有旋转，因此我们可以遵循使用便捷函数 wcs_from_fiducial 的 示例。

未来应该会有一些便捷方法，用于转换常见的 FITS WCS，并优雅地处理（不幸的是相当常见的）不一致性。

为了方便起见，我们将使用一个 astropy 的 WCS 对象来获取世界坐标系统（WCS）信息。这个特定的例子揭示了在 FITS 文件中非常常见的一个问题：冗余且可能不一致的 WCS 信息。在这种情况下，文件同时使用了 PC 矩阵方法和 CD 矩阵方法。我们将把 PC 矩阵视为“真相”。

fitswcs = WCS(header1)  # 使用header1创建WCS对象

print(fitswcs)  # 打印WCS对象以查看其内容

从WCS（世界坐标系统）中获取一些值，以提高下方单元格的可读性。

# 从fitswcs.wcs中提取CRVAL（坐标参考值）
crval1, crval2 = fitswcs.wcs.crval

# 从fitswcs.wcs中提取单位并转换为单位对象
cunit1, cunit2 = [u.Unit(cu) for cu in fitswcs.wcs.cunit]

# 从fitswcs.wcs中提取PC矩阵（像素坐标变换矩阵）
pcmatrix = fitswcs.wcs.pc

# 输出单位1和单位2
cunit1, cunit2

要从天空中的一个指向创建一个世界坐标系统（WCS），至少需要将天空坐标和投影传递给该函数。

# 创建一个SkyCoord对象，使用给定的天球坐标值和单位，框架为ICRS
fiducial = coord.SkyCoord(crval1*cunit1, crval2*cunit2, frame='icrs')

# 创建一个切片到天空的模型，使用TAN投影
tan = models.Pix2Sky_TAN()

创建一个坐标变换的流程。在这种情况下，首先应用平移（shift），然后进行重缩放（rescaling）。函数 wcs_from_fiducial 将这些操作添加到天空投影（sky projection）之前。

trans = models.Shift(-crval1) & models.Shift(-crval2) |\  # 创建一个平移模型，分别对crval1和crval2进行负向平移
        models.Scale(-pcmatrix[0, 0]) & models.Scale(pcmatrix[1, 1])  # 创建一个缩放模型，分别对pcmatrix的元素进行负向和正向缩放

# 从给定的fiducial点生成WCS（世界坐标系统）对象，使用正切投影和指定的变换
wcsobj = wcs_from_fiducial(fiducial, projection=tan, transform=trans)

# 输出生成的WCS对象
wcsobj

从ASDF头部中删除现已过时的FITS WCS关键字。它简洁得令人愉悦，可以说几乎没有任何有用的信息。

fits_wcs_keywords = [
    'CRPIX1', 'CRVAL1', 'CTYPE1', 'CD1_1', 'CD2_1', 'CRPIX2', 'CRVAL2',  # WCS相关的关键字列表
    'CTYPE2', 'CD1_2', 'CD2_2', 'WCSAXES', 'PC1_1', 'PC2_2', 'CDELT1',  # 包含坐标参考点和单位等信息
    'CDELT2', 'CUNIT1', 'CUNIT2', 'LONPOLE', 'LATPOLE', 'RADESYS']      # 用于描述坐标系的其他参数

# 从树结构中移除指定的WCS关键字，如果不存在则返回None
[ff.tree.pop(old_kw, None) for old_kw in fits_wcs_keywords]

将wcs对象添加到树中。gwcs包负责将该对象序列化为ASDF的相关机制。

ff.tree['wcs'] = wcsobj  # 将WCS对象赋值给ff.tree字典中的'wcs'键

t = tree_to_table(ff.tree)  # 将树结构转换为表格格式

t.show_in_notebook()  # 在笔记本中显示表格
Contents