202.2. Visit images

202.2. Visit images¶

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-21
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To understand the format and metadata of a visit_image.

LSST data products: visit_images

Packages: lsst.daf.butler, lsst.rsp, lsst.afw.display

Credit: Originally developed 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¶

A visit_image is a fully processed, calibrated, science-ready astronomical image from a single detector (a single CCD), obtained at a given time in a single band.

Metadata for visits (observations) and visit_images are available in the Visit and CcdVisit TAP tables, respectively.

• butler dataset type: visit_image
• number of images: 15972

Related tutorials: The 100-level tutorials on how to use the Butler, TAP, SIA, Firefly, and the cutout service demonstrate how to query for and retrieve images, how to create small image cutouts ("stamps"), and how to manipulate the Firefly display interface. The 200-level tutorials on the photometric and astrometric calibration, and the point spread function, describe these aspects of the images in more detail.

1.1. Import packages¶

From the LSST Science Pipelines import the packages for the Butler, the TAP service, 2-dimensional sky geometry, and for image display. Also import standard astropy and numpy packages.

In [1]:

from lsst.daf.butler import Butler, Timespan
from lsst.rsp import get_tap_service
import lsst.geom
import lsst.afw.display as afwDisplay
from astropy.time import Time
import numpy as np

1.2. Define parameters and functions¶

Instantiate the Butler.

In [2]:

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

Create an instance of the TAP service, and assert that it exists.

In [3]:

service = get_tap_service("tap")
assert service is not None

Set afwDisplay to use Firefly and open the Firefly display tab.

In [4]:

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

2. Data access¶

It is recommended to use the Rubin data Butler to retrieve image data within Jupyter Notebooks, but they are also available via the TAP and SIA services (Table Access Protocol and Simple Image Access).

Recommended access method: Butler

2.1. Butler¶

Show that the dimensions for the visit_image dataset type are the filter, instrument, observation date, visit, and detector, and that the returned image will be of type ExposureF. The LSST Science Pipelines tools for image manipulation and visualization work best with the ExposureF format. This format includes pixel data and metadata.

In [5]:

butler.get_dataset_type('visit_image')

Out[5]:

DatasetType('visit_image', {band, instrument, day_obs, detector, physical_filter, visit}, ExposureF)

Show that only instrument, detector, and visit are required to uniquely identify a visit_image.

In [6]:

butler.get_dataset_type('visit_image').dimensions.required

Out[6]:

{instrument, detector, visit}

2.1.1. Demo query¶

Query for visit_images that overlap coordinates RA, Dec near the center of the ECDFS field, were obtained with the r-band filter, and a within timespan near the start of the LSSTComCam campaign.

In [7]:

ra = 53.076
dec = -28.110
band = "r"
time1 = Time("2024-11-09T00:00:00.0", format="isot", scale="tai")
time2 = Time(60624.0, format="mjd", scale="tai")
timespan = Timespan(time1, time2)
del time1, time2

In [8]:

dataset_refs = butler.query_datasets("visit_image",
                                     where="band.name = :band AND \
                                     visit.timespan OVERLAPS :timespan AND \
                                     visit_detector_region.region OVERLAPS POINT(:ra, :dec)",
                                     bind={"band": band, "timespan": timespan,
                                           "ra": ra, "dec": dec},
                                     order_by=["visit.timespan.begin"])
print(len(dataset_refs))

18

In [9]:

del ra, dec, band, timespan

2.1.2. Get an image¶

The visit_image can be retrieved from the Butler using the (a) the dataset_refs; (b) the dataId; or (c) the visit and detector numbers.

(a) Use the dataset reference¶

For the first dataset reference returned by the query, get the corresponding visit_image.

In [10]:

ref = dataset_refs[0]
visit_image = butler.get(ref)

See that it is of type ExposureF.

In [11]:

visit_image

Out[11]:

<lsst.afw.image._exposure.ExposureF at 0x7d2286d349f0>

Display the image in the Firefly window. Set the mask transparency to 100 (fully transparent, i.e., do not show the mask).

In [12]:

afw_display.mtv(visit_image)
afw_display.setMaskTransparency(100)

Clean up.

In [13]:

del ref, visit_image

(b) Use the dataId¶

