prose icon indicating copy to clipboard operation
prose copied to clipboard

Rejecting saturated stars based on custom value

Open jpdeleon opened this issue 6 months ago • 1 comments

I am trying to reduce LCO/MuSCAT3 & 4 data using prose. It has 4 different cameras used simultaneously and each one has different saturation level which depends on the mode of operation. The saturation values saved in the header are somehow wrong so I figured to simply re-define the saturation value by instantiating the telescope class.

saturation = 50_000 #example value
telescope = Telescope(saturation=saturation)
ref = FITSImage(ref_filename, telescope=telescope)

However, is it recommended to instantiate the telescope class inside an aperture photometry sequence? I thought it would be simpler if somehow there is an additional saturation argument in PointSourceDetection block but I couldn't work it out. Alternatively, I created a custom block to reject stars given saturation level and overwrite the sources detected in the previous block. For example

class RejectSaturatedSources(Block):
    def __init__(self, saturation_level=None, name=None):
        """Filter out saturated sources from the detected sources.

        Parameters
        ----------
        saturation_level : float, optional
            The level above which a star is considered saturated.
            If None, will use image.saturation_level if available.
        name : str, optional
            Name of the block, by default None
        """
        super().__init__(name=name)
        self.saturation_level = saturation_level

    def run(self, image):
        if len(image.sources) == 0:
            return

        # Determine saturation level
        sat_level = self.saturation_level
        if sat_level is None and hasattr(image, 'saturation_level'):
            sat_level = image.saturation_level

        if sat_level is None:
            raise ValueError("No saturation level provided or found in image")

        # Filter out sources with peak values above saturation level
        non_saturated = [s for s in image.sources if s.peak < sat_level]
        image.sources = Sources(non_saturated, type=image.sources.type)

and then it can be inserted after the PointSourceDetection block

calibration = Sequence(
        [
            blocks.Trim(),
            blocks.PointSourceDetection(),  
            RejectSaturatedSources(saturation), 
            blocks.Cutouts(),  
            blocks.MedianEPSF(),  
            blocks.psf.Moffat2D(),  
            blocks.CentroidQuadratic(),  
            blocks.AperturePhotometry(),
            blocks.AnnulusBackground(),  
        ]
    )
calibration.run(ref, show_progress=False)

What is the correct way to do this? On the other hand, in case I need to define a .telescope file in built-ins, should I create one .telescope file per camera?

jpdeleon avatar Jun 16 '25 06:06 jpdeleon

Thanks for your question @jpdeleon. I think this way is totally valid. Indeed, you would have to define a .telescope file for each camera.

lgrcia avatar Jun 19 '25 23:06 lgrcia