105.4. Synthetic source injection

105.4. Synthetic source injection¶

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

Learning objective: Use LSST Science Pipelines tools to inject synthetic sources into images.

LSST data products: visit_image

Packages: lsst.source.injection, lsst.daf.butler, lsst.afw.display

Credit: Originally developed by Jeff Carlin in collaboration with Lee Kelvin and the Rubin Community Science Team. Much of the material is based on the documentation for the LSST source_injection package. 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¶

This notebook demonstrates the use of the source_injection package from the LSST Science Pipelines to inject artificial sources (stars and galaxies) into images using the images' measured point-spread function, WCS, and photometric calibration. Sources can be injected into images produced during various stages of Pipelines processing, including visit-level images (e.g., visit_images), any dataset with a datasetType of Exposure (e.g., post_isr_images), and coadd images (e.g., deep_coadds). The driver that both creates and injects synthetic sources into images is based on Galsim, enabling injection of many types of sources, including stars, a variety of parameterized galaxy models, and postage stamp images.

This notebook teaches some basics of using the source_injection package to insert artificial sources into images. Read the extensive documentation for more details and examples. Furthermore, the Galsim documentation provides much more detail on the types of sources that can be injected.

Related tutorials: Detect and measure sources, Forced photometry

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery). Tools from Astropy will be used to read a FITS image an work with Tables.

From the lsst package, import the butler, the lsst.source.injection package to inject synthetic sources into images, and image display functions from the LSST Science Pipelines (pipelines.lsst.io).

In [1]:

import matplotlib.pyplot as plt
from lsst.daf.butler import Butler
import lsst.afw.display as afwDisplay

from astropy.io import fits
from astropy.table import Table, vstack

from lsst.source.injection import generate_injection_catalog
from lsst.source.injection import VisitInjectConfig, VisitInjectTask

1.2. Define parameters and functions¶

Set afwDisplay to use matplotlib.

In [2]:

afwDisplay.setDefaultBackend("matplotlib")

Instantiate the Butler pointing to the DP1 repository and collection, and assert it exists.

In [3]:

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

2. Retrieve an image¶

2.1. Select a visit_image¶

Identify images overlapping a particular sky position and select one to inject sources into.

Query the butler for $r$-band visit_images from detector number 4 that overlap a given RA, Dec position, and select one to use for source injection.

In [4]:

ra = 53.076
dec = -28.110
band = 'r'

In [5]:

query = f"band='{band}' AND detector=4 AND \
          visit_detector_region.region OVERLAPS POINT({ra}, {dec})"

visit_img_refs = butler.query_datasets('visit_image',
                                       where=query)

In [6]:

visit_img = butler.get(visit_img_refs[0])

2.1.1. Extract the image center and dimensions¶

To generate synthetic sources to inject into the image, the coordinates and size of its bounding box are necessary.

Use the image's WCS and bounding box to extract the central (RA, Dec) coordinate and image size (in degrees), and print those to the screen.

In [7]:

cen = visit_img.wcs.pixelToSky(visit_img.getBBox().getCenter())
imsize = visit_img.getBBox().getDimensions()*visit_img.wcs.getPixelScale().asDegrees()

print('Image center in (RA, Dec): ', cen)
print('Size of visit_image in degrees: ', imsize)

Image center in (RA, Dec):  (53.1379468189, -28.0732505763)
Size of visit_image in degrees:  (0.22652, 0.22252)

3. Make a catalog of synthetic sources¶

Having retrieved a visit_image to inject into, now set up a simple synthetic source catalog.

3.1. Make a catalog of galaxies and stars¶

The source_injection package provides tools to create catalogs of synthetic sources. Here, use the generate_injection_catalog function to create the sources to inject.

Note that the "inject_size" is selected to be 0.1 degrees, or slightly smaller than the size of the image as determined above (inject_size is a radius, so it equals 0.2/2 degrees).

In [8]:

inject_size = 0.2/2

3.1.1. Galaxies¶

The simplest form of a galaxy in Galsim is parameterized by a Sersic model: $I(r) = I_e~{\rm exp}\{-b_n [(\frac{r}{r_e})^{1/n}-1]\}$, which defines the shape of the galaxy's light profile as a function of radius (r) in terms of the "Sersic index" (n) and the "half-light radius" ($r_e$). (Note that $n = 1$ corresponds to an exponential profile.)

The above equation results in a circular galaxy. To create elongated (elliptical) Sersic profiles, additionally specify the minor-to-major axis ratio (q) with a rotation angle (beta).

generate_injection_catalog will automatically generate sources with all possible permutations of the parameters you provide. For example, the cell below specifies "number=2" to create 2 synthetic galaxies, but then specifies a single magnitude (mag), 3 values of Sersic index (n), 3 values of axis ratio (q), 2 values of beta, and a single value for half_light_radius. This should result in $2*1*3*3*2*1 = 36$ different combinations of those parameters.

Finally, provide minimum and maximum RA and Dec coordinates in degrees. In this case, generate_injection_catalog will randomly select positions within those limits.

