Astropy.model 2DGaussian拟合FITS图像参数超限错误求助

Hey there, let's break down what's causing that error and fix your code step by step! That TypeError happens because you've accidentally created a Gaussian model with way more parameters than your number of data points—here's how to fix it:

Core Issues in Your Original Code #

• Incorrect Gaussian2D Initialization: You passed the entire data array as the first argument to Gaussian2D, which made Astropy create a model with thousands of amplitude parameters (one for every pixel) instead of a single 2D Gaussian with 5 standard parameters (amplitude, x/y mean, x/y stddev). This is why you got the "N must not exceed M" error—your model had 65541 parameters, but only 65536 data points.
• Wrong Input Dimensions for Fitting: 2D fitting requires gridded (2D) x/y coordinates, not 1D arrays. Your code used 1D x and y arrays, which the fitter couldn't map to your 2D image data.
• Unnecessary Manual Gaussian Calculation: You tried to fit your hand-crafted eq instead of the actual FITS image data—the fitter needs the raw image values to optimize the Gaussian model against.
• Minor Syntax/Import Errors: Missing imports for Circle, incorrect matplotlib import, and undefined ax variable would have broken your plotting even if the fitting worked.

Corrected Code #

Here's the fixed version that properly fits a 2D Gaussian to your FITS image:

import numpy as np
import astropy.io.fits as fits
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

# Load and preprocess FITS data
data = fits.getdata('AzTECC100.fits')
med = np.median(data)
data = data - med
# Extract the 2D slice from your data cube
data_2d = data[0, 0, :, :]

# Initialize Gaussian model with sensible starting parameters
# Find the brightest pixel to set initial center
max_val = np.max(data_2d)
y0, x0 = np.where(data_2d == max_val)
# Pick the first bright pixel if there are multiple with the same max
x0 = x0[0]
y0 = y0[0]
# Use global std dev as initial guess for sigma (or use a local region for better guess)
init_sigma = np.std(data_2d)

# Correctly initialize 2D Gaussian: amplitude, x_mean, y_mean, x_stddev, y_stddev
g_init = models.Gaussian2D(amplitude=max_val, x_mean=x0, y_mean=y0,
                           x_stddev=init_sigma, y_stddev=init_sigma)

# Create 2D grid coordinates for fitting
y_grid, x_grid = np.mgrid[:data_2d.shape[0], :data_2d.shape[1]]

# Run the fitter
fit_w = fitting.LevMarLSQFitter()
g_fit = fit_w(g_init, x_grid, y_grid, data_2d)

# Plot results
fig, ax = plt.subplots(figsize=(8, 5))
# Show the background-subtracted image
ax.imshow(data_2d, origin='lower', cmap='viridis')
# Overlay contour of the fitted Gaussian (half-max level)
ax.contour(x_grid, y_grid, g_fit(x_grid, y_grid), colors='white', levels=[max_val/2])
# Mark original peak and fitted center
ax.scatter(x0, y0, s=100, c='red', marker='x', label='Original Brightest Pixel')
ax.scatter(g_fit.x_mean.value, g_fit.y_mean.value, s=100, c='blue', marker='o', label='Fitted Gaussian Center')
# Add circle around fitted center
circle = Circle((g_fit.x_mean.value, g_fit.y_mean.value), 4, facecolor='none', edgecolor='blue', linewidth=2)
ax.add_patch(circle)

ax.legend()
plt.show()

# Print fitted parameters for reference
print("Fitted Gaussian Parameters:")
print(f"Amplitude: {g_fit.amplitude.value:.2f}")
print(f"X Center: {g_fit.x_mean.value:.2f}")
print(f"Y Center: {g_fit.y_mean.value:.2f}")
print(f"X Std Dev: {g_fit.x_stddev.value:.2f}")
print(f"Y Std Dev: {g_fit.y_stddev.value:.2f}")

Key Fixes Explained #

• Proper Gaussian Initialization: We now create a standard 2D Gaussian with only 5 adjustable parameters, which is way fewer than your number of data points.
• Gridded Coordinates: Using np.mgrid creates 2D arrays where each pixel has a corresponding (x,y) coordinate pair, which the fitter needs to map the model to your image.
• Fitting to Real Data: We pass the actual background-subtracted image (data_2d) to the fitter instead of your manual Gaussian calculation.

Bonus Tip #

Astropy v0.3.2 is extremely outdated (released in 2014!). If you can, upgrade to the latest version—you'll get better documentation, more fitting tools, and bug fixes that will make your workflow smoother.

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