如何将含光谱信息的XY坐标映射数据转换为字典（或替代结构）并实现光谱高斯拟合与振幅可视化？

Hey there! Let's tackle this problem step by step. Here's a comprehensive solution tailored to your needs, including data transformation, Gaussian fitting, and visualization:

You've already got the pandas import down—great start! Let's convert that tabular data into structures that scale well for 400+ pixels and 512-wave spectra.

Option 1: Dictionary of (X,Y) Tuples to Numpy Arrays #

This matches the structure you requested, and it automatically handles any number of pixels/spectrum lengths:

import pandas as pd
import numpy as np

# Your existing data import code
IN_PATH = r'PATH_TO_FILE'
FNAME = r'\FILENAME.txt'
data = pd.read_csv(IN_PATH + FNAME, sep='\t', skiprows=1)
data.columns = ["X", "Y", "Wave", "Intensity"]

# Build the spectrum dictionary
spectrum_dict = {}
for (x, y), group in data.groupby(['X', 'Y']):
    # Convert each group's Wave/Intensity pairs to a 2D numpy array
    spectrum = group[['Wave', 'Intensity']].to_numpy()
    spectrum_dict[(x, y)] = spectrum

This will capture every unique (X,Y) coordinate and store its full spectrum as a (n,2) array (where n=512 for your real data).

Option 2: 3D Numpy Array (Recommended for Large Datasets) #

Numpy arrays are faster for numerical operations like fitting, and they simplify plotting. Here's how to create a grid-aligned array:

# Get sorted unique coordinates to form a grid
unique_x = np.sort(data['X'].unique())
unique_y = np.sort(data['Y'].unique())
num_x, num_y = len(unique_x), len(unique_y)
spectrum_length = 512  # Use len(group) if your spectrum lengths vary (though you said they're fixed)

# Map coordinates to array indices
x_to_idx = {x: i for i, x in enumerate(unique_x)}
y_to_idx = {y: j for j, y in enumerate(unique_y)}

# Initialize empty array (shape: [Y rows, X columns, spectrum length, 2 (Wave/Intensity)])
spectrum_array = np.zeros((num_y, num_x, spectrum_length, 2), dtype=np.float64)

# Fill the array with spectrum data
for (x, y), group in data.groupby(['X', 'Y']):
    i = x_to_idx[x]
    j = y_to_idx[y]
    spectrum_array[j, i] = group[['Wave', 'Intensity']].to_numpy()

Next, let's fit each spectrum to a Gaussian distribution and extract the peak amplitude. We'll use scipy.optimize.curve_fit for this:

First, define the Gaussian function:

from scipy.optimize import curve_fit

def gaussian(x, amp, mean, std):
    """Gaussian function: amplitude, mean (center), standard deviation"""
    return amp * np.exp(-(x - mean)**2 / (2 * std**2))

Fitting with the Dictionary Approach #

amplitude_dict = {}
for (x, y), spectrum in spectrum_dict.items():
    wave = spectrum[:, 0]
    intensity = spectrum[:, 1]
    
    # Initial guess for fitting (max intensity, mid-wave value, rough spread)
    initial_guess = [np.max(intensity), np.mean(wave), np.std(wave)]
    
    try:
        # Fit the Gaussian to the data
        popt, _ = curve_fit(gaussian, wave, intensity, p0=initial_guess)
        amplitude_dict[(x, y)] = popt[0]  # Store the amplitude parameter
    except RuntimeError:
        # Handle cases where fitting fails (e.g., noisy spectra)
        amplitude_dict[(x, y)] = np.nan

Fitting with the Numpy Array Approach #

This is more efficient for large datasets:

amplitude_array = np.full((num_y, num_x), np.nan)

for j in range(num_y):
    for i in range(num_x):
        spectrum = spectrum_array[j, i]
        wave = spectrum[:, 0]
        intensity = spectrum[:, 1]
        
        initial_guess = [np.max(intensity), np.mean(wave), np.std(wave)]
        
        try:
            popt, _ = curve_fit(gaussian, wave, intensity, p0=initial_guess)
            amplitude_array[j, i] = popt[0]
        except RuntimeError:
            amplitude_array[j, i] = np.nan

Finally, plot the maximum amplitude for each (X,Y) coordinate using matplotlib:

For the Numpy Array Approach (Simpler) #

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 8))
# Use actual X/Y values for accurate axis labels
plt.pcolormesh(unique_x, unique_y, amplitude_array, cmap='viridis')
plt.colorbar(label='Gaussian Peak Amplitude')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('2D Map of Spectrum Gaussian Amplitudes')
plt.show()

For the Dictionary Approach #

First reshape the data into a grid, then plot:

# Extract coordinates and amplitudes
xs, ys = zip(*amplitude_dict.keys())
amps = list(amplitude_dict.values())

# Reshape into grid (assuming your coordinates form a regular grid)
xs = np.array(xs).reshape(num_y, num_x)
ys = np.array(ys).reshape(num_y, num_x)
amps = np.array(amps).reshape(num_y, num_x)

plt.figure(figsize=(10, 8))
plt.pcolormesh(xs, ys, amps, cmap='viridis')
plt.colorbar(label='Gaussian Peak Amplitude')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('2D Map of Spectrum Gaussian Amplitudes')
plt.show()

Notes for Scaling #

• The numpy array approach is better for 400+ pixels—it's faster for fitting and plotting.
• Add more error handling if needed (e.g., checking for NaN values in spectra before fitting).
• If your spectra have multiple peaks, you can modify the Gaussian function to a multi-peak version (e.g., sum of 2-3 Gaussians) and adjust the initial guess accordingly.

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