In [9]:

my_injection_catalog_galaxies = generate_injection_catalog(
    ra_lim=[cen.getRa().asDegrees()-inject_size, cen.getRa().asDegrees()+inject_size],
    dec_lim=[cen.getDec().asDegrees()-inject_size, cen.getDec().asDegrees()+inject_size],
    number=2,
    seed='3210',
    source_type="Sersic",
    mag=[15.0],
    n=[1, 2, 4],
    q=[0.9, 0.5, 0.1],
    beta=[31.0, 144.0],
    half_light_radius=[5.0],
)

lsst.source.injection.utils._generate_injection_catalog INFO: Generated an injection catalog containing 36 sources: 18 combinations repeated 2 times.

3.1.2. Stars¶

Using a similar method, create a catalog of stars to be injected. Specify 7 values of magnitude, with 7 instances, which will result in 49 stars.

In [10]:

my_injection_catalog_stars = generate_injection_catalog(
    ra_lim=[cen.getRa().asDegrees()-inject_size, cen.getRa().asDegrees()+inject_size],
    dec_lim=[cen.getDec().asDegrees()-inject_size, cen.getDec().asDegrees()+inject_size],
    number=7,
    seed='432',
    mag=[19.0, 19.5, 20.0, 20.5, 21.0, 21.5, 22.0],
    source_type="Star",
)

lsst.source.injection.utils._generate_injection_catalog INFO: Generated an injection catalog containing 49 sources: 7 combinations repeated 7 times.

3.1.3. Combine catalogs¶

Before combining the catalogs, offset the injection_id values for the star catalog by 1000 to avoid overlapping with the galaxy ids.

Use the vstack method from Astropy's Table class to combine the galaxy and star catalogs.

In [11]:

my_injection_catalog_stars['injection_id'] += 1000

inject_cat = vstack([my_injection_catalog_galaxies, my_injection_catalog_stars])

Inspect the catalog.

In [12]:

inject_cat

Out[12]:

Table length=85

injection_id	ra	dec	source_type	mag	n	q	beta	half_light_radius
int64	float64	float64	str6	float64	int64	float64	float64	float64
0	53.228239896616735	-28.106616457094464	Sersic	15.0	1	0.9	31.0	5.0
1	53.109489896616736	-27.98315966697101	Sersic	15.0	1	0.9	31.0	5.0
10	53.046989896616736	-28.01032016079817	Sersic	15.0	1	0.9	144.0	5.0
11	53.040739896616735	-28.15846830894632	Sersic	15.0	1	0.9	144.0	5.0
20	53.18448989661673	-28.11649300030434	Sersic	15.0	1	0.5	31.0	5.0
21	53.12823989661673	-28.091801642279652	Sersic	15.0	1	0.5	31.0	5.0
30	53.09073989661673	-28.00291275339076	Sersic	15.0	1	0.5	144.0	5.0
31	53.15948989661673	-28.00538188919323	Sersic	15.0	1	0.5	144.0	5.0
40	53.05323989661673	-28.06217201265002	Sersic	15.0	1	0.1	31.0	5.0
...	...	...	...	...	...	...	...	...
1054	53.10965263826026	-28.004947737655698	Star	21.5	--	--	--	--
1055	53.21277763826026	-27.99260205864335	Star	21.5	--	--	--	--
1056	53.16277763826026	-28.08149094753224	Star	21.5	--	--	--	--
1060	53.22527763826026	-28.148157614198908	Star	22.0	--	--	--	--
1061	53.19090263826026	-28.11852798456928	Star	22.0	--	--	--	--
1062	53.11277763826026	-28.05926872531002	Star	22.0	--	--	--	--
1063	53.14715263826026	-28.101244033951996	Star	22.0	--	--	--	--
1064	53.17527763826026	-28.037046503087797	Star	22.0	--	--	--	--
1065	53.22840263826026	-28.167910700618663	Star	22.0	--	--	--	--
1066	53.04715263826026	-28.14568847839644	Star	22.0	--	--	--	--

3.2. Add a postage stamp image entry¶

In addition to parameterized sources of many types, a postage-stamp image (for example, an output image from a simulation) can also be injected. This can be any image; the only requirement is that it must be in FITS format and have a valid WCS (see a fun example here). This notebook demonstrates injection of an SDSS g-band image of the galaxy PGC 038749 downloaded from the NASA/IPAC Extragalactic Database (NED).

Load the image, then display it to see what it depicts:

In [13]:

stamp_path = '/rubin/cst_repos/tutorial-notebooks-data/data/'
stamp_img_hdu = fits.open(stamp_path + 'PGC_038749_I_g_bbl2011.fits')

In [14]:

fig = plt.figure()
plt.subplot()
im = plt.imshow(stamp_img_hdu[0].data, cmap='gray', vmin=-20.0, vmax=100, origin='lower')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Figure 1: The SDSS g-band image of PGC 038749 displayed in grayscale, with axes in pixel coordinates.

