MRS立方体分析

MRS立方体分析#

用例： 从IFU立方体中提取空间-光谱特征并测量其属性。

数据： 模拟的MIRI MRS AGB星光谱。

工具： specutils, astropy, scipy。

跨仪器： NIRSpec, MIRI。

文档： 本笔记本是STScI更大后处理数据分析工具生态系统的一部分，可以直接从JDAT Notebook GitHub目录中下载。

模拟来源： MIRISim

管道版本： JWST管道

注意： 本笔记本包含使用MIRISim（https://wiki.miricle.org//bin/view/Public/MIRISim_Public）获得的MIRI模拟数据立方体，并通过JWST管道（https://jwst-pipeline.readthedocs.io/en/latest/）处理，代表晚M型星的点源光谱。

介绍#

本笔记本分析了一颗星，其尘埃SED对应于Kraemer等人（2002）和Sloan等人（2003）的ISO SWS光谱W Per，以覆盖MRS光谱范围5-28微米。JWST光谱立方体的分析需要提取感兴趣的空间-光谱特征并测量其属性。

这是第二个笔记本，将使用specutils进行数据分析。具体而言，它将拟合一个模型光球/黑体到光谱中。然后，它将计算每个尘埃和分子特征的质心、线积分通量和当量宽度。

本笔记本利用在第一个笔记本（JWST_Mstar_dataAnalysis_runpipeline.ipynb）中创建的减缩数据，尽管您不必运行该笔记本即可使用本笔记本。第一个笔记本中创建的所有数据都可以在本笔记本中下载。

待办事项：#

  • 制作函数以使用不同的提取方式从数据立方体中提取光谱。

  • 用恒星光球模型替换光谱中光球部分的黑体拟合。

  • 确保在质心、线积分通量和当量宽度的计算中正确传播了误差。（科学作者）

  • 在specutils框架内制作简单函数，以拟合连续体并测量宽固态和分子特征的质心、线积分通量和当量宽度。

导入包#

# 导入有用的python包
import numpy as np  # 导入numpy库，用于数值计算

# 导入在笔记本中内联显示图像的包
import matplotlib.pyplot as plt  # 导入matplotlib库中的pyplot模块，用于绘图

%matplotlib inline  # 设置Jupyter Notebook以内联方式显示图像

# 设置一般绘图选项
params = {
    'legend.fontsize': '18',  # 图例字体大小
    'axes.labelsize': '18',    # 坐标轴标签字体大小
    'axes.titlesize': '18',    # 坐标轴标题字体大小
    'xtick.labelsize': '18',    # x轴刻度标签字体大小
    'ytick.labelsize': '18',    # y轴刻度标签字体大小
    'lines.linewidth': 2,       # 线条宽度
    'axes.linewidth': 2,        # 坐标轴线宽度
    'animation.html': 'html5'   # 动画输出格式设置为html5
}

plt.rcParams.update(params)  # 更新matplotlib的绘图参数
plt.rcParams.update({'figure.max_open_warning': 0})  # 禁用最大打开图形警告

# 导入astropy包

from astropy import units as u  # 导入单位模块
from astropy.io import ascii  # 导入ASCII输入输出模块
from astropy.nddata import StdDevUncertainty  # 导入标准偏差不确定性模块
from astropy.utils.data import download_file  # 导入下载文件的工具

# 处理一维光谱

from specutils import Spectrum1D  # 导入一维光谱类
from specutils.manipulation import box_smooth, extract_region  # 导入平滑和提取区域的函数
from specutils.analysis import line_flux, centroid, equivalent_width  # 导入分析光谱的函数
from specutils.spectra import SpectralRegion  # 导入光谱区域类
from specutils import SpectrumList  # 导入光谱列表类

from jdaviz import Specviz  # 导入Specviz可视化工具
from jdaviz import Cubeviz  # 导入Cubeviz可视化工具

# 显示视频

from IPython.display import YouTubeVideo  # 导入显示YouTube视频的工具

# 使用Pickle保存和加载对象

import pickle  # 导入pickle模块，用于对象的序列化和反序列化

def save_obj(obj, name):
    # 定义保存对象的函数，接受对象和文件名作为参数
    with open(name, 'wb') as f:  # 以二进制写入模式打开文件
        pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)  # 使用pickle将对象序列化并写入文件

