python 2d-gaussian fitting

#目的：对一个二维的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