Cosmic Ray rejection for small data sets

[jjencson]

I am working with MOSFIRE data sets were we typically only have 6-10 frames (3-5 at a single dither position). The DRP defaults to using median combine for cosmic ray rejection in this case, which does not seem to work well enough. I can lower the threshold from >5 to >3 images for min/max CCR, but then I am throwing away ~ half my integration time. Has anyone had any success trying other modifications to do a better job, or implementing another method of CCR?

[KeckDRP]

I have had similar problems reducing data containing only a small numbers of observations. It's definitely a problem we need to work on and we will try to come up with a plan as soon as possible. We welcome suggestions on methods to improve the algorithms if you have any experience or ideas that we could try.

On Apr 10, 2015, at 2:49 AM, jjencson notifications@github.com wrote:

I am working with MOSFIRE data sets were we typically only have 6-10 frames (3-5 at a single dither position). The DRP defaults to using median combine for cosmic ray rejection in this case, which does not seem to work well enough. I can lower the threshold from >5 to >3 images for min/max CCR, but then I am throwing away ~ half my integration time. Has anyone had any success trying other modifications to do a better job, or implementing another method of CCR?

— Reply to this email directly or view it on GitHub.

[csteidel]

This is an area where the current DRP can definitely be improved; I have modified my local version of Background.py so that it has the following changes:

As long as the number of frames is >= 5, the "Sigclip CRR" algorithm is used. I also changed the code that estimates the std to exclude only the top and bottom value at each pixel, rather than the top and bottom 2 values. Note that the advantage here is that pixels are excluded only if they are outliers (the top and bottom are excluded only for the calculation of the mean and std, and in most cases they will be included in the stacked average.)

In my version, the min/max CRR is never used.

Note that none of the algorithms work well if the sky background is changing rapidly during the observation.

There are a number of algorithms out there for identifying CRs morphologically, which could be implemented during the relevant stage. CRR would also be improved by identifying the CR events after subtracting an estimated background from each individual frame first.

In any case, I have attached the section of code with the changes I have been using, which at least improves the situation when the number of frames is >= 5.

Appended is the relevant section of code (in Background.py)

---

Cosmic ray rejection code begins here. This code construction the

# electrons and itimes arrays.
if len(files) >= 5:
    print "Sigclip CRR"
    srt = np.argsort(electrons, axis=0, kind='quicksort')
    shp = el_per_sec.shape
    sti = np.ogrid[0:shp[0], 0:shp[1], 0:shp[2]]

    electrons = electrons[srt, sti[1], sti[2]]
    el_per_sec = el_per_sec[srt, sti[1], sti[2]]
    itimes = itimes[srt, sti[1], sti[2]]

    # Construct the mean and standard deviation by dropping the top and bottom two 
    # electron fluxes. This is temporary.
    mean = np.mean(el_per_sec[1:-1,:,:], axis = 0)
    std = np.std(el_per_sec[1:-1,:,:], axis = 0)

    drop = np.where( (el_per_sec > (mean+std*4)) | (el_per_sec < (mean-std*4)) )
    print "dropping: ", len(drop[0])
    electrons[drop] = 0.0
    itimes[drop] = 0.0

    electrons = np.sum(electrons, axis=0)
    itimes = np.sum(itimes, axis=0)
    Nframe = len(files) 

elif len(files) > 5:
    print "WARNING: Drop min/max CRR"
    srt = np.argsort(el_per_sec,axis=0)
    shp = el_per_sec.shape
    sti = np.ogrid[0:shp[0], 0:shp[1], 0:shp[2]]

    electrons = electrons[srt, sti[1], sti[2]]
    itimes = itimes[srt, sti[1], sti[2]]

    electrons = np.sum(electrons[1:-1,:,:], axis=0)
    itimes = np.sum(itimes[1:-1,:,:], axis=0)

    Nframe = len(files) - 2

else:
    print "WARNING: CRR through median filtering"
    for i in xrange(len(files)):
        el = electrons[i,:,:]
        it = itimes[i,:,:]
        el_mf = scipy.signal.medfilt(el, 5)

        bad = np.abs(el - el_mf) / np.abs(el) > 10.0
        el[bad] = 0.0
        it[bad] = 0.0

        electrons[i,:,:] = el
        itimes[i,:,:] = it

    electrons = np.sum(electrons, axis=0)
    itimes = np.sum(itimes, axis=0)
    Nframe = len(files)