<p>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 <a href="http://photutils.readthedocs.io/en/stable/detection.html" rel="nofollow noreferrer">example script</a> and to the <a href="http://photutils.readthedocs.io/en/stable/api/photutils.DAOStarFinder.html" rel="nofollow noreferrer">documentation</a> of DAOStarFinder. And here the example script explicitly:</p>

<pre><code>&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)
</code></pre>

CC BY-SA 4.0
2018-06-27T11:24:19.587Z
2018-06-27T11:24:19.587Z
NeStack
<p>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:</p>

<ol>
<li>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.</li>
</ol>

<p>If that doesn't work:</p>

<ol start="2">
<li><p>Get the mean value for each pixel across the time series, reject all pixels that are more that $x$ (%) away from it.</p></li>
<li><p>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.</p></li>
</ol>

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

CC BY-SA 4.0
2018-06-27T11:01:10.900Z
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.
NeStack
CC BY-SA 4.0
2018-06-18T17:48:09.873Z
2018-06-18T17:48:09.873Z
Alex
<p>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</p>

<p>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)?</p>

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

<hr>

<p>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</p>

<p><a href="https://i.sstatic.net/mgTV7.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/mgTV7.png" alt="In the left panel everything is normal, in the right panel, in the top right corner you see one of the artefacts, having a value below a saturation treshhold"></a></p>

CC BY-SA 4.0
2018-06-18T13:44:54.830Z
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
CC BY-SA 4.0
2018-06-18T17:45:21.800Z
Sounds like you're asking for a 3-D kernel, i.e. x,y,time, to apply in your search.
Carl Witthoft
CC BY-SA 4.0
2018-06-19T21:48:15.573Z
@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?
NeStack
CC BY-SA 4.0
2018-06-20T06:50:56.357Z
@Nestak Well... googling for "python astropy psf" has shown me that there is a PSF matching function in the photutils package of astropy.
Mick
CC BY-SA 4.0
2018-06-27T11:26:19.310Z
@Mick Thanks, googling solved the issue, indeed :) Based on your suggestion I wrote my own answer to the question
NeStack
CC BY-SA 4.0
2018-06-27T11:55:45.770Z
Cool! I'm glad you found a solution, +1 for sharing it.
Mick
<python><light-curve><fits-header>
Python: Find brightest pixel in fits-file and ignore brighter artefacts
NeStack
# Python: Find brightest pixel in fits-file and ignore brighter artefacts

(Asked by: username_0 on 2018-06-18)

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

[![In the left panel everything is normal, in the right panel, in the top right corner you see one of the artefacts, having a value below a saturation treshhold](https://i.sstatic.net/mgTV7.png)](https://i.sstatic.net/mgTV7.png)

**Tags:** python, light-curve, fits-header
> **Comments:**
> **username_1**: 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.  
> **username_2**: Sounds like you're asking for a 3-D kernel, i.e. x,y,time, to apply in your search.  
> **username_0**: @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?  
> **username_1**: @Nestak Well... googling for "python astropy psf" has shown me that there is a PSF matching function in the photutils package of astropy.  
> **username_0**: @Mick Thanks, googling solved the issue, indeed :) Based on your suggestion I wrote my own answer to the question  
> **username_1**: Cool! I'm glad you found a solution, +1 for sharing it.  

---

## 2 answers

### Answer by username_0 on 2018-06-27. Score: 1 | Accepted answer

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 by username_3 on 2018-06-18. Score: 0

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.

> **Comments:**
> **username_0**: 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.  

---