For the first dataset reference returned by the query, get the dataId.

In [14]:

ref = dataset_refs[0]
use_dataId = ref.dataId
use_dataId

Out[14]:

{instrument: 'LSSTComCam', detector: 0, visit: 2024110800246, band: 'r', day_obs: 20241108, physical_filter: 'r_03'}

Use the defined dataId to retrieve the corresponding visit_image.

In [15]:

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

Option to display the visit_image.

In [16]:

# afw_display.mtv(visit_image)

Clean up.

In [17]:

del ref, use_dataId, visit_image

(c) Use the visitId and detector¶

For the third dataset reference returned by the query, get the dataId and extract the visit and detector numbers.

In [18]:

ref = dataset_refs[2]
dataId = ref.dataId
use_visit = dataId.get('visit')
use_detector = dataId.get('detector')
print(use_visit, use_detector)

2024110800250 4

Use the visit and detector to get the corresponding visit_image.

Even if the Butler is not used to query for the images (e.g., if the SIA or TAP service is used to find images matching a set of constraints), images can be retrieved from the Butler by passing the visit and detector.

In [19]:

visit_image = butler.get('visit_image', visit=use_visit,
                         detector=use_detector)

Option to display the visit_image.

In [20]:

# afw_display.mtv(visit_image)

Clean up.

In [21]:

del ref, dataId, use_visit, use_detector, visit_image

2.2. SIA (Simple Image Access)¶

The Simple Image Access (SIA) service provides a standardized model for image metadata, and the capability to query and retrieve image datasets. Use of SIA, and data retrieval via the access_url in the SIA results table, is demonstrated in the 100-level tutorials.

To return only visit_images from a SIA query, set:

• calib_level = 2
• dpsubtype = 'lsst.visit_image'

Since 'lsst.visit_image' is the only data product subtype included for calibration level 2, SIA queries only need to specifify one or the other.

See Section 2.1.2. for how to retrieve images via the butler using the lsst_visit and lsst_detector columns in the SIA results table.

2.3. TAP (Table Access Protocol)¶

The Table Access Protocol (TAP) service provides a standardized model for accessing image metadata and catalog data with the Astronomical Data Query Language (ADQL). Use of TAP, and data retrieval via the access_url in the TAP results table, is demonstrated in the 100-level tutorials.

2.3.1. ObsCore table¶

Image metadata is stored in the ObsCore table.

To return only visit_images from a TAP query, set:

• Calibration level 2.
• Dataproduct subtype: lsst.visit_image

See Section 2.1.2. for how to retrieve images via the butler using the lsst_visit and lsst_detector columns in the TAP results table.

2.3.2. CcdVisit table¶

Image metadata for individual visit_images is also stored in the CcdVisit table, which is access via TAP.

Unlike the SIA model and the ObsCore table, the CcdVisit table includes measured image properties such as the size of the PSF, the zeropoint, and the sky background.

Query the CcdVisit table for i-band detector images whose corners contain the target coordinates. Include a constraint that the measured seeing be less than 0.8 arcseconds. Order the results by expMidptMJD (the MJD at the midpoint of the exposure).

In [22]:

query = "SELECT visitId, detector, expMidptMJD, seeing, band " \
        "FROM dp1.CcdVisit "\
        "WHERE CONTAINS(POINT('ICRS', 53.13, -28.10), " \
        "POLYGON('ICRS', llcra,llcdec, ulcra,ulcdec, urcra,urcdec, lrcra,lrcdec)) = 1 " \
        "AND band = 'i' " \
        "AND seeing < 0.8 " \
        "ORDER BY expMidptMJD ASC"

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

Job phase is COMPLETED

In [23]:

assert job.phase == 'COMPLETED'
results = job.fetch_result().to_table()
assert len(results) == 3
results

Out[23]:

Table length=3

visitId	detector	expMidptMJD	seeing	band
		d	arcsec
int64	int64	float64	float64	object
2024120500065	4	60650.20580252321	0.7925209258075652	i
2024120500070	5	60650.208053304406	0.7913529349439603	i
2024120500093	5	60650.220256660876	0.7788051213021917	i

The method for retrieving and displaying the first visit_image in the list above is the same as is demonstrated in Section 2.1.2.

