Use imexamine automatically given an image and a list of x,y coordinates

[Gabriel-p]

I'd like to use imexamine to emulate IRAF's psfmeasure task, which is useful to estimate the average FWHM of an image.

To do this, I need to pass to imexam an image, and a list of x,y star center coordinates. Then I have to emulate pressing the a key on each of those x,y coordinates, and log the results either to screen or a file.

Is this possible? Or must I press the a key manually star after star to obtain the estimated FWHMs?

[sosey]

nope, you can just feed the data and x,y positions in and save the information.... You probably don't need the full aperture photometry though, yes? You can save the aperture photometry similar to the example below, but if you're mostly interested in the gauss fit results, you can use the Gaussian line fit function in either x or y ( attached to the j and k keys) and return the fit directly, then save however you like. It could look something like this:

from imexam.imexamine import Imexamine
from astropy.io import fits

#get your data into an array
data=fits.getdata(filename)

#get an object up with your data attached
plots=Imexamine()
plots.set_data(data)

If you have a text file with columns of x, y, star id, then one option is to read the information and save the results as a dictionary, you might also try an astropy table, here's a simple read and save, just as a quick example

starlist=list()
with open(filename) as fh:
    for line in fh:
        starlist.append(tuple(map(float,line.strip().split())))

Now you have a list of tuples which are the float values of x,y,star id You can continue and feed this to the plots object and save to a dictionary. I'll be a little obtuse here and just save a dictionary indexed by star id with a list of values in order x, y, std, mean, amp, fwhm If you were doing this in an astropy table you could do math on the columns

#make a function to calculate the fwhm
g = lambda x: x * np.sqrt(8.0 * np.log(2.))

or you can import the function from imexam, note that this function currently returns a tuple (fwhm,fwhm) because it supports accepting the stddev from x and y assumes they are equal when only 1 value is given

from imexam.math_helper import gfwhm
fwhm=gfwhm(stddev)[0]  #to get just the x fwhm

#continuing on ...

results=dict()
for star in starlist:
    x,y,sid=star
    gauss=plots.line_fit(x,y,genplot=False)
    results[sid]=[x, y, gauss.stddev, gauss.mean, gauss.amplitude, g(gauss.stddev)]

Now you have a dictionary, indexed by the star id with all the information, you can use this to make plots, for example plot of the x vs fwhm to see how it changes across your image, or calculate the average fwhm for all the stars... you can save the results in any form you like :-)

As long as centering is on, line_fit and column_fit will first center the object using a gaussian 2D fit and then calculate the 1D fit at that center in either the x or y direction.

Does this example help? I believe all/most of the functions accept genplot=False and instead return useful information.

[Gabriel-p]

Thank you very much for your detailed answer @sosey! A couple of questions if you don't mind:

1. Importing imexam, even if it is not used, makes the plt.show() window crash (or immediately close). Try the following small script for example. In my case, the plot window never appears (I think it appears and is closed so quickly I can't even see it).

   import matplotlib.pyplot as plt
   from imexam.imexamine import Imexamine
   
   plt.scatter([2., 3.], [2., 3.])
   plt.show()
   

   I'm using Python 2.7.12 and imexam v0.6.2 (installed via AstroConda).

2. Why does line_fit requires the x,y coordinates as integers? Doesn't it make more sense to accept floats?

3. The line_fit Notes state: "If centering is True in the parameter set". Where is this "centering" parameter? It is not listed as a parameter from line_fit.

4. Calling plots.line_fit(x, y, genplot=False) writes the xc, yc values to screen. It is a minor annoyance, but I'd like to turn this off if possible.

5. If I understand correctly, plots.line_fit() returns the 1D FWHM of a star in the x axis. To obtain the y axis FWHM I'd have to use

   gauss_y = plots.column_fit(x, y, genplot=False)
   fwhm_y = gfwhm(gauss_y.stddev)[0]
   

   Is this correct?

6. Sometimes I get negative FWHMs using the method described in your answer. What is the meaning of such values?

[sosey]

I'm in the middle of a conference this week @Gabriel-p, but I will get to these questions as soon as I can. In the mean time, if you are using a conda virtual environment, will you dump it to a text file and send me the result so I can repilcate it? conda env export > environment.yml

[Gabriel-p]

Of course @sosey, take as much time as you need.

Here's the environment.yml.zip file. Let me know if there's anything else you need.

In the meantime, I'll add another question:

7- What FWHM does the a key return? The averaged x, y (line, column) FWHMs? The FWHM of a 2D fitted Gaussian?

[sosey]

Thankyou. The "a" key is returning photometry from photutils, but I use the sigma from the 2D gaussian fit which recenters the star to calculate the FWHM

You can call this function yourself, imexam.math_helper.gauss_center, it takes a data array which is assumes is centered on the object and then calculates a better 2d center with the fit.

3. the line fit parameters, these are the user configurable dictionaries of parameters for each function, not the actual parameter inputs to the function, I'll try and make this more clear in the documentation

In [4]: plots.line_fit_pars

{'background': [False, 'Solve for background? [bool]'],
 'center': [True, 'Recenter around the local max'],  <------------center using 2D gaussian
 'func': ['Gaussian1D', 'function for fitting (see available)'],
 'function': ['line_fit'],
 'logx': [False, 'log scale x-axis?'],
 'logy': [False, 'log scale y-axis?'],
 'order': [1, 'Polynomial order to fit, 1=line'],
 'pointmode': [True, 'plot points instead of lines? [bool]'],
 'rplot': [15, 'Plotting radius in pixels'],
 'title': [None, 'Title of the plot'],
 'width': [10.0, 'Background  width in pixels'],
 'xlabel': ['Column', 'The string for the xaxis label'],
 'ylabel': ['Flux', 'The string for the yaxis label']}

[Gabriel-p]

Thanks @sosey. Updating status of the questions:

1. Still an issue.

2. Resolved. I see that it expects the column/line number, which is why it needs an integer.

3. Resolved.

4. I could "fix" this wrapping calls to column_fit or line_fit around

   sys.stdout = open(os.devnull, "w")
   ...
   sys.stdout = sys.__stdout__
   

   It feels hackish though.

5. Unless I understood something wrong, this is correct.

6. Here's a comparision of line_fit and column_fit FWHM values, with those returned by psfmeasure for the exact same group of stars:

   # psfmeasure  gfwhm_line_fit  gfwhm_column_fit
   4.684   -18.212 34.701
   4.992   202.577 7.016
   0.906   39.609  -585.165
   1.018   23.748  -36.779
   0.809   -1.603  50.768
   5.428   6.345   5.861
   4.948   0.853   4.783
   5.076   -5.210  5.225
   4.352   3.209   -4.828
   5.268   1.266   0.450
   4.392   -5.525  4.550
   5.268   1.266   0.450
   4.496   -5.193  6.898
   5.588   5.948   5.174
   5.588   5.948   5.174
   5.068   -5.683  4.860
   5.096   5.544   5.346
   5.096   5.544   5.346
   5.024   -0.357  -5.041
   4.052   -5.909  -4.805
   5.024   -0.357  -5.041
   5.024   -0.357  -5.041
   4.052   -5.909  -4.805
   5.5 -5.429  5.551
   2.244   13339.176   628.546
   5.26    -5.435  5.176
   5.716   7.686   -2.709
   0.969   1.051   5.259
   4.82    0.723   -4.977
   3.856   -5.134  4.826
   3.414   4.721   4.228
   3.594   0.648   0.689
   4.376   5.480   7.224
   5.04    -83.771 -5.033
   4.628   5.032   4.890
   4.956   -5.244  -5.252
   4.908   6.064   -5.513
   4.908   6.064   -5.513
   4.672   -0.698  4.755
   0.993   1.405   1.056
   4.756   -5.871  -5.100
   5.3 5.808   -6.084
   5.032   6.905   -0.034
   4.864   -5.313  0.890
   4.84    -5.099  4.864
   0.941   1.067   0.181
   5.292   -5.882  0.643
   4.716   5.091   4.874
   4.308   -5.071  -4.741
   4.476   4.918   4.941
   6.02    -5.945  -5.896
   1.244   33.442  0.261
   5.276   5.743   -5.184
   4.712   5.154   4.880
   2.796   4.779   -4.766
   5.232   6231.627    -5.071
   4.912   4.910   5.102
   4.82    5.272   0.873
   4.788   5.600   5.024
   4.032   -0.401  4.818
   

   I understand that x,y FWHMs can differ a bit, but I don't understand the negative values, or when line_fit and column_fit differ by a lot.

7. "The "a" key is returning photometry from photutils", could you expand a bit on how the FWHM is calculated with a? Is this a 2D Gaussian (or Moffat, etc) fit to the data?

[sosey]

@Gabriel-p Some more answers:

2) Why does line_fit requires the x,y coordinates as integers? Doesn't it make more sense to accept floats?

This is to accommodate the imexam loop. The integer pixel location of the pointer is used to make a subarray cut of the data to localize the fit, then with centering on, a better floating point center is found assuming that the object is close to the center of the smaller array. This gives a little more flexibility for people visually interacting with an image and clicking on things within a few pixels of the objects true center. In reality, you could pass a float number to line_fit, but it will get converted to make this smaller array from the data ...actually in looking at the code I see a bug here too where the int isn't getting applied, so I'm going to fix that :) So my first suggestion is to pass a float or int to line_fit and with centering on you should be fine. Secondarily, I would suggest using the astropy models directly, you can take a look at the code in the math_helper.gauss_center function, ignore the first part that guesses the center, and use the rest of the code to do the 2d fits for yourself.

4) Calling plots.line_fit(x, y, genplot=False) writes the xc, yc values to screen. It is a minor annoyance, but I'd like to turn this off if possible.

I think this is part of the logging that I have enabled. Currently a log message at debug level is sent to the screen, you could set the logging to a higher level. From the full object, if you are using the connect() the loggins sets up and there is a setlog() function which you can pass the logging level to, setlog(level=logging.WARNING) I think would turn off the print to screen. But this logging function isn't attached to the Imexamine() object. Let me think about the best way to do this a little bit.

If I understand correctly, plots.line_fit() returns the 1D FWHM of a star in the x axis. To obtain the y axis FWHM I'd have to use column_fit

yup, these are the 1D, line or column fit functions

5) gauss_y = plots.column_fit(x, y, genplot=False)
fwhm_y = gfwhm(gauss_y.stddev)[0]
Is this correct?

that should be correct

6) Sometimes I get negative FWHMs using the method described in your answer. What is the meaning of such values?

Most likely the centering is off and the model fit is barfing, often getting good starting conditions can be tricky and data dependent. I have them set to try and cover the average dataset. Since you are doing something very specific I think you should consider just using the astropy models directly like discussed above? If you run into issues there, you can send me more questions, or ask on the astropy models page.

[sosey]

I had the above typed up when your update came in, lemme think on the psfmeasure differences a bit and I'll get back to you