def load_obj(name):
    # 定义加载对象的函数，接受文件名作为参数
    with open(name, 'rb') as f:  # 以二进制读取模式打开文件
        return pickle.load(f)  # 使用pickle从文件中反序列化并返回对象

def checkKey(dict, key):  # 定义一个函数，检查字典中是否存在指定的键

    if key in dict.keys():  # 如果键在字典的键中
        print("Present, ", end=" ")  # 打印“存在”，不换行
        print("value =", dict[key])  # 打印该键对应的值
        return True  # 返回True，表示键存在
    else:  # 如果键不在字典中
        print("Not present")  # 打印“不存在”
        return False  # 返回False，表示键不存在

import warnings  # 导入警告模块

warnings.simplefilter('ignore')  # 忽略所有警告

# 检查是否存在Pipeline 3处理后的数据，如果不存在，则下载数据

import os  # 导入操作系统模块

import urllib.request  # 导入用于请求URL的模块

import tarfile  # 导入处理tar文件的模块

# 检查文件是否存在
if os.path.exists("combine_dithers_all_exposures_ch1-long_s3d.fits"):
    print("Pipeline 3 Data Exists")  # 如果文件存在，打印信息
else:
    # 数据下载链接
    url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/reduced.tar.gz'
    
    # 从URL下载数据并保存为reduced.tar.gz
    urllib.request.urlretrieve(url, './reduced.tar.gz')

    base_extract_to = os.path.abspath(".")  # 获取当前目录的绝对路径

    # 打开并安全地提取tar归档中的文件
    with tarfile.open('./reduced.tar.gz', "r:gz") as tar:
        for member in tar.getmembers():  # 遍历tar文件中的每个成员
            # 计算文件提取后的绝对路径
            member_path = os.path.abspath(os.path.join(base_extract_to, member.name))

            # 检查文件路径是否在基础提取目录内
            if member_path.startswith(base_extract_to):
                # 仅提取安全的文件，直接到基础目录
                tar.extract(member, path=base_extract_to)
            else:
                print(f"Skipped {member.name} due to potential security risk")  # 如果路径不安全，打印跳过信息

    # 移动提取的fits文件到当前目录
    os.system('mv reduced/*fits .')

从笔记本1中读取12个光谱，使用 SpectrumList 并创建一个主光谱。

MIRI MRS（中红外仪器中分辨率光谱）数据通常会为每个通道（Channel）和子带（Sub-Band）关联一个单独的1D光谱和一个3D立方体。

# 读取所有的X1D光谱到一个光谱列表，并合并成一个单一的Spectrum1D对象

ddd = '.'  # 设置当前目录为数据读取的路径

splist = SpectrumList.read(ddd)  # 从指定路径读取光谱列表

wlallorig = []  # 初始化原始波长列表
fnuallorig = []  # 初始化原始光度列表
dfnuallorig = []  # 初始化原始不确定度列表

# 遍历光谱列表中的每个光谱
for bndind in range(len(splist)):
    # 遍历当前光谱的每个波长点
    for wlind in range(len(splist[bndind].spectral_axis)):
        wlallorig.append(splist[bndind].spectral_axis[wlind].value)  # 添加波长值
        fnuallorig.append(splist[bndind].flux[wlind].value)  # 添加光度值
        dfnuallorig.append(splist[bndind].uncertainty[wlind].array)  # 添加不确定度值

# 将原始列表转换为NumPy数组
wlallarr = np.array(wlallorig)  # 波长数组
fnuallarr = np.array(fnuallorig)  # 光度数组
dfnuallarr = np.array(dfnuallorig)  # 不确定度数组

srtind = np.argsort(wlallarr)  # 获取波长数组的排序索引

# 根据排序索引重新排列波长和光度
wlall = wlallarr[srtind]  # 排序后的波长
fnuall = fnuallarr[srtind]  # 排序后的光度

# 开发者注释：我们设置了0.0001的误差，但不确定度需要修正
dfnuall = (0.0001) * np.ones(len(fnuall))  # 初始化不确定度为0.0001

Specviz 对 SpectrumList 的可视化#

您可以在 Jupyter notebook 中使用 Specviz 可视化光谱列表。

视频 1：#

此 Specviz 教学演示来自 STScI 的官方 YouTube 频道，提供了对 Specviz 的介绍。

vid = YouTubeVideo("zLyRnfG32Bo")  # 创建一个YouTubeVideo对象，传入视频ID

display(vid)  # 显示该YouTube视频

读取 SpectrumList（12 个独特光谱）#

