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-09-02
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).
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¶
Define parameters to use colorblind-friendly colors with matplotlib.
plt.style.use('seaborn-v0_8-colorblind')
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
Set afwDisplay to use matplotlib.
afwDisplay.setDefaultBackend("matplotlib")
Instantiate the Butler pointing to the DP1 repository and collection, and assert it exists.
butler = Butler("dp1", collections="LSSTComCam/DP1")
assert butler is not None
ra = 53.076
dec = -28.110
band = 'r'
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,
order_by='visit.timespan.begin')
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.
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).
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.
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.
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.
my_injection_catalog_stars['injection_id'] += 1000
inject_cat = vstack([my_injection_catalog_galaxies, my_injection_catalog_stars])
Inspect the catalog.
inject_cat
| 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 |
| 41 | 53.203239896616736 | -27.995505345983354 | Sersic | 15.0 | 1 | 0.1 | 31.0 | 5.0 |
| 50 | 53.171989896616736 | -27.988097938575944 | Sersic | 15.0 | 1 | 0.1 | 144.0 | 5.0 |
| 51 | 53.15323989661673 | -28.151060901538912 | Sersic | 15.0 | 1 | 0.1 | 144.0 | 5.0 |
| 60 | 53.218864896616736 | -28.160937444748786 | Sersic | 15.0 | 2 | 0.9 | 31.0 | 5.0 |
| 61 | 53.05948989661673 | -28.138715222526564 | Sersic | 15.0 | 2 | 0.9 | 31.0 | 5.0 |
| 70 | 53.196989896616735 | -28.143653494131502 | Sersic | 15.0 | 2 | 0.9 | 144.0 | 5.0 |
| 71 | 53.071989896616735 | -28.12143127190928 | Sersic | 15.0 | 2 | 0.9 | 144.0 | 5.0 |
| 80 | 53.21573989661674 | -27.980690531168538 | Sersic | 15.0 | 2 | 0.5 | 31.0 | 5.0 |
| 81 | 53.118864896616735 | -28.072048555859897 | Sersic | 15.0 | 2 | 0.5 | 31.0 | 5.0 |
| 90 | 53.078239896616736 | -27.97328312376113 | Sersic | 15.0 | 2 | 0.5 | 144.0 | 5.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1042 | 53.06590263826026 | -28.153095885803847 | Star | 21.0 | -- | -- | -- | -- |
| 1043 | 53.07215263826026 | -27.990132922840882 | Star | 21.0 | -- | -- | -- | -- |
| 1044 | 53.13465263826026 | -28.049392182100142 | Star | 21.0 | -- | -- | -- | -- |
| 1045 | 53.10340263826026 | -28.108651441359402 | Star | 21.0 | -- | -- | -- | -- |
| 1046 | 53.05340263826026 | -27.997540330248288 | Star | 21.0 | -- | -- | -- | -- |
| 1050 | 53.15027763826026 | -28.11112057716187 | Star | 21.5 | -- | -- | -- | -- |
| 1051 | 53.04090263826026 | -28.074083540124835 | Star | 21.5 | -- | -- | -- | -- |
| 1052 | 53.16590263826026 | -28.05186131790261 | Star | 21.5 | -- | -- | -- | -- |
| 1053 | 53.09090263826026 | -27.985194651235943 | Star | 21.5 | -- | -- | -- | -- |
| 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:
stamp_path = '/rubin/cst_repos/tutorial-notebooks-data/data/'
stamp_img_hdu = fits.open(stamp_path + 'PGC_038749_I_g_bbl2011.fits')
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.)
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:
# 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).
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 thevisit_imagethat is input for source injection will be edited in-place. Inputting a clone makes it possible to compare the output image to the originalvisit_imagelater in this notebook.
Run the source injection task and extract the "output_exposure" and "output_catalog":
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.
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.scale('linear', min=-20, max=150)
display0.mtv(visit_img.image)
plt.title('calexp image')
plt.sca(ax[1])
display1 = afwDisplay.Display(frame=fig)
# display1.scale('linear', 'zscale')
display1.scale('linear', min=-20, max=150)
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()
