Astropy.model 2dgaussian Issue

I'm trying to a 2D Gaussian to a fits image in order to find the brightest object in the image. I am trying to follow the 1D example here http://docs.astropy.org/en/v0.3.2/modelin

Solution 1:

There are several mistakes in your code but the one responsible for this particular crash is that you do not create the model correctly: the first argument to w = models.Gaussian2D() should be the amplitude and not data! Data have nothing to do with a "model" which is a theoretical dependence/function that you plan to fit to data.

Using data similarly generated as in https://stackoverflow.com/a/50560755/8033585, specifically,

y, x = np.indices((51,51))
data = 3 * np.exp(-0.7 * ((x - 24)**2 + (y - 26)**2))

and also commenting out lines that do not work, here is a modification of your code that should work (note: I increased the radius of the circle by 30x and intensity plot is on the log scale):

import astropy.io.fits as fits
import os
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt

#data = fits.getdata('AzTECC100.fits')
med = np.median(data) 
# data = data - med# data = data[0,0,:,:] # NOT SURE THIS IS NEEDED!

fit_w = fitting.LevMarLSQFitter()

y0, x0 = np.unravel_index(np.argmax(data), data.shape)
sigma = np.std(data)
amp = np.max(data)

w = models.Gaussian2D(amp, x0, y0, sigma, sigma)

yi, xi = np.indices(data.shape)

g = fit_w(w, xi, yi, data)
print(w)
model_data = g(xi, yi)

fig, ax = plt.subplots()
eps = np.min(model_data[model_data > 0]) / 10.0# use logarithmic scale for sharp Gaussians
ax.imshow(np.log(eps + model_data), label='Gaussian') 
circle = Circle((g.x_mean.value, g.y_mean.value),
                30 * g.x_stddev.value, facecolor ='none',
                edgecolor = 'red', linewidth = 1)

ax.add_patch(circle)
plt.show()

Printed model parameters are:

Model: Gaussian2D
Inputs: ('x', 'y')Outputs: ('z',)
Model set size: 1Parameters:
    amplitude x_mean y_mean     x_stddev        y_stddev    theta
    --------- ------ ------ --------------- --------------- -----
          3.024.026.00.08811846273940.08811846273940.0

---

NOTES:

You will find a great deal of examples of using astropy.modeling and an important discussion about the effect of the "scale" of the image data (order of magnitude) here: https://github.com/astropy/astropy/issues/6269 According to that issue, it may be helpful to "normalize" your image (rescale intensity to bring it into a "reasonable" range). Also follow links to notebooks that may provide further examples, such as https://github.com/ibusko/experiments_specviz/blob/master/experiment_normalization.ipynb Unfortunately these examples are for 1D data but can be easily expanded to 2D data.