# 打开Specviz以显示光谱数据
specviz = Specviz()  # 创建Specviz对象

specviz.show()  # 显示Specviz界面

从上面加载光谱列表。请注意，列表中的第一个光谱会自动显示。您需要在“DATA”下拉菜单中启用其余光谱，然后点击工具栏中的“Home”按钮，并相应地缩放我们的图表，以查看其他光谱。#

# 从上面加载光谱列表
specviz.load_data(splist)  # 加载光谱数据到SpecViz中

Cubeviz 可视化#

您还可以使用 Cubeviz 在 Jupyter notebook 中可视化图像。

视频 2：#

此 Cubeviz 教学演示来自 STScI 的官方 YouTube 频道，提供了对 Cubeviz 的介绍。

vid = YouTubeVideo("HMSYwiH3Gl4")  # 创建一个YouTubeVideo对象，传入视频ID

display(vid)  # 显示该视频对象

cubeviz = Cubeviz()  # 创建一个Cubeviz实例，用于可视化立方体数据

cubeviz.show()  # 显示Cubeviz界面

开发者注释。需要选择不同于米的单位。 https://jira.stsci.edu/browse/JDAT-1792#

# 在这里，我们将数据加载到Cubeviz应用程序中以进行可视化检查。

# 在这种情况下，我们只查看一个通道，因为与Specviz不同，Cubeviz一次只能加载一个数据立方体。

ch1short_cubefile = 'combine_dithers_all_exposures_ch1-long_s3d.fits'  # 定义通道1的文件名

cubeviz.load_data(ch1short_cubefile)  # 加载数据立方体到Cubeviz中

接下来，您需要定义一个特定于AGB星的像素区域子集。您可以使用区域工具按钮，在AGB星周围绘制一个圆形区域，约在像素 x=20, y=30 的位置。

# 现在从光谱查看器中提取光谱

# 需要展示如何使用光谱提取插件并计算平均光谱而不是总和光谱

try:
    # 尝试获取仅包含AGB星的光谱数据（子集1，总和）
    spec_agb = cubeviz.get_data('Spectrum (Subset 1, sum)')  # AGB星仅
    print(spec_agb)  # 打印AGB星的光谱数据
    spec_agb_exists = True  # 标记AGB星光谱数据存在

except Exception:
    # 如果没有选择子集，则捕获异常
    print("There are no subsets selected.")  # 打印错误信息
    spec_agb_exists = False  # 标记AGB星光谱数据不存在
    # 获取整个视场的光谱数据（总和）
    spec_agb = cubeviz.get_data('Spectrum (sum)')  # 整个视场
    print(spec_agb)  # 打印整个视场的光谱数据

开发者说明：由于 Cubeviz 仅能同时显示一个数据立方体，因此您无法在当前时间提取完整的光谱。因此，您应该使用上面定义的光谱（’spec’）。#

创建一维光谱对象#

wav = wlall * u.micron  # 波长：微米

fl = fnuall * u.Jy  # Fnu：杰伊

efl = dfnuall * u.Jy  # 误差通量：杰伊

# 创建一个一维光谱对象

spec = Spectrum1D(spectral_axis=wav, flux=fl,

                  uncertainty=StdDevUncertainty(efl))  # 不确定性使用标准偏差不确定性

# 对光谱应用5像素的箱形平滑

spec_bsmooth = box_smooth(spec, width=5)  # 使用box_smooth函数对光谱进行平滑处理，宽度为5像素

# 绘制光谱和光谱平滑后的图像以检查特征

plt.figure(figsize=(8, 4))  # 创建一个8x4英寸的图形

plt.plot(spec.spectral_axis, spec.flux, label='Source')  # 绘制原始光谱

plt.plot(spec.spectral_axis, spec_bsmooth.flux, label='Smoothed')  # 绘制平滑后的光谱

plt.xlabel('Wavelength (microns)')  # 设置x轴标签为波长（微米）

plt.ylabel("Flux ({:latex})".format(spec.flux.unit))  # 设置y轴标签为通量，单位为光谱单位

plt.ylim(-0.05, 0.15)  # 设置y轴的范围

# 叠加原始输入光谱以进行比较

origspecfile = fn = download_file(  # 下载原始光谱数据文件
    'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/63702662.txt', cache=True)

origdata = ascii.read(origspecfile)  # 读取下载的光谱数据

wlorig = origdata['col1']  # 获取原始光谱的波长数据

# 原始数据以mJy为单位，转换为Jy以便与管道输出进行比较

fnujyorig = origdata['col2'] * 0.001  # 将原始光谱的通量从mJy转换为Jy

plt.plot(wlorig, fnujyorig, '.', color='grey',  # 绘制原始输入光谱，使用灰色点表示
         markersize=1, label='Original Input')

plt.legend(frameon=False, fontsize='medium')  # 添加图例，不带边框，字体大小为中等

plt.tight_layout()  # 自动调整子图参数以给图形留出足够的空间

plt.show()  # 显示图形

plt.close()  # 关闭当前图形

可视化分析上述从12个单独光谱创建的单光谱1D对象#

您可以在Specviz中可视化提取的光谱。

# 打开Specviz以显示光谱数据
specviz = Specviz()  # 创建Specviz实例

specviz.show()  # 显示Specviz界面

# 加载AGB星体的光谱数据
specviz.load_data(spec_agb)

# 导入必要的库
import numpy as np  # 导入NumPy库用于数值计算
import matplotlib.pyplot as plt  # 导入Matplotlib库用于绘图
from specviz import Spectrum  # 从specviz库导入Spectrum类

# 创建一个新的光谱对象
spectrum = Spectrum()  # 实例化Spectrum对象

# 生成模拟光谱数据
wavelength = np.linspace(4000, 7000, 1000)  # 创建波长范围从4000到7000的数组
flux = np.random.normal(size=wavelength.size)  # 生成随机的光谱强度数据

# 将数据添加到光谱对象中
spectrum.add_data(wavelength, flux)  # 将波长和强度数据添加到光谱对象

# 平滑光谱数据
smoothed_flux = np.convolve(flux, np.ones(10)/10, mode='same')  # 使用卷积平滑光谱数据

# 绘制原始和光滑后的光谱
plt.figure(figsize=(10, 5))  # 设置绘图窗口大小
plt.plot(wavelength, flux, label='原始光谱', alpha=0.5)  # 绘制原始光谱
plt.plot(wavelength, smoothed_flux, label='平滑光谱', color='red')  # 绘制平滑后的光谱
plt.xlabel('波长 (Å)')  # 设置x轴标签
plt.ylabel('强度')  # 设置y轴标签
plt.title('光谱平滑示例')  # 设置图表标题
plt.legend()  # 显示图例
plt.show()  # 显示绘图

数据分析#

分析从上述 spectrumlist 列表创建的 spectrum1d 对象（spec）。

拟合连续谱 - 找到最佳拟合模板（恒星光球模型或黑体）#

对于具有光球成分的AGB星，拟合恒星光球模型或黑体到光谱的短波长端。

开发者备注：理想情况下希望使用一组Phoenix模型来拟合光球。#

我认为 template_comparison 可能是一个不错的函数，可以与已设置为与 pysynphot 接口的Phoenix模型一起使用。

目前暂时切换到黑体。

# 导入必要的库
import numpy as np  # 导入NumPy库用于数值计算
import matplotlib.pyplot as plt  # 导入Matplotlib库用于绘图
from astropy.modeling import models, fitting  # 导入Astropy库中的模型和拟合工具

# 创建一个黑体模型函数
def blackbody_model(wavelengths, temperature):
    # 使用Planck公式计算黑体辐射
    return (2.0 * 6.626e-34 * 3.0e8) / (wavelengths**5 * (np.exp((6.626e-34 * 3.0e8) / (wavelengths * 1.38e-23 * temperature)) - 1))

# 生成示例数据
wavelengths = np.linspace(1e-7, 3e-6, 100)  # 生成从100nm到3000nm的波长数据
true_temperature = 5800  # 设定真实温度为5800K
spectrum = blackbody_model(wavelengths, true_temperature)  # 计算真实的黑体光谱

# 添加一些噪声以模拟真实观测数据
noise = np.random.normal(0, 0.05, spectrum.shape)  # 生成高斯噪声
observed_spectrum = spectrum + noise  # 将噪声添加到光谱中

# 绘制观测光谱
plt.figure(figsize=(10, 6))  # 设置绘图大小
plt.plot(wavelengths, observed_spectrum, label='Observed Spectrum', color='blue')  # 绘制观测光谱
plt.xlabel('Wavelength (m)')  # 设置x轴标签
plt.ylabel('Intensity')  # 设置y轴标签
plt.title('Observed Spectrum with Noise')  # 设置图表标题
plt.legend()  # 显示图例
plt.show()  # 显示绘制的图形

