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
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:

&gt;&gt;&gt; from astropy.stats import sigma_clipped_stats
&gt;&gt;&gt; from photutils import datasets
&gt;&gt;&gt; from photutils import DAOStarFinder

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

&gt;&gt;&gt; print(sources)