Difference between deep_coadd and a custom coadd

Science
, ,
#1

Custom coadd built from the same visits as the DP1 deep_coadd does not match it — how should two coadds be compared?

Hi everyone,

I am testing the custom coadd workflow from tutorial 105.6 “Create custom coadd images” on the Rubin Science Platform (DP1, LSST Science Pipelines v30.0.11, large container). As a sanity check, instead of selecting a subset of visits, I rebuilt the coadd using exactly the same input visits that were used for the official DP1 deep_coadd, expecting to get (nearly) the same image. However, when I subtract the two images, they are clearly not identical.

What I did

  • Tract / patch / band: tract = 2394, patch = 3, band = 'i', skymap = 'lsst_cells_v1'
  • Input visits: the 27 visits from deep_coadd.getInfo().getCoaddInputs().visits
  • Config overrides (as in the tutorial):
pipeline.addConfigOverride('makeDirectWarp', 'useVisitSummaryPsf', False)
pipeline.addConfigOverride('makeDirectWarp', 'useVisitSummaryPhotoCalib', False)
pipeline.addConfigOverride('makeDirectWarp', 'useVisitSummaryWcs', False)
pipeline.addConfigOverride('makeDirectWarp', 'connections.calexp_list', 'visit_image')
pipeline.addConfigOverride('makeDirectWarp', 'doSelectPreWarp', False)
  • Query string:
query_string = f"tract = {my_tract} AND patch = {my_patch} AND visit IN {my_visits_tupleString} " \
               f"AND skymap = 'lsst_cells_v1' AND band = 'i'"

Then I retrieved both images and compared them:

deepCoadd_mine = out_butler.get('deep_coadd_predetection', tract=my_tract, patch=my_patch,
                                band='i', instrument='LSSTComCam', skymap='lsst_cells_v1',
                                collections=executor.quantum_graph.metadata["output_run"])
deepCoadd_rubin = butler.get('deep_coadd', tract=my_tract, patch=my_patch,
                             band='i', instrument='LSSTComCam', skymap='lsst_cells_v1')

diff = deepCoadd_mine.clone()
diff.image.array[:] = deepCoadd_mine.image.array - deepCoadd_rubin.image.array

Results

Side by side, both coadds look very similar (same display scale):

coadd_comparison
coadd_comparison1485×643 145 KB

Simple image comparison

But a simple image comparison at the array level (mine - rubin, displayed between the 1st and 99th percentiles) shows:

coadd_difference
coadd_difference693×636 55.5 KB

Alard-Lupton Substraction

We also tried to implement a difference function using ´AlardLuptonSubtractTask´. The approach was adapted from this forum thread, where the DP1 ´deep_coadd´ is used as the template and my custom coadd as the science image (See code below.) . But it also gives some residuals:

coadd_difference_alardlupton
coadd_difference_alardlupton803×636 21.6 KB

Questions

  1. Is there a recommended method or tool in the Science Pipelines to compare two coadds of the same tract/patch beyond a direct pixel subtraction? Apart from our sanitiy check that our custom coadd reproduces de official one, we would like to compare custom coadds with injected images against the provided deep coadds, for instance.
  2. Is it expected that the two coadds differ?

Thanks!