visit_image = butler.get('visit_image',
                         visit=2024120500065,
                         detector=4)
afw_display.mtv(visit_image)

In [24]:

job.delete
del query, results

3. Pixel data¶

The visit_images have three planes of pixel data: image, variance, and mask.

Retrieve the second visit_image from the Butler query in Section 2.1.1.

In [25]:

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

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

In [26]:

visit_image = butler.get(ref)

Display the visit_image in frame 1, and set the mask tranparency to 0 (not transparent; to show the mask).

In [27]:

afw_display.mtv(visit_image)
afw_display.setMaskTransparency(0)

3.1. Image plane¶

Sky pixel data in units of nJy (nanojanskys).

The background has been subtracted so values can be negative.

Extract the image data from the visit_image, and convert it to an array.

In [28]:

image = visit_image.getImage()
image_array = image.array

Option to show the image and the image_array.

In [29]:

# image

In [30]:

# image_array

Print statistics on the image_array.

In [31]:

print(image_array.min())
print(image_array.mean())
print(image_array.max())

-158.54198
6.979102
67914.836

Display the image in Firefly frame 2.

In [32]:

afw_display = afwDisplay.Display(frame=2)
afw_display.mtv(image)

Mouse-over display frame 2. As the image contains only pixel data, the image is displayed without a WCS: the coordinates are in pixels, not RA, Dec (as in frame 1).

3.2. Variance plane¶

Uncertainty (noise) in the sky pixel data in units of nJy^2 (nanoJanskys squared).

Values are always positive.

Extract the variance and variance_array.

In [33]:

variance = visit_image.getVariance()
variance_array = variance.array

Print statistics on the variance_array.

In [34]:

print(variance_array.min())
print(variance_array.mean())
print(variance_array.max())

163.38919
301.5471
61834.23

Display the variance.

In [35]:

afw_display = afwDisplay.Display(frame=3)
afw_display.mtv(variance)

3.3. Mask plane¶

An integer bitmask of representative flag values that indicate processing status or issues, similar to the SDSS bitmasks.

Extract the mask.

In [36]:

mask = visit_image.getMask()
mask_array = mask.array

Mask keys (names) are defined by the mask and interpreted as colors by afw_display.

Print the list of mask keys and color.

In [37]:

for mask_key, bit in mask.getMaskPlaneDict().items():
    print('{}: {}'.format(mask_key,
                          afw_display.getMaskPlaneColor(mask_key)))

BAD: red
CLIPPED: blue
CR: magenta
CROSSTALK: cyan
DETECTED: blue
DETECTED_NEGATIVE: cyan
EDGE: yellow
INEXACT_PSF: magenta
INJECTED: yellow
INJECTED_TEMPLATE: orange
INTRP: green
ITL_DIP: red
NOT_DEBLENDED: green
NO_DATA: orange
REJECTED: blue
SAT: green
SAT_TEMPLATE: cyan
SENSOR_EDGE: magenta
STREAK: green
SUSPECT: yellow
UNMASKEDNAN: yellow
VIGNETTED: red

Display the mask in Firefly frame 4.

In [38]:

afw_display = afwDisplay.Display(frame=4)
afw_display.mtv(mask)

The colors are different because the Firefly display interprets the pixel values as flux values (mouse-over and see).

Set the mask transparency to opaque to see the same colors as in frame 1.

In [39]:

afw_display.setMaskTransparency(0)

Set all mask key names to transparent except the "DETECTED" mask bit.

In [40]:

afw_display.setMaskTransparency(100)
afw_display.setMaskTransparency(0, 'DETECTED')

Option to print a list of the unique mask array values, the number of pixels with that value, and the overlapping masks for pixels with that value.

In [41]:

# values, counts = np.unique(mask_array, return_counts=True)
# for value, count in zip(values, counts):
#     print(value, '('+str(count)+')', mask.interpret(value))
# del values, counts

Clean up.

In [42]:

del image, image_array, variance, variance_array, mask, mask_array
del ref, dataset_refs

4. Metadata and attributes¶

4.1. Header¶

The informational equivalent of a FITS image header is the metadata.

Get the visit_image metadata.

In [43]:

metadata = visit_image.getMetadata()

Option to display the long list of metadata.

In [44]:

# metadata

Convert the metadata to a python dictionary.

In [45]:

md_dict = metadata.toDict()

Show any metadata key that contains the string "UNIT".

In [46]:

for key in md_dict.keys():
    if key.find('UNIT') >= 0:
        print(key)

LSST ISR UNITS
LSST ISR OVERSCAN SERIAL UNITS
LSST ISR READNOISE UNITS
LSST ISR OVERSCAN PARALLEL UNITS
BUNIT

Print the "BUNIT" from the dictionary (the flux units).

Warning: As of June 2025, there was an error in the header that reported the visit_image flux units as 'adu' when they are 'nJy'. If 'adu' is returned below, then DM-51270 is still in progress.

In [47]:

md_dict['BUNIT']

Out[47]:

'nJy'

Clean up.

In [48]:

del metadata, md_dict

4.2. Visit information¶

Information about the observation associated with this visit_image, such as pointing, airmass, and weather.

Get the visit info.

In [49]:

info = visit_image.getInfo().getVisitInfo()

In [50]:

print(info.getBoresightAirmass())
print(info.getBoresightParAngle())
print(info.getDate())
print(info.getExposureTime())

1.03979582779195
1.78014 rad
DateTime("2024-11-09T06:14:14.959995994", TAI)
30.0

Get weather information.

In [51]:

weather = info.getWeather()

In [52]:

print(weather.getAirPressure())
print(weather.getAirTemperature())
print(weather.getHumidity())

74135.0
14.375
17.1749992370605

In [53]:

del info, weather

4.3. World Coordinate System (WCS)¶

The reference system that maps pixel coordinate (x, y) to sky coordinate (RA, Dec).

Get the WCS and show it.

In [54]:

wcs = visit_image.getWcs()
wcs

Out[54]:

Non-standard SkyWcs (Frames: PIXELS, NORMEDPIXELS, INTERMEDIATE0, IWC, SKY):
Sky Origin: (53.1413782856, -28.1312112069)
Pixel Origin: (1031.71, 2283.96)
Pixel Scale: 0.200266 arcsec/pixel

Convert pixel coordinates to sky coordinates, and vice versa.

In [55]:

coord = lsst.geom.SpherePoint(53.13*lsst.geom.degrees, -28.10*lsst.geom.degrees)
xy = wcs.skyToPixel(coord)

In [56]:

print(xy)

(1538.7, 1983.4)

In [57]:

xy2 = lsst.geom.Point2D(1021.41, 1962.7)
coord2 = wcs.pixelToSky(xy2)

In [58]:

print(coord2)

(53.1214755393, -28.1277984229)

In [59]:

del xy, xy2, coord, coord2

4.4. Bounding box (corners)¶

The bounding box defines the extent (corners) of the image.

Get the bounding box.

In [60]:

bbox = visit_image.getBBox()

In [61]:

bbox

Out[61]:

Box2I(corner=Point2I(0, 0), dimensions=Extent2I(4072, 4000))

Print the visit_image corners in pixels and sky coordinates.

In [62]:

corners = (lsst.geom.Point2D(bbox.beginX, bbox.beginY),
           lsst.geom.Point2D(bbox.beginX, bbox.endY),
           lsst.geom.Point2D(bbox.endX, bbox.endY),
           lsst.geom.Point2D(bbox.endX, bbox.beginY))
for corner in corners:
    print(corner, wcs.pixelToSky(corner))

(0, 0) (52.9863939585, -28.1588635486)
(0, 4000) (53.2325427696, -28.2083585147)
(4072, 4000) (53.2893308170, -27.9874270049)
(4072, 0) (53.0436604449, -27.9380370324)

Print the area of the bounding box.

In [63]:

print(bbox.area, 'square pixels')
scale = wcs.getPixelScale().asArcminutes()
area = bbox.area * scale**2
print(np.round(area, 2), 'square arcminutes')

16288000 square pixels
181.46 square arcminutes

Check that the first set of pixel values are within the bounding box, and that the second set is not.

In [64]:

assert bbox.contains(100, 100) is True
print(bbox.contains(100, 100))
assert bbox.contains(10000, 10000) is False
print(bbox.contains(10000, 10000))

True
False

Check that these sky coordinates are not within the bounding box.

In [65]:

coord = lsst.geom.SpherePoint(63.0*lsst.geom.degrees, -38.0*lsst.geom.degrees)
xy = wcs.skyToPixel(coord)
assert bbox.contains(xy.getX(), xy.getY()) is False
print(bbox.contains(xy.getX(), xy.getY()))

False

Clean up.

In [66]:

del wcs, bbox, corners, scale, area, coord, xy

4.5. Photometric calibration¶

The conversion between instrumental flux and calibrated flux (nanoJanskys, AB magnitudes).

Calibrated fluxes for all of the sources detected in the visit_image are already listed in the Source table. But if a source's flux is remeasured from the visit_image, then use the photometric calibraton.

Get the photometric calibration.

In [67]:

photcalib = visit_image.getPhotoCalib()

In [68]:

photcalib

Out[68]:

PhotoCalib(spatially constant with mean: 1 error: 0.000388662)

The flux units of the visit_image are nanoJanskys. The flux can be converted to AB magnitude with $m = -2.5 \log(f) + 31.4$.

Use the instFluxToMagnitude method associated with the photoCalib to convert an instrumental flux of 10000 nJy to AB magnitude.

In [69]:

print(photcalib.instFluxToMagnitude(10000))

21.4

In [70]:

del photcalib

4.6. Point spread function (PSF)¶

The PSF is the empirical estimate of the shape of a point source as a function of position in the image. The PSF encapsulates all of the atmospheric and optical effects that blur point sources into their measured shapes.

Get the PSF.

In [71]:

psf = visit_image.getPsf()

At pixel coordinates x, y = 100, 100, compute the PSF and display it in frame 5.

In [72]:

xy = lsst.geom.Point2D(100, 100)
psf_image = psf.computeImage(xy)
afw_display = afwDisplay.Display(frame=5)
afw_display.mtv(psf_image.convertF())

Compute the shape of the PSF. Print the size in $\sigma$ and in FWHM (full-width at half-max), and the second moments (all units are pixels).

In [73]:

psf_shape = psf.computeShape(xy)
psf_sigma = psf_shape.getDeterminantRadius()
psf_fwhm = psf_sigma * 2.0*np.sqrt(2.0*np.log(2.0))
print(f"PSF sigma: {psf_sigma}, PSF FWHM: {psf_fwhm}")
print("PSF shape (ixx, iyy, ixy): ",psf_shape.getIxx(), psf_shape.getIyy(), psf_shape.getIxy())

PSF sigma: 1.7530962025790353, PSF FWHM: 4.1282260787007505
PSF shape (ixx, iyy, ixy):  2.9027495411187836 3.2633505857877685 0.16502109974125848

Clean up.

In [74]:

del psf, xy, psf_image, psf_shape, psf_sigma, psf_fwhm

5. Detected sources¶

In general the TAP service is the recommended access mechanism for catalog data.

However, in the case where all the sources detected in a given image are desired, the Butler offers a quick way to retrieve them.

Get the visit and detector numbers for the visit_image.

In [75]:

visit = visit_image.getInfo().getVisitInfo().id
detector = visit_image.getInfo().getDetector().getId()
print(visit, detector)

2024110800247 4

Get all sources in this visit_image.

In [76]:

source = butler.get('source', visit=visit)
print(len(source))

27446

Warning: The source table contains all sources in the 9 detectors of the visit.

Print the number of sources in the detector for this visit_image.

In [77]:

tx = np.where(source['detector'] == detector)[0]
print(len(tx))

3492

Option to show the source table.

In [78]:

# source[tx]

Redisplay the visit_image in frame 1 with no mask.

In [79]:

afw_display = afwDisplay.Display(frame=1)
afw_display.mtv(visit_image)
afw_display.setMaskTransparency(100)

Overplot sources in the visit_image with orange circles.

In [80]:

with afw_display.Buffering():
    for i in tx:
        afw_display.dot('o', source[i]['x'], source[i]['y'],
                        size=20, ctype='orange')

Clean up.

In [81]:

del visit_image
del visit, detector, source, tx

In [ ]:

In [ ]: