Python光谱绘图异常：25组观测数据叠加出现背景干扰

我有25组光谱观测数据，每组包含29个阶（order，行）和1973个像素（pixel，列，原始2047个已剔除坏点和异常值）。已完成质心修正（Barycentric correction），需要绘制指定阶的波长（Tw）与强度（Tint）曲线并叠加25组数据。多数数据能正常叠加，但出现异常背景图，无法定位原因。

原始代码（存在多处潜在问题） #

import numpy as np
import matplotlib.pyplot as plt
import h5py
import scipy.io as sio
from scipy.interpolate import interp1d

# Corrected the number of pixels
num_pixels = 1973

# Corrected the number of orders
num_orders = 29

# Load data from the "Prepared.mat" file
prepared_mat_file_path = '/Users/xxx/Desktop/Star/Prepared.mat'
data = sio.loadmat(prepared_mat_file_path)

# Specify the path to the HDF5 file
hdf5_file_path = '/Users/xxx/Desktop/Star/bcorr_fib4_all.mat'

# Open the HDF5 file for reading
with h5py.File(hdf5_file_path, 'r') as hdf5_file:
    # Access the dataset you want
    bcorr_fib4_all = hdf5_file['bcorr_fib4_all'][:]

# Access the data from "Prepared.mat"
# Trim the data to 1973 pixels
ALLfib4spec = data['ALLfib4spec'][:, :num_pixels, :]
nALLfib4specbNc = data['nALLfib4specbNc'][:, :num_pixels, :]
ALLfib4specwave = data['ALLfib4specwave'][:, :num_pixels, :]

# Generate some example data (replace with your actual data)
# Note: You should replace this random data with your actual data.
wave = np.random.rand(num_orders, num_pixels)
bcorr_fib4_all = np.random.rand(num_observations)


# Barycentric correction
# Note: Ensure that the barycentric correction is being applied correctly.

# Initialize ALLwv with the appropriate dimensions
ALLwv = np.zeros((num_orders, num_pixels, num_observations))

for s in range(ALLfib4spec.shape[2]):
    tWVout = ALLfib4specwave[:, :, s]
    for ord in range(tWVout.shape[0]):
        tWVout[ord, :] = shift_vel(tWVout[ord, :], bcorr_fib4_all[s] * -1)
    # dimensions are: [orders pixels observations]
    ALLwv[0:tWVout.shape[0], :, s] = tWVout  # Assign values to the corresponding orders in ALLwv

# Create the log step wavelength axis so each pixel is the same velocity size
c = 299792458
logstep = 1e-6
a = 8
velstep = (np.exp(a + logstep) - np.exp(a)) / np.exp(a) * c
lw = np.arange(np.log(min(Tw[:, 1]) + 2), np.log(max(Tw[:, 1]) - 2), logstep)

# Define the order you want to plot (e.g., order 16)
order_to_plot = 16

# Extract wavelengths for the chosen order
Tw = np.transpose(ALLwv[order_to_plot, :, :])
Tint = np.transpose(nALLfib4specbNc[order_to_plot, :, :])

# Number of observations
num_observations = Tint.shape[0]

# Create a plot
plt.figure(figsize=(10, 6))

# Loop through observations and plot each one
for observation_index in range(num_observations):
    wavelengths = Tw[observation_index, :]
    intensities = Tint[observation_index, :]

    # Plot the observation
    plt.plot(wavelengths, intensities, label=f'Observation {observation_index + 1}')

# Add labels and title
plt.xlabel('Wavelength')
plt.ylabel('Intensity')
plt.title(f'Wavelength vs. Intensity for Order {order_to_plot}')
plt.legend()  # Add a legend to distinguish observations
plt.show()

核心问题排查与修正 #

1. 变量定义顺序错误

原始代码中num_observations在使用后才定义，直接导致生成随机数据时报错，需提前从数据中提取观测数量：

# 提前获取观测数量
num_observations = ALLfib4spec.shape[2]

2. 缺失质心修正核心函数

代码调用了shift_vel但未实现，需补充多普勒位移计算逻辑：

def shift_vel(wavelength, vel):
    c = 299792458  # 光速（m/s）
    return wavelength * (1 + vel / c)  # 基于多普勒效应的波长偏移

3. 无效的对数波长轴代码

lw计算用到了未定义的Tw，且逻辑与后续绘图无关，建议删除或重构。

4. 数据维度与赋值错误

修正ALLwv的赋值逻辑，避免修改原始数据：

for s in range(num_observations):
    tWVout = ALLfib4specwave[:, :, s].copy()  # 复制数据避免修改原始值
    for ord in range(num_orders):
        tWVout[ord, :] = shift_vel(tWVout[ord, :], bcorr_fib4_all[s] * -1)
    ALLwv[:, :, s] = tWVout

5. 绘图优化（解决异常背景）

• 过滤异常强度值：用np.clip限制强度范围
• 降低曲线透明度：减少重叠导致的背景混乱
• 调整图例位置：避免遮挡曲线

修正后可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
import h5py
import scipy.io as sio
from scipy.interpolate import interp1d

# 配置参数
num_pixels = 1973
num_orders = 29

# 加载数据
prepared_mat_file_path = '/Users/xxx/Desktop/Star/Prepared.mat'
data = sio.loadmat(prepared_mat_file_path)

hdf5_file_path = '/Users/xxx/Desktop/Star/bcorr_fib4_all.mat'
with h5py.File(hdf5_file_path, 'r') as hdf5_file:
    bcorr_fib4_all = hdf5_file['bcorr_fib4_all'][:]

# 裁剪并提取数据
ALLfib4spec = data['ALLfib4spec'][:, :num_pixels, :]
nALLfib4specbNc = data['nALLfib4specbNc'][:, :num_pixels, :]
ALLfib4specwave = data['ALLfib4specwave'][:, :num_pixels, :]

# 获取观测数量
num_observations = ALLfib4spec.shape[2]

# 实现质心修正的波长偏移函数
def shift_vel(wavelength, vel):
    c = 299792458
    return wavelength * (1 + vel / c)

# 初始化修正后的波长数组
ALLwv = np.zeros((num_orders, num_pixels, num_observations))

# 应用质心修正
for s in range(num_observations):
    tWVout = ALLfib4specwave[:, :, s].copy()
    for ord in range(num_orders):
        tWVout[ord, :] = shift_vel(tWVout[ord, :], bcorr_fib4_all[s] * -1)
    ALLwv[:, :, s] = tWVout

# 指定绘制阶数
order_to_plot = 16

# 提取目标阶数据
Tw = np.transpose(ALLwv[order_to_plot, :, :])
Tint = np.transpose(nALLfib4specbNc[order_to_plot, :, :])

# 绘图
plt.figure(figsize=(12, 6))
for idx in range(num_observations):
    wavelengths = Tw[idx, :]
    intensities = Tint[idx, :]
    # 过滤异常强度值
    intensities = np.clip(intensities, 0, np.percentile(intensities, 99))
    # 降低透明度减少重叠干扰
    plt.plot(wavelengths, intensities, alpha=0.6, label=f'观测{idx+1}')

plt.xlabel('波长')
plt.ylabel('强度')
plt.title(f'第{order_to_plot}阶：波长-强度叠加曲线')
plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left')
plt.tight_layout()
plt.show()

内容的提问来源于stack exchange，提问作者AlbPet