105.2. Detect and measure sources¶

105_2_Detect_and_measure_sources

105.2. Detect and measure sources¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.1.1
Last verified to run: 2025-06-25
Repository: github.com/lsst/tutorial-notebooks

Learning objective: How to run source detection, deblending, and measurement pipeline tasks.

LSST data products: visit_image

Packages: lsst.meas.algorithms, lsst.meas.base

Credit: Originally developed by Alex Drlica-Wagner and Imran Hasan in the context of the LSST Stack Club and later refined by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

The lsst.meas.algorithms and lsst.meas.base packages provide the high-level pipeline tasks useful for source detection, deblending, and measurement. Of particular importance for this tutorial are the SourceDetectionTask, the SourceDeblendTask, and the SingleFrameMeasurementTask.

Documentation related to the source detection, deblending, and measurement tasks:

This tutorial provides a brief introduction to running the LSST Science Pipelines source detection, deblending, and measurement tasks.

Related tutorials: The first 100-level tutorial in this series provides a basic introduction to how to use the LSST Science Pipelines; later tutorials in this series demonstrate other image analysis techniques available with the LSST Science Pipelines.

1.1. Import packages¶

From the LSST Science Pipelines, import modules for butler data access (lsst.daf.base and lsst.daf.butler), for image display (lsst.afw), for image geometry (lsst.geom), for table handling (lsst.afw.table), and for various tasks for source detection (lsst.meas.algorithms.detection), deblending (lsst.meas.deblender), and measurement (lsst.meas.base). In addition, import the matplotlib module.

In [1]:

import lsst.daf.base as dafBase
from lsst.daf.butler import Butler
import lsst.afw.image as afw_image
import lsst.afw.display as afw_display
import lsst.afw.table as afwTable
import lsst.geom as geom
from lsst.meas.algorithms.detection import SourceDetectionTask
from lsst.meas.deblender import SourceDeblendTask
from lsst.meas.base import SingleFrameMeasurementTask

1.2. Define parameters and functions¶

Create an instance of the butler with the repository and collection for DP1, and assert that it exists.

In [2]:

butler = Butler("dp1", collections="LSSTComCam/DP1")
assert butler is not None

Define the approximate central coordinates of the Extended Chandra Deep Field South (ECDFS) field (in degrees).

In [3]:

ra = 53.13
dec = -28.10

Set the visualization backend for the LSST Application Framework (AFW) display system to be Firefly, and open frame 1.

In [4]:

afw_display.setDefaultBackend('firefly')
display = afw_display.Display(frame=1)

2. Source detection, deblending, and measurement¶

Learn how to perform source detection, deblending, and measurement on a visit_image.

The LSST Science Pipelines' source detection, deblending, and measurement tasks were imported in the first code cell.

2.1. Select a visit image¶

Using the approximate coordinates of the center of the ECDFS field as declared in Section 1.2, identify the first r-band visit_image in order of time and date of observation (visit.timespan.begin).

First, prepare the query and bind parameters.

In [5]:

query = "band.name = band AND " \
        "visit_detector_region.region OVERLAPS POINT(ra, dec)"
print(query)

bind_params = {"band": 'r', "ra": ra, "dec": dec}
print(bind_params)

band.name = band AND visit_detector_region.region OVERLAPS POINT(ra, dec)
{'band': 'r', 'ra': 53.13, 'dec': -28.1}

Next, perform the query for datasets on the butler. Order the results by visit.timespan.begin. Limit the number of visits returned to 1. Save the dataId.

In [6]:

dataset_refs = butler.query_datasets("visit_image",
                                     where=query,
                                     bind=bind_params,
                                     order_by=["visit.timespan.begin"],
                                     limit=1)
assert len(dataset_refs) == 1

ref = dataset_refs[0]
dataId = ref.dataId
print(dataId)

del dataset_refs

{instrument: 'LSSTComCam', detector: 4, visit: 2024110800247, band: 'r', day_obs: 20241108, physical_filter: 'r_03'}

Retrieve the visit_image for this dataId from the butler.

In [7]:

visit_image = butler.get('visit_image', dataId=dataId)

Clear the DETECTED information from the mask plane, as source detection is re-done below.

In [8]:

visit_image.mask.removeAndClearMaskPlane('DETECTED')

2.2. Source table schema¶

The schema describes the output properties that will be measured for each source. The schema needs to be passed to all of the tasks, each of which will add columns to the schema when it runs.

Create a minimal schema for a source table.

In [9]:

schema = afwTable.SourceTable.makeMinimalSchema()
print(schema)

Schema(
    (Field['L'](name="id", doc="unique ID"), Key<L>(offset=0, nElements=1)),
    (Field['Angle'](name="coord_ra", doc="position in ra/dec"), Key<Angle>(offset=8, nElements=1)),
    (Field['Angle'](name="coord_dec", doc="position in ra/dec"), Key<Angle>(offset=16, nElements=1)),
    (Field['L'](name="parent", doc="unique ID of parent source"), Key<L>(offset=24, nElements=1)),
)

Add coord_raErr and coord_decErr to the minimal schema to avoid warnings from sourceMeasurementTask in Section 2.5, which expects these columns to exist.

In [10]:

raerr = schema.addField("coord_raErr", type="F")
decerr = schema.addField("coord_decErr", type="F")

Create a container which will be used to record metadata about algorithm execution.

In [11]:

algMetadata = dafBase.PropertyList()

2.3. Configure the tasks¶

Each task possesses an associated configuration class. The properties of these configuration classes can be determined from the classes themselves.

Option: uncomment the following line to view the help documentation for the SourceDetectionTask configuration. Replace SourceDetectionTask with another task name to view its help documentation.

In [12]:

