#目的：对一个二维的data做二维高斯拟合
#需加载的包，另外还需要numpy等
from astropy.io import fits
from astropy.modeling import Fittable2DModel, Parameter, models, fitting
from astropy.wcs import WCS
import astropy.units as u
#读入XXXX.fits
fitsname = ‘XXXX.fits’
FITS_1 = fits.open(fitsname)
data = FITS_1[0].data
hdr = FITS_1[0].header
wcs = WCS(m0_fits[0].header)
#二维高斯模型&输入参数初值和限制范围：
g_init = models.Gaussian2D(amplitude=amp_start, x_mean=x_mu_start, y_mean=y_mu_start, x_stddev=x_sig_start, y_stddev=y_sig_start)
g_init.amplitude.fixed=True #固定amp值
g_init.x_mean.bounds=[x_mu_start-5,x_mu_start+5] #限制x_mu变化范围，下同
g_init.y_mean.bounds=[y_mu_start-5,y_mu_start+5]
g_init.x_stddev.bounds=[0.1,5]
g_init.y_stddev.bounds=[0.1,5]
#建立一个用于拟合的数据块（fit_data，选择3倍rms以上的信号用于拟合）
ny, nx = data.shape
fit_y, fit_x = np.mgrid[ :ny, :nx]
fit_data = data.copy()
fit_data[fit_data < 3*rms] = 0
#2d拟合开始，使用fitting.LevMarLSQFitter()
fit = fitting.LevMarLSQFitter()
#fitting.LevMarLSQFitter.supported_constraints
fitted_model = fit(g_init, fit_y, fit_x, fit_data)
#拟合参数：
Intensity = fitted_model.amplitude.value
x_fwhm = fitted_model.x_stddev.value*2.*np.sqrt(2.*np.log(2.))#pixel
y_fwhm = fitted_model.y_stddev.value*2.*np.sqrt(2.*np.log(2.))#pixel
x_center = fitted_model.x_mean.value
y_center = fitted_model.y_mean.value
theta = fitted_model.theta.value
#查看参数
print(fitted_model)
P.S. wcs->pix & pix->wcs
co = SkyCoord(‘XXhXXmXXs’, ‘+YYdYYmYYs’)
wcs.world_to_pixel(co) #wcs->pix
wcs.wcs_pix2world([[y_pix,x_pix]], 0) #pix->wcs