# Python: Find brightest pixel in fits-file and ignore brighter artefacts

I have 4000 fits-files images from SPITZER, each containing the same star over a time period. Each image is a photmetry cube of 32 pixel * 32 pixel * 64 consequent time frames. I need to look through all the single frames and order the pixel by brightness. Not hard to do with Python, but sometimes I have in the frames short-lived events/artefacts which are brighter than the brightest pixel of the star (see below). These artefacts last only over one or two consequent time frames

How can I easily correct for the artefacts? Is there an astropy/pyfits package that does this? Or some kind of bad pixel flagging in SPITZER (I know that Kepler-images have flagging of bad pixel)?

Optional: It would be convenient to have the code in vectorized form, avoiding loops, so that I save computtional time

Fig: In the left frame everything is normal, in the right frame, in the top right corner you see one of the artefacts, having a value below a saturation treshhold

• The point spread function of the star should be consistent across the dataset, and artefacts will likely have different PSFs to the star. Instead of just picking out the brightest pixel in each frame (which is what I assume from the description), identify the star from its PSF and then select the peak values within that domain.
– Mick
Jun 18 '18 at 13:44
• Sounds like you're asking for a 3-D kernel, i.e. x,y,time, to apply in your search. Jun 18 '18 at 17:45
• @Mick Yes, I am indeed so far looking for the brightest pixel in the frames. And you are right, it will be better to look for the PSF, but I see no easy way to code this. Is there some function in astropy or pyfits that would recognize the PSF of the star and give me the coordinates of the PSF-center? Jun 19 '18 at 21:48
• @Nestak Well... googling for "python astropy psf" has shown me that there is a PSF matching function in the photutils package of astropy.
– Mick
Jun 20 '18 at 6:50
• @Mick Thanks, googling solved the issue, indeed :) Based on your suggestion I wrote my own answer to the question Jun 27 '18 at 11:26

Following Mick's comment, I found an astropy-function - DAOStarFinder - that does what I was looking for. It scans an image for light-sources and one can set threshhold for a detection, e.g. the "roundness", "sharpness" or brightness of an object have to be between [-0.05,0.05], [0.6,0.8] or >5*standard_deviation, respectively. I must say, that I have the impression the roundness and sharpness measure don't work consistently and sometimes values for the same object change enormously from one frame to the other. But DAOStarFinder gives you also a rough estimate of the flux density of the object, which I could use for sorting the sources and filtering out the relevant one. Here is a link to example script and to the documentation of DAOStarFinder. And here the example script explicitly:

>>> from astropy.stats import sigma_clipped_stats
>>> from photutils import datasets
>>> from photutils import DAOStarFinder

>>> data = hdu.data[0:400, 0:400]
>>> mean, median, std = sigma_clipped_stats(data, sigma=3.0, iters=5)
>>> daofind = DAOStarFinder(fwhm=3.0, threshold=5.*std)
>>> sources = daofind(data - median)

>>> print(sources)


Apart from the solution suggested in the comments (PSF fit), which is probably the proper way of doing further analyses, several solutions come to mind:

1. Check if the artefacts are in fact saturated pixels (sometimes the bias is subtracted, so you may need to check for saturation value minus bias). Ignore any pixels with that value.

If that doesn't work:

1. Get the mean value for each pixel across the time series, reject all pixels that are more that $x$ (%) away from it.

2. For each pixel across your time series, check if its brightness changes by more than $x$ (%). Note previous value, reject as long as subsequent values are also above the threshold.

All these methods will have problems if the bright pixel is at the position of the star.

• I was thinking of doing what you suggest, implementing some conditions and if clauses, which my python code has to check. The problem would be, that this will slow done my code, processing 260 000 data points, already running a few hours on a computational server. So I found an astropy function - DAOStarFinder - that does the work. Jun 27 '18 at 11:01