# 拟合黑体模型
# 创建一个黑体模型实例
bb_model = models.BlackBody(temperature=5000)  # 初始温度设为5000K
fitter = fitting.LevMarLSQFitter()  # 使用Levenberg-Marquardt拟合器

# 拟合观测光谱
fitted_model = fitter(bb_model, wavelengths, observed_spectrum)  # 拟合模型到观测数据

# 绘制拟合结果
plt.figure(figsize=(10, 6))  # 设置绘图大小
plt.plot(wavelengths, observed_spectrum, label='Observed Spectrum', color='blue')  # 绘制观测光谱
plt.plot(wavelengths, fitted_model(wavelengths), label='Fitted Blackbody Model', color='red')  # 绘制拟合的黑体模型
plt.xlabel('Wavelength (m)')  # 设置x轴标签
plt.ylabel('Intensity')  # 设置y轴标签
plt.title('Fitted Blackbody Model to Observed Spectrum')  # 设置图表标题
plt.legend()  # 显示图例
plt.show()  # 显示绘制的图形

spectra = specviz.get_spectra()  # 从Specviz获取光谱数据

a = checkKey(spectra, "BB1")  # 检查光谱数据中是否存在"BB1"键

if a is True:  # 如果存在"BB1"键

    # 从Specviz中提取黑体拟合数据
    blackbody = spectra["BB1"]  # 将"BB1"光谱数据赋值给blackbody变量

else:  # 如果不存在"BB1"键
    print("No Blackbody")  # 输出提示信息

    # 下载黑体光谱数据文件
    fn = download_file(
        'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/blackbody.fits', cache=True)

    blackbody = Spectrum1D.read(fn)  # 从下载的文件中读取光谱数据并赋值给blackbody变量

# 删除当前目录中任何现有的输出文件

if os.path.exists("blackbody.fits"):  # 检查文件是否存在

    os.remove("blackbody.fits")  # 如果存在，则删除该文件

else:  # 如果文件不存在
    print("The blackbody.fits file does not exist")  # 输出提示信息

# 如果需要，可以保存文件。否则保持注释状态。

# blackbody.write('blackbody.fits')  # 将黑体数据写入名为 'blackbody.fits' 的文件

# 将 blackbody.flux 重命名为 ybest
ybest = blackbody.flux  # 将黑体辐射的 flux 属性赋值给 ybest

开发者说明：此时，需要使用为MRS编写的特殊函数将提取的12个1D光谱拼接在一起。#

# 绘制光谱及其在短波长区域的模型拟合

plt.figure(figsize=(8, 4))  # 创建一个8x4英寸的图形

plt.plot(spec.spectral_axis, spec.flux, label='Source')  # 绘制源光谱

plt.plot(spec.spectral_axis, ybest, label='BB')  # 绘制黑体模型拟合

plt.xlabel('Wavelength (microns)')  # 设置x轴标签为波长（微米）

plt.ylabel("Flux ({:latex})".format(spec.flux.unit))  # 设置y轴标签为通量，并格式化单位

plt.title("Spectrum with blackbody fit")  # 设置图形标题

plt.legend(frameon=False, fontsize='medium')  # 添加图例，不带边框，字体大小为中等

plt.tight_layout()  # 自动调整子图参数以给图形留出足够的空间

plt.ylim(-0.05, 0.15)  # 设置y轴范围

plt.show()  # 显示图形

plt.close()  # 关闭当前图形

# 现在减去黑体光谱并绘制底层尘埃连续谱

plt.figure(figsize=(8, 4))  # 创建一个8x4英寸的图形

plt.plot(spec.spectral_axis, spec.flux.value -  # 绘制尘埃光谱
         ybest.value, color='purple', label='Dust spectra')  # 使用紫色绘制尘埃光谱

plt.axhline(0, color='r', linestyle='dashdot', alpha=0.5)  # 绘制y=0的水平线，红色，虚线，透明度为0.5

plt.xlabel('Wavelength (microns)')  # 设置x轴标签为波长（微米）

plt.ylabel("Flux ({:latex})".format(spec.flux.unit))  # 设置y轴标签为通量，并格式化单位

plt.title("Continuum-subtracted spectrum")  # 设置图形标题