def diff_AlardLupton(templateExposure, scienceExposure, warp=True):
    """
    The idea was take from:
    https://community.lsst.org/t/issues-with-image-subtraction/10429
    """
    import numpy as np
    import lsst.afw.table as afwTable
    import lsst.daf.base as dafBase
    from lsst.meas.algorithms.detection import SourceDetectionTask
    from lsst.meas.deblender import SourceDeblendTask
    from lsst.meas.base import SingleFrameMeasurementTask
    from lsst.ip.diffim import AlardLuptonSubtractTask, AlardLuptonSubtractConfig
    
    # Remove the bit mask
    print("===> Cleaning masks")
    template_masks = set(templateExposure.mask.getMaskPlaneDict())
    scienceExposure_masks = set(scienceExposure.mask.getMaskPlaneDict())
    common_masks = template_masks & scienceExposure_masks

    for exp in [templateExposure, scienceExposure]:
        for mask in list(exp.mask.getMaskPlaneDict().keys()):
            if mask not in common_masks:
               print(f"masks: {mask}")
               try:
                   exp.mask.removeAndClearMaskPlane(mask)
               except Exception:
                   pass
    
    # Make a schema
    # Identified sources on the scienceExposure. 
    # This catalog is used to select sources to perform 
    # the AL PSF matching on stamp images around them.
    schema = afwTable.SourceTable.makeMinimalSchema()
    schema.addField("coord_raErr", type="F")
    schema.addField("coord_decErr", type="F")
    schema.addField("detect_isPrimary", type="F")  # Flag
    schema.addField("sky_source", type="Flag")
    algMetadata = dafBase.PropertyList()

    # Detection task
    detectConfig = SourceDetectionTask.ConfigClass()
    detectConfig.thresholdValue = 5
    detectConfig.thresholdType = "stdev"
    detectTask = SourceDetectionTask(schema=schema, config=detectConfig)

    # here because schema is modify by detectTask.run
    deblendTask = SourceDeblendTask(schema=schema)
    
    measConfig = SingleFrameMeasurementTask.ConfigClass()
    measTask = SingleFrameMeasurementTask(schema=schema, config=measConfig, algMetadata=algMetadata)

    print("===> Run detection")
    tab = afwTable.SourceTable.make(schema)
    result = detectTask.run(tab, scienceExposure)
    sources = result.sources

    # Deblend + Measurement
    print("===> Run Deblend and Measurement")
    deblendTask.run(scienceExposure, sources)
    measTask.run(measCat=sources, exposure=scienceExposure)
    sources = sources.copy(True)

    # Sky sources flag
    print("===> Selecting sky sources")
    select_sky_sources(sources, schema)
    
    # Warp the templateExposure to match with scienceExposure (imagen + PSF)
    if warp:
        print("===> Warping")
        warp_templateExposure = warp_img(
            ref_img=scienceExposure,
            img_to_warp=templateExposure)
    else:
        warp_templateExposure = templateExposure

    # Subtraction task 
    print("===> Subtracting")
    # https://pipelines.lsst.io/py-api/lsst.ip.diffim.AlardLuptonSubtractTask.html#lsst.ip.diffim.AlardLuptonSubtractTask.run
    config = AlardLuptonSubtractConfig()
    subtractTask = AlardLuptonSubtractTask(config=config)
    
    result = subtractTask.run(
        template=warp_templateExposure,
        science=scienceExposure,
        sources=sources
    )
    diff = result.difference
    return diff


def select_sky_sources(sources, schema, max_sky=500):
    import random

    sky_source_key = schema["sky_source"].asKey()
    selected = []
    for record in sources:
        try:
            flux = record.get("base_PsfFlux_instFlux")
            flux_err = record.get("base_PsfFlux_instFluxErr")

            snr = flux / flux_err if flux_err > 0 else 0

            psf_flag = record.get("base_PsfFlux_flag")
            saturated = record.get("base_PixelFlags_flag_saturated")
            bad = record.get("base_PixelFlags_flag_bad")
            edge = record.get("base_PixelFlags_flag_edge")
            nchild = record.get("deblend_nChild")
            
            is_sky = (
                abs(snr) < 20 and      # 2 stricto
                not saturated and
                not bad and
                not edge and
                nchild == 0
            )
        except Exception:
            is_sky = False

        record.set(sky_source_key, is_sky)
        if is_sky:
            selected.append(record)

    print(f"Sky sources seleccionadas: {len(selected)}")

    # FALLBACK SI NO HAY
    if len(selected) == 0:
        print("No sky sources → usando fallback aleatorio")
        n = len(sources)
        fallback_idx = random.sample(range(n), max(5, int(0.2 * n)))

        for i, record in enumerate(sources):
            record.set(sky_source_key, i in fallback_idx)
    elif len(selected) > max_sky:
        # limitar número
        keep = set(random.sample(range(len(selected)), max_sky))
        for i, record in enumerate(selected):
            if i not in keep:
                record.set(sky_source_key, False)

def warp_img(ref_img, img_to_warp, warping_kernel="lanczos5"):
    """
    Warp an exposure (image + PSF) onto the coordinate system of another.
    """
    import copy
    import lsst.afw.math as afwMath
    import lsst.afw.geom as afwGeom
    import lsst.meas.algorithms as measAlg

    # warp imagen
    config = afwMath.Warper.ConfigClass()
    config.warpingKernelName = warping_kernel
    warper = afwMath.Warper.fromConfig(config)

    bbox = ref_img.getBBox()
    warpedExp = warper.warpExposure(ref_img.wcs, img_to_warp, destBBox=bbox)
    warpedExp = copy.deepcopy(warpedExp)

    # Warp PSF
    psf = img_to_warp.getPsf()
    if psf is not None:
        xyTransform = afwGeom.makeWcsPairTransform(img_to_warp.getWcs(),
                                                   ref_img.getWcs())
        warped_psf = measAlg.WarpedPsf(psf, xyTransform)
        warpedExp.setPsf(warped_psf)
    return warpedExp

diff = diff_AlardLupton(deepCoadd_rubin, deepCoadd_mine, warp=True)

plt.figure(figsize=(8,8))
plt.imshow(diff.image.array, origin="lower", cmap="gray",
           vmin=np.percentile(diff.image.array, 0.1),
           vmax=np.percentile(diff.image.array, 99.9))
plt.colorbar()
plt.title("Difference of deep_coadd made by me and the one rubin gives using AlardLupton adaptation.")
plt.show()
2 Likes