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

Question

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

Answer 1

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

>>> hdu = datasets.load_star_image()    
>>> 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)

Answer 2

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:

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

3. 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.