plt.legend(frameon=False, fontsize='medium')  # 添加图例，不带边框，字体大小为中等

plt.tight_layout()  # 自动调整子图参数以给图形留出足够的空间

plt.ylim(-0.05, 0.15)  # 设置y轴范围

plt.show()  # 显示图形

plt.close()  # 关闭当前图形

现在我们有了尘埃连续谱，想要寻找特征并测量它们的属性。#

想要找到：

  • 等效宽度 (Equivalent width)

  • 等效通量 (Equivalent flux)

  • 光学深度 (Optical depth)

  • 重心 (Centroids) = 在两侧通量各占一半的波长

作为一个例子，让我们关注无定形硅酸盐的10微米区域。#

方法 - 重复使用

  • 拟合一个样条曲线 (spline) 到减去光球连续谱的光谱中，排除该拟合中的特征。

  • 将光谱裁剪到该波长区域，因为样条曲线的大小现在与光谱的完整波长范围不同。

  • 制作一个减去连续谱的光谱和一个归一化的连续谱。

  • 将通量的单位从焦耳每秒每平方米每波长 (Jy) 转换为瓦特每平方米每波长 (W/m²/wavelength)，以便在进行线积分后使用更好的单位。

  • 确定特征线通量，单位为瓦特每平方米 (W/m²) 和特征重心。使用减去连续谱的光谱。

  • 确定特征的等效宽度。使用归一化的连续谱。

  • 确保误差已正确传播。

  • 将这些结果存储在一个表格中。

  • 光谱中通常存在多个分子和尘埃特征。对每个特征重复上述步骤。

注意

这似乎是一种冗长的方法。有没有更简单的方式？

例如，一个工具可以接受四个波长，使用从 lam0 到 lam1 和 lam2 到 lam3 的数据拟合一条线，然后

将减去连续谱的光谱用于从 lam1 到 lam2 的线积分，并进行误差传播，这对于尘埃特征需要多次进行。

但是在当前的 spectra1d 框架下，这需要手动编写许多步骤，并且在处理两个特征后就变得非常繁琐，更不用说20个以上的特征了。对于带有不确定性的综合线重心和提取的等效宽度，也需要类似的框架。

# 从光谱中减去连续背景，并在新的specviz实例中绘制

bbsub_spectra = spec - ybest.value     # 减去连续背景后的光谱 - 仅包含尘埃

specviz = Specviz()  # 创建一个Specviz对象，用于光谱数据的可视化

specviz.show()  # 显示Specviz界面

# 加载数据到specviz中
specviz.load_data(bbsub_spectra)  # 将bbsub_spectra数据加载到specviz工具中

# 导入必要的库
import numpy as np  # 用于数值计算
import matplotlib.pyplot as plt  # 用于绘图
from scipy.optimize import curve_fit  # 用于曲线拟合

# 定义多项式函数
def polynomial(x, *coeffs):  # 定义多项式函数，coeffs为多项式系数
    return sum(c * x**i for i, c in enumerate(coeffs))  # 计算多项式值

# 生成示例数据
x_data = np.linspace(0, 10, 100)  # 生成0到10之间的100个点
y_data = polynomial(x_data, 1, -2, 3) + np.random.normal(0, 1, x_data.size)  # 生成带噪声的多项式数据

# 定义两个光谱区域的索引
region1_indices = (x_data < 5)  # 第一区域：x小于5
region2_indices = (x_data >= 5)  # 第二区域：x大于等于5

# 拟合第一区域
popt_region1, pcov_region1 = curve_fit(polynomial, x_data[region1_indices], y_data[region1_indices], p0=[1, -1, 1])  # 拟合第一区域的多项式
fit_region1 = polynomial(x_data[region1_indices], *popt_region1)  # 计算拟合值

# 拟合第二区域
popt_region2, pcov_region2 = curve_fit(polynomial, x_data[region2_indices], y_data[region2_indices], p0=[1, -1, 1])  # 拟合第二区域的多项式
fit_region2 = polynomial(x_data[region2_indices], *popt_region2)  # 计算拟合值

# 绘制结果
plt.figure(figsize=(10, 6))  # 设置图形大小
plt.scatter(x_data, y_data, label='Data', color='gray', alpha=0.5)  # 绘制原始数据点
plt.plot(x_data[region1_indices], fit_region1, label='Fit Region 1', color='blue')  # 绘制第一区域的拟合曲线
plt.plot(x_data[region2_indices], fit_region2, label='Fit Region 2', color='red')  # 绘制第二区域的拟合曲线
plt.xlabel('X-axis')  # 设置X轴标签
plt.ylabel('Y-axis')  # 设置Y轴标签
plt.title('Polynomial Fit to Two Spectral Regions')  # 设置图形标题
plt.legend()  # 显示图例
plt.show()  # 显示图形

spectra = specviz.get_spectra()  # 从Specviz获取光谱数据

a = checkKey(spectra, "PolyFit")  # 检查光谱数据中是否包含"PolyFit"键

if a is True:  # 如果包含"PolyFit"键

    # 从Specviz中提取多项式拟合数据
    poly = spectra["PolyFit"]  # 将多项式拟合数据赋值给poly

else:  # 如果不包含"PolyFit"键
    print("No Polyfit")  # 输出提示信息

    fn = download_file(  # 下载多项式拟合文件
        'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/poly.fits', cache=True)

    poly = Spectrum1D.read(fn)  # 读取下载的FITS文件并赋值给poly

# 删除当前目录中任何现有的输出文件

if os.path.exists("poly.fits"):  # 检查文件是否存在

    os.remove("poly.fits")  # 如果存在，则删除该文件

else:  # 如果文件不存在

    print("The poly.fits file does not exist")  # 输出提示信息

# 如果需要，可以保存文件。否则保持注释状态。

# poly.write('poly.fits')  # 将多项式数据写入名为 'poly.fits' 的文件

# 在8.0 - 8.1和14.9 - 15.0微米之间拟合局部连续谱
# （即排除该谱线本身）

sw_region = 8.0   # lam0，局部连续谱的起始波长
sw_line = 8.1     # lam1，谱线的起始波长
lw_line = 14.9    # lam2，谱线的结束波长
lw_region = 15.0  # lam3，局部连续谱的结束波长

# 缩放到谱线复合区域并提取
line_reg_10 = SpectralRegion([(sw_region*u.um, lw_region*u.um)])  # 定义谱线区域
line_spec = extract_region(bbsub_spectra, line_reg_10)  # 从背景减去的光谱中提取谱线区域
polysub = extract_region(poly, line_reg_10)  # 从多项式中提取谱线区域
line_y_continuum = polysub.flux  # 获取局部连续谱的flux

# -----------------------------------------------------------------

# 生成去除连续谱和归一化的光谱
line_spec_norm = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux /
                            line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))  # 归一化光谱
line_spec_consub = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux -
                              line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))  # 去除连续谱的光谱

# -----------------------------------------------------------------

# 绘制尘埃特征和局部连续谱拟合区域
plt.figure(figsize=(8, 4))  # 设置图形大小

plt.plot(line_spec.spectral_axis, line_spec.flux.value,
         label='Dust spectra 10 micron region')  # 绘制尘埃光谱

plt.plot(line_spec.spectral_axis, line_y_continuum, label='Local continuum')  # 绘制局部连续谱

plt.xlabel('Wavelength (microns)')  # x轴标签
plt.ylabel("Flux ({:latex})".format(spec.flux.unit))  # y轴标签
plt.title(r"10$\mu$m feature plus local continuum")  # 图形标题
plt.legend(frameon=False, fontsize='medium')  # 图例设置
plt.tight_layout()  # 自动调整子图参数
plt.show()  # 显示图形
plt.close()  # 关闭图形

# -----------------------------------------------------------------

# 绘制去除连续谱的10微米特征
plt.figure(figsize=(8, 4))  # 设置图形大小

plt.plot(line_spec.spectral_axis, line_spec_consub.flux, color='green',
         label='continuum subtracted')  # 绘制去除连续谱的光谱

plt.xlabel('Wavelength (microns)')  # x轴标签
plt.ylabel("Flux ({:latex})".format(spec.flux.unit))  # y轴标签
plt.title(r"Continuum subtracted 10$\mu$m feature")  # 图形标题
plt.tight_layout()  # 自动调整子图参数
plt.show()  # 显示图形
plt.close()  # 关闭图形

# 导入Specviz库并创建一个Specviz实例
specviz = Specviz()

# 显示Specviz界面
specviz.show()

# 加载经过连续体减法处理的光谱数据，并为其指定标签
specviz.load_data(line_spec_consub, data_label='Continuum Subtraction')

# 加载经过归一化处理的光谱数据，并为其指定标签
specviz.load_data(line_spec_norm, data_label='Normalized')

