[AstroPy] Gaussian models fitter non-functional?

[Erik Bray]

On 10/19/2015 05:56 PM, John Wainwright wrote:
> Just on the topic of PSF fitting, is there any plan to implement Moffat PSF
> models?  I'm trying to replicate some of the PSF fitting
> <http://pixinsight.com/doc/tools/DynamicPSF/DynamicPSF.html> capabilities found
> in PixInsight <http://pixinsight.com/> which provides excellent Moffat modeling
> support.

Well, there are are already these implementations:

http://docs.astropy.org/en/latest/api/astropy.modeling.functional_models.Moffat1D.html#astropy.modeling.functional_models.Moffat1D

http://docs.astropy.org/en/latest/api/astropy.modeling.functional_models.Moffat2D.html#astropy.modeling.functional_models.Moffat2D

So you can check how they compare.  You can also always add your own version:

http://docs.astropy.org/en/latest/modeling/new.html



>> On Oct 19, 2015, at 11:45 AM, Jeff Mangum <jmangum at nrao.edu
>> <mailto:jmangum at nrao.edu>> wrote:
>>
>> Thanks for the responses.  Adam Ginsburg and I iterated a bit on this problem
>> and found that the issue was the fact that the image I was fitting contained
>> nans.  Once I removed the nans the gaussian fit routine works.  Note that if
>> the fitter in astropy.models is not able to handle images with nans, then
>> perhaps it should issue a warning.
>>
>> -- Jeff
>>
>> On 10/19/15 11:00 AM, Erik Bray wrote:
>>> On 10/16/2015 09:47 AM, Giuliano Iorio wrote:
>>>>
>>>> Hi Jeff,
>>>> this is the kind of staff that usually happens if the position of the
>>>> initial guess is too far away the right one. Did you plot the contour of the
>>>> fitted gaussian on the image?
>>>> It can be also that the image is too noisy around the target to fit. I
>>>> usually use only a little subarray around the target, then I search the
>>>> position of the maximum and I set the initial guess of xo and yo to its
>>>> position. After I normalize all the subarray to the value of the maximum and
>>>> cut all the signal below a certain threshold (usually 0.1).
>>>> In this way I am able to obtain a robust fit also in very noisy map.
>>>
>>> I agree pretty much with all of Giuliano's advice.  And again it might help
>>> to exclude data outside a smaller subarray.  Without actually seeing the data
>>> it's hard to provide much other advice.
>>>
>>> Have you tried a different optimizer?  There are a couple other non-linear
>>> fitters in the package and some have different sensitivities to initial
>>> conditions.  But it's certainly not "non-functional".
>>>
>>> Erik
>>>
>>>>> Hi Jeff--
>>>>>
>>>>> The fitted model will not be in g_init (this remains unchanged), but in g.
>>>>>   If you take a look in g, you should be able to see the fitted results.
>>>>>
>>>>> Cheers
>>>>> Steve
>>>>>
>>>>> ----- Original Message -----
>>>>> From: "Jeff Mangum" <jmangum at nrao.edu <mailto:jmangum at nrao.edu>>
>>>>> To: astropy at scipy.org <mailto:astropy at scipy.org>
>>>>> Sent: Friday, October 16, 2015 3:11:56 PM
>>>>> Subject: [AstroPy] Gaussian models fitter non-functional?
>>>>>
>>>>> Hello,
>>>>>
>>>>> Long-time listener but first-time caller.  Have been trying to use the
>>>>> 2D Gaussian fitting capabilities within astropy.models to fit a 2D
>>>>> Gaussian to an image (in FITS format).  Here is my script:
>>>>>
>>>>> import numpy as np
>>>>> import radio_beam
>>>>> from astropy.modeling import models, fitting
>>>>> from astropy import units as u
>>>>> from astropy.io import fits
>>>>> import astropy.coordinates as coord
>>>>> import astropy.wcs as wcs
>>>>> import matplotlib
>>>>> import pylab as pl
>>>>>
>>>>> # The following is the position of the peak in both J1 and J2, so use
>>>>> # it as the starting guess
>>>>> ra0 = '13h15m03.50s'
>>>>> dec0 = '24d37m08.2s'
>>>>>
>>>>> radeg = coord.Angle(ra0, unit=u.hour).degree
>>>>> decdeg = coord.Angle(dec0, unit=u.degree)
>>>>>
>>>>> # Now open the image
>>>>> hdulist = fits.open('IC860CbandCarrayH2COJ1_moment0.fits')
>>>>> #hdulist = fits.open('IC860KubandDarrayH2COJ2_moment0.fits')
>>>>> w = wcs.WCS(hdulist[0].header, hdulist)
>>>>>
>>>>> yy,xx = np.indices(hdulist[0].data.shape)
>>>>>
>>>>> x0,y0 = w.wcs_world2pix(radeg, decdeg, 1)
>>>>>
>>>>> # Now fit 2D gaussian
>>>>> g_init = models.Gaussian2D(amplitude=1., x_mean=x0, y_mean=y0,
>>>>> x_stddev=1, y_stddev=1)
>>>>> fit_g = fitting.LevMarLSQFitter()
>>>>> g = fit_g(g_init, xx, yy, hdulist[0].data)
>>>>>
>>>>>
>>>>> When I look at g_init (which is where I believe the fit results should
>>>>> be), I see that the fit has not moved from its initial settings.  It
>>>>> appears that the fitter just did not run.  With Adam Ginsburg's help I
>>>>> have a workaround using gaussfitter, but would certainly like to make
>>>>> this work with astropy.models.  Thanks!
>>>>>
>>>>> -- Jeff