Add a postage stamp entry into the injection catalog, using vstack to add an astropy Table containing a single row. Specify its position, magnitude, and "source_type = Stamp," in addition to the filename of the image to inject. (It is OK to leave unnecessary columns empty.)

In [15]:

my_injection_catalog_stamp = Table(
    {
        'injection_id': [9999],
        'ra': [cen.getRa().asDegrees()+0.01],
        'dec': [cen.getDec().asDegrees()-0.01],
        'source_type': ['Stamp'],
        'mag': [14.8],
        'stamp': [stamp_path + 'PGC_038749_I_g_bbl2011.fits'],
    }
)

inject_cat = vstack([inject_cat, my_injection_catalog_stamp])

Option to view the injection catalog:

In [16]:

# inject_cat

4. Inject sources into an image¶

In Section 3, a catalog specifying the synthetic sources to inject was created.

Now run the task from source_injection that injects sources into the visit_image that was retrieved earlier.

Pass the point spread function (PSF), photometric calibration object, and the WCS (World Coordinate System) associated with the visit_image that was retrieved earlier to the injection task so that sources can be injected using the properties of the image itself.

4.1. Instantiate the injection classes¶

First, instantiate the VisitInjectConfig class. The VisitInjectConfig class is where configuration of the injection task occurs, allowing for modifications to be made to how the task operates.

Then instantiate the VisitInjectTask, using inject_config as the configuration argument.

NOTE: For injections into other dataset types, use the appropriate option from the following list:

from lsst.source.injection import CoaddInjectConfig, CoaddInjectTask

from lsst.source.injection import ExposureInjectConfig, ExposureInjectTask

from lsst.source.injection import VisitInjectConfig, VisitInjectTask

Both ExposureInject and VisitInject accept dimensions of ("instrument", "visit", "detector"). They differ in that ExposureInject takes a post_isr_image exposure (i.e., a raw image that has had "instrument signature removal," or ISR, applied, and no further processing) as an input, while VisitInject expects to operate on a visit_image (an image that has been astrometrically and photometrically calibrated to an external source).

In [17]:

inject_config = VisitInjectConfig()

inject_task = VisitInjectTask(config=inject_config)

4.2. Run the source_injection task¶

Execute the run method of the inject_task.

The run method needs the following inputs:

• the input injection catalog
• the input visit_image
• the PSF of the input exposure
• the WCS information
• the photometric calibration information

The PSF, WCS, and photo_calib inputs are associated with the visit_image that was already retrieved.

The inject task provides two outputs:

• the output exposure with sources injected
• the output source injection catalog

The output source injection catalog is identical to the input, but with two additional columns (x and y) denoting the pixel coordinates of these sources. Note that this catalog is NOT the science catalog containing the full suite of LSST Science Pipelines outputs. To get that, this source injected image will need to be processed by additional Science Pipelines tasks.

Note: it is best to use a clone of the input visit_image. This is because the visit_image that is input for source injection will be edited in-place. Inputting a clone makes it possible to compare the output image to the original visit_image later in this notebook.

Run the source injection task and extract the "output_exposure" and "output_catalog":

In [18]:

injected_output = inject_task.run(
    injection_catalogs=inject_cat,
    input_exposure=visit_img.clone(),
    psf=visit_img.psf,
    photo_calib=visit_img.photoCalib,
    wcs=visit_img.wcs,
)
injected_exposure = injected_output.output_exposure
injected_catalog = injected_output.output_catalog

lsst.visitInjectTask INFO: Retrieved 86 injection sources from 86 HTM trixels.

lsst.visitInjectTask INFO: Catalog cleaning removed 0 of 86 sources; 86 remaining for catalog checking.

lsst.visitInjectTask INFO: Catalog checking flagged 0 of 86 sources; 86 remaining for source generation.

lsst.visitInjectTask INFO: Adding INJECTED and INJECTED_CORE mask planes to the exposure.

lsst.visitInjectTask INFO: Generating 86 injection sources consisting of 3 unique types: Sersic(36), DeltaFunction(49), Stamp(1).

lsst.visitInjectTask INFO: Injected 86 of 86 potential sources. 0 sources flagged and skipped.

4.3. Compare the images before and after injection¶

Display the images in separate panels with matplotlib, and compare them to confirm that the sources have been injected.

Optionally, uncomment the commented-out "display1.mtv" cell to zoom in on the injected stamp image.

In [19]:

plot_injected_img = injected_exposure.clone()

fig, ax = plt.subplots(1, 2, figsize=(10, 6), dpi=150)

plt.sca(ax[0])
display0 = afwDisplay.Display(frame=fig)
display0.scale('linear', 'zscale')
display0.mtv(visit_img.image)
plt.title('calexp image')

plt.sca(ax[1])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(plot_injected_img.image)
# To zoom on the PGC 038749 stamp:
# display1.mtv(plot_injected_img.image[1700:2200, 1950:2450])
plt.title('injected_calexp image')

plt.tight_layout()
plt.show()

In [ ]: