Mask for visit_image cutout is not aligned with image

Support Data Products
,
#1

I am working with some visit_image’s and getting cutouts around my object of interest, I noticed in one of the cutouts that the associated mask was wrong (clearly from a totally different part of the image) meanwhile other cutout masks lined up perfectly. I need the masks to flag bad pixels etc, so I might have to just take the full image then manually crop on the full visit_image.mask.arraywhich does seem to be lined up correctly with the visit_image.image.array. I’ve attached a script to reproduce the issue, sorry its not really a “minimal example” but I am still very new so I didn’t want to mess up and lose the one I found that’s misaligned.

from lsst.rsp import get_tap_service
service = get_tap_service("tap")

from lsst.daf.butler import Butler
butler = Butler("dp1", collections="LSSTComCam/DP1")
assert butler is not None

import lsst.geom as geom

import matplotlib.pyplot as plt
import numpy as np
from astropy import units

ra_sn, dec_sn = 58.335054, -48.750303
query = (
    "SELECT ra, dec, diaObjectId "
    "FROM dp1.DiaObject "
    f"WHERE CONTAINS (POINT('ICRS', ra, dec), CIRCLE('ICRS',{ra_sn}, {dec_sn}, {5/60/60})) = 1 "
)

job = service.submit_job(query)
job.run()
job.wait(phases=["COMPLETED", "ERROR"])
print("Job phase is", job.phase)
job.raise_if_error()
assert job.phase == "COMPLETED"

df_sn = job.fetch_result().to_table()
print(f"\nFound {len(df_sn)} source(s)")
diaobjectid = df_sn["diaObjectId"][0]
print(df_sn)

query = (
    "SELECT fsodo.coord_ra, fsodo.coord_dec, fsodo.diaObjectId, fsodo.visit, fsodo.detector, fsodo.band, vis.skyRotation, "
    "fsodo.tract, fsodo.patch, fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, fsodo.psfFlux, fsodo.psfFluxErr, vis.expMidptMJD "
    "FROM dp1.ForcedSourceOnDiaObject as fsodo "
    "JOIN dp1.Visit as vis ON vis.visit = fsodo.visit "
    f"WHERE fsodo.diaObjectId = {diaobjectid}"
)

job = service.submit_job(query)
job.run()
job.wait(phases=["COMPLETED", "ERROR"])
print("Job phase is", job.phase)
job.raise_if_error()
assert job.phase == "COMPLETED"

df_exposure = job.fetch_result().to_table()
print(f"\nFound {len(df_exposure)} sources")
print(df_exposure)

i = -1 # set to 5 and the mask is aligned correctly!

band = df_exposure["band"][i]
visit = df_exposure["visit"][i]
detector = df_exposure["detector"][i]
tract = df_exposure["tract"][i]
patch = df_exposure["patch"][i]
mjd = df_exposure["expMidptMJD"][i]
dataId = {"band": band, "visit": visit, "detector": detector, "tract": tract, "patch": patch}

# Get calexp, template and sources from butler
visit_image = butler.get("visit_image", dataId=dataId)

# Full image, the mask appears visually aligned
plt.imshow(np.log(visit_image.image.array), origin="lower") 
plt.show()
plt.imshow(visit_image.mask.array, origin="lower") 
plt.show()

# Cutout, the mask is misaligned
spherePoint = geom.SpherePoint(ra_sn*geom.degrees, dec_sn*geom.degrees)
Radius = 0.01 * units.deg
cutout = visit_image.getCutout(spherePoint, geom.Extent2I(256))

plt.imshow(cutout.image.array, origin="lower") 
plt.show()
plt.imshow(cutout.mask.array, origin="lower") 
plt.show()

image
image389×756 63.1 KB

#2

Ok, slight correction. I noticed that some of the mask values are gigantic and were washing out other mask values. So I made a figure that just identifies every non-zero pixel in the mask. Now it looks like the mask is in fact aligned (that’s good!) but it is flagging things that don’t exist also. The mask in my first comment appears to be trying to identify a bright star when in fact there are no bright stars in the image.

image

#3

Note that those mask objects are bitmasks, not booleans; one bit can mean something entirely different from another. My best guess is that you’re seeing the CROSSTALK plane, which says where we applied a crosstalk correction for a star somewhere else on the same detector.

There are some methods on the Mask object for inspecting the list of mask planes and the bits their correspond to, and I think there’s a tutorial notebook in the works on that subject (cc @plazas). But the best way to view them IMO is to use the portal (i.e. Firefly), e.g via lsst.afw.display.

1 Like
#4

Ah I see, is there some documentation that describes what each of the flags means? I’m not sure if this crosstalk flag means I should try to ignore that data or not. I am working on an analysis that will involve fitting the pixel data and so I need to know which bits indicate a bad pixel value. The portal might let me see the mask for diagnostic purposes, but I’ll need an automated way to identify which bits are “bad” and build them into a mask.

#5

Hi @ConnorStone, so far we have one notebook tutorial that provides an introduction to the mask plane: 202.2. Visit images. That link will take you to an html file of the executed notebook; to open the executable ipynb file in the Notebook Aspect of the RSP, click on “Tutorials” and select it from the drop-down menu.

I don’t think that all pixels flagged as potentially affected by crosstalk always need to be ignored, no.

As Jim mentions, a new tutorial that provides a deeper discussion of the mask planes, and provides guidance for their use, is currently in development.

Take a look at tutorial notebook 202.2, and if that is sufficient for your current needs, mark this reply post as the solution? And if it is not – that’s ok, we can keep the conversation on masks going here in this thread.

#6

This is a function we have used in the past for extracting selected mask planes.
Could be useful so I put it here.

import numpy as np
from lsst.afw.image import ExposureF


def get_mask(image: ExposureF, mask_names: str | list[str]) -> np.ndarray:
    """Extract a boolean mask array from an LSST image for specified mask planes.
    
    Args:
        image: LSST image object with mask information
        mask_names: Single mask plane name or list of mask plane names
        
    Returns:
        Boolean array. True indicates pixels flagged in any of the specified mask planes
    """
    mask = image.getMask()
    if isinstance(mask_names, str):
        mask_names = [mask_names]

    mask_array = mask.getArray()
    out = np.zeros_like(mask_array)
    for mask_name in mask_names:
        target_bit = mask.getMaskPlane(mask_name)
        out |= (mask_array & (2 ** target_bit)) != 0
    return out

You can find mask plane names in Section 3.3 of the notebook Melissa linked in her last message.

Bests!

2 Likes
#7

@MelissaGraham thanks for the clarification, I think I can use this information for now, and I’ll be very interested when the more in depth tutorial comes out!

Also thanks @deppep that will help me pick out what I want from the bitmask!

2 Likes