# 导入必要的库
import numpy as np  # 导入NumPy库用于数值计算
import matplotlib.pyplot as plt  # 导入Matplotlib库用于绘图
from specviz import Specviz  # 从specviz模块导入Specviz类

# 创建Specviz实例
specviz = Specviz()  # 实例化Specviz对象

# 加载数据
data = np.load('spectrum_data.npy')  # 从.npy文件中加载光谱数据
specviz.load_data(data)  # 将数据加载到Specviz中

# 绘制光谱
plt.figure()  # 创建一个新的图形
plt.plot(data)  # 绘制光谱数据
plt.title('Spectrum')  # 设置图形标题
plt.xlabel('Wavelength')  # 设置x轴标签
plt.ylabel('Flux')  # 设置y轴标签
plt.show()  # 显示图形

# 测量线
line_positions = [5000, 6000, 7000]  # 定义需要测量的线的位置（单位：Å）
for position in line_positions:  # 遍历每个线的位置
    plt.axvline(x=position, color='r', linestyle='--')  # 在指定位置绘制垂直线

# 显示测量结果
plt.title('Measured Lines in Spectrum')  # 设置图形标题
plt.show()  # 显示图形

# 在笔记本中分析10微米线的替代方法。计算线通量；线心；等效宽度

line_centroid = centroid(  # 计算线心
    line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))  # 使用给定的光谱区域

line_flux_val = line_flux(  # 计算线通量
    line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))  # 使用给定的光谱区域

equivalent_width_val = equivalent_width(line_spec_norm)  # 计算等效宽度

# 黑客方法将线通量值转换为更常规的单位

# 由于光谱具有混合单位：f_nu+lambda，因此这是必要的
line_flux_val = (line_flux_val * u.micron).to(u.W * u.m**-2 * u.micron,  # 转换单位
                                              u.spectral_density(line_centroid)) / u.micron  # 使用线心进行光谱密度转换

print("Line_centroid: {:.6} ".format(line_centroid))  # 打印线心
print("Integrated line_flux: {:.6} ".format(line_flux_val))  # 打印集成线通量
print("Equivalent width: {:.6} ".format(equivalent_width_val))  # 打印等效宽度

# 计算10微米特征的光学深度

# 通过对比线谱的通量与连续谱的通量计算光学深度
tau = -(np.log(line_spec.flux.value / line_y_continuum.value))

# 创建一个新的光谱对象，包含光谱轴和计算得到的光学深度
optdepth_spec = Spectrum1D(spectral_axis=line_spec.spectral_axis,

                           flux=tau*(u.Jy/u.Jy))  # 将光学深度的单位设置为Jy

# 绘制10微米区域的光学深度与波长的关系图

plt.figure(figsize=(10, 6))  # 创建一个10x6英寸的图形

plt.plot(optdepth_spec.spectral_axis, optdepth_spec.flux)  # 绘制光学深度谱的波长与光学深度的关系

plt.xlabel("Wavelength ({:latex})".format(spec.spectral_axis.unit))  # 设置x轴标签为波长，并格式化单位

plt.ylabel('Tau')  # 设置y轴标签为光学深度（Tau）

plt.tight_layout()  # 自动调整子图参数，使之填充整个图像区域

plt.show()  # 显示绘制的图形

plt.close()  # 关闭当前图形

注意 此时请重复 所有 上述步骤，以隔离固态特征，例如在大约 13.3 微米处的橄榄石特征。

这是展示使用Cubeviz和Specviz对MIRI MRS光谱进行一些基本分析的笔记本的结尾。还有更多的分析是可能的。

额外资源#

关于此笔记本#

作者： Olivia Jones，项目科学家，英国天文台（UK ATC）。

更新时间： 2020-08-11

更新时间： 2021-09-06，由 B. Sargent，STScI 科学家，空间望远镜科学研究所（Space Telescope Science Institute）更新（添加了 MRS 模拟数据）。

更新时间： 2021-12-12，由 O. Fox，STScI 科学家更新（在笔记本中添加了黑体和多项式拟合）。

更新时间： 2024-10-29，由 C. Pacifici，STScI 数据科学家更新，适应 Jdaviz 4.0（仍需更新视频）。

看起来您没有提供任何内容需要翻译。如果您有特定的Markdown内容需要翻译，请将其粘贴在这里，我将很高兴为您提供帮助。

页面顶部

Contents