# help(SourceDetectionTask.ConfigClass())

Set the basic configuration parameters and instantiate each task in turn.

SourceDetectionTask detects sources.

Set the configuration for the detection threshold to 4 (the default for LSST catalogs is 5) in order to detect sub-threshold (faint) sources. Set the detection threshold type to be stdev, which means the threshold is applied to the image's standard deviation. In other words, this configuration will detect sources at $4\sigma$.

In [13]:

config = SourceDetectionTask.ConfigClass()
config.thresholdValue = 4

config.thresholdType = "stdev"
sourceDetectionTask = SourceDetectionTask(schema=schema, config=config)
del config

SourceDeblendTask deblends overlapping sources into "parent" and "child" sources.

No configuration is needed for this tutorial.

In [14]:

sourceDeblendTask = SourceDeblendTask(schema=schema)

SingleFrameMeasurementTask measures the properties of the deblended sources.

No special configuration is needed for this tutorial.

In [15]:

config = SingleFrameMeasurementTask.ConfigClass()
sourceMeasurementTask = SingleFrameMeasurementTask(schema=schema,
                                                   config=config,
                                                   algMetadata=algMetadata)
del config

Create a SourceTable called tab using the minimal schema. It will hold the output of the analysis in the following section.

In [16]:

tab = afwTable.SourceTable.make(schema)

2.3.1. Re-configure tasks and explore methods¶

The configuration parameters cannot be changed after the task is constructed. To change a configuration parameter, redefine it and then also redefine the task.

For example, to change the detection threshold value used by source detection:

config.thresholdValue = 10
sourceDetectionTask = SourceDetectionTask(schema=schema, config=config)

Option: To explore a task, use the help function on the task or any of its methods.

In [17]:

# help(SourceDetectionTask)

In [18]:

# help(SourceDetectionTask.getPsf)

Option: To see a pop-up list of all methods for a given task, uncomment this line, place the cursor after the period, and press tab. Then re-comment the line. Do not execute the cell.

In [19]:

# SourceDetectionTask.

Now, as schema is no longer needed, delete it.

In [20]:

del schema

2.4. Source detection¶

Run source detection.

In [21]:

result = sourceDetectionTask.run(tab, visit_image)

lsst.sourceDetection INFO: Setting factor for negative detections equal to that for positive detections: 1.000000

lsst.sourceDetection INFO: Detected 4385 positive peaks in 3226 footprints and 125 negative peaks in 114 footprints to 4 +ve and 4 -ve sigma

lsst.sourceDetection INFO: Resubtracting the background after object detection

Print the number of positive peaks detected.

In [22]:

result.numPosPeaks

Out[22]:

Create sources to hold the detected sources from the results.

In [23]:

sources = result.sources
print(len(sources))

With thresholdValue set to 4, roughly 3200 sources are detected.

Had the thresholdValue been set instead to 10, there would have been only roughly 1200 sources detected.

Option: Display sources to see that most of the columns are NaN at first.

In [24]:

# sources

2.5. Source deblending and measurement¶

Run the SourceDeblendTask and SingleFrameMeasurementTask.

In [25]:

sourceDeblendTask.run(visit_image, sources)

lsst.sourceDeblend INFO: Deblending 3226 sources

lsst.sourceDeblend INFO: Deblended: of 3226 sources, 391 were deblended, creating 1084 children, total 4310 sources

In [26]:

sourceMeasurementTask.run(measCat=sources, exposure=visit_image)

lsst.measurement INFO: Measuring 4310 sources (3226 parents, 1084 children)

Ensure that the sources catalog is contiguous (sequential) in memory by making a copy.

In [27]:

sources = sources.copy(True)

Option: Display the sources catalog as an astropy table, which has a nice, human-readable output format.

In [28]:

# sources.asAstropy()

Option: Write the sources to a FITS table, or save the visit_image as a FITS file, in the current directory.

In [29]:

# sources.writeFits("NB_105_4_sources.fits")
# visit_image.writeFits("NB_105_4_visit_image.fits")

2.6. Visualize sources on an image cutout¶

Overplot the detected sources on a cutout image of the visit_image using afwDisplay.

There are several ways to create image cutouts. This tutorial demonstrates the use of the Factory method for a visit_image.

The arguments that must be passed to the Factory method are:

visit_image : the image to make the cutout from
bbox   : the desired bounding box of the cutout
origin : the image pixel origin, local to the cutout array
deep   : if True, copy the data rather than passing by reference

Define the bbox region for a cutout that is centered on pixel coordinate 1700, 2100 and is 400 pixels wide and 400 pixels high.

In [30]:

x, y = 1700, 2100
w, h = 400, 400
xmin, ymin = x - w//2, y - h//2
print(xmin, ymin)
bbox = geom.Box2I()
bbox.include(geom.Point2I(xmin, ymin))
bbox.include(geom.Point2I(xmin + w, ymin + h))
del x, y, w, h, xmin, ymin

1500 1900

Option: Display the bbox information.

In [31]:

# bbox

An alternative statement to create bbox is:

bbox = geom.Box2I(geom.Point2I(xmin, ymin), geom.Extent2I(w, h))

Generate the cutout.

In [32]:

cutout = visit_image.Factory(visit_image, bbox, origin=afw_image.LOCAL, deep=False)

Display the cutout, and set mask transparency to 100%.

In [33]:

display.image(cutout)
display.setMaskTransparency(100)

Mark detected sources with orange circles.

In [34]:

with display.Buffering():
    for s in sources:
        x = s.getX()
        y = s.getY()
        if (x > 1500) & (x < 1900) & (y > 1900) & (y < 2300):
            display.dot('o', x, y, size=20, ctype='orange')
        del x, y

In [ ]: