103.2. Simple Image Access (SIA)

103.2. Simple Image Access (SIA) service¶

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: How to use the SIA2 service to access image data.

LSST data products: deep_coadd, visit_image

Packages: lsst.rsp.utils, lsst.rsp.service, pyvo.dal, lsst.afw.image

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¶

Simple Image Access (SIA) is a protocol of the International Virtual Observatory Alliance (IVOA). It provides a standardized model for image metadata, and the capability to query and retrieve image datasets. Learn more in the IVOA SIA documentation.

Why use the SIA service?

SIA is an image access service that is broadly-used and is relatively easy to get started with.

When to use the SIA service.

Use the SIA service to find, retrieve, and display whole images, based on their spatial, temporal, and spectral coverage. SIA is good for quick access to images' pixel and header data, for display and examination by eye.

If only a small part of the image is needed, use the image cutouts tool instead.

For science cases that involve source detection or image analysis with the LSST Science Pipelines, the butler should instead be used to find and retrieve images.

The SIA Version 2 service is used in this tutorial (SIA2).

Related tutorials: Other 100-level tutorials demonstrate how to access images via other services, such as ObsTAP and the Butler. The 100-level tutorials on image display have more on using Firefly. The 200-level tutorials on coadded and visit images have more on these data products (e.g., the pixel data and metadata).

1.1. Import packages¶

Import common scientific analysis packages numpy and astropy.

Import LSST Science Pipelines packages for image display, image type conversion, and utilities for remote data access.

Import pyvo packages for working with the SIA2Service.

In [1]:

import numpy as np
from astropy.coordinates import SkyCoord
from astropy.time import Time
import lsst.afw.display as afwDisplay
from lsst.afw.image import ExposureF
from lsst.rsp.service import get_siav2_service
from lsst.rsp.utils import get_pyvo_auth
from pyvo.dal.adhoc import DatalinkResults

1.2. Define parameters and functions¶

Define search coordinates right ascension (target_ra) and declination (target_dec), in degrees, to use in all queries. These coordinates are near the center of the DP1 ECDFS field.

In [2]:

target_ra = 53.076
target_dec = -28.110

Set afwDisplay to use Firefly, and define afw_display to show images in frame 1.

In [3]:

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

2. Instantiate the SIA2 service¶

Instantiate the SIA2 service and assert that it exists.

In [4]:

sia_service = get_siav2_service("dp1")
assert sia_service is not None

3. Query for images¶

It is recommended to tightly constrain image queries, so that they return only the image data products needed for a given scientific analysis.

The main constraints on SIA2 Service queries are:

• position / coordinate / region overlap
• image type (calibration level: raw, processed, coadded)
• band (filter, spectral coverage)
• time (date and time of acquisition)

Find a list of all possible constraints, and their formats, in the SIA2Query documentation.

3.1. Position¶

Query by position.

To query only by position is not recommended because it would return every image that overlaps that region, including raw, direct, difference, and coadded images. Typically, for scientific analyses, a specific image type or band is desired and should be included in the query, as demonstrated in following sections.

For these examples, set the maximum number of records (maxrec) to 3.

In [5]:

use_maxrec = 3

The position is a region, defined as a 3-, 4-, or $\geq$6-element tuple:

• circle: (ra, dec, radius)
• range: (long1, long2, lat1, lat2)
• polygon: (ra, dec, ra, dec, ra, dec, ...)

Define pos as a circle using a tuple of right ascension, declination, and radius in degrees.

In [6]:

circle = (target_ra, target_dec, 0.05)

Run the SIA2 query.

In [7]:

results = sia_service.search(pos=circle, maxrec=use_maxrec)

Show the results as a table.

Notice: Some exposures of ECDFS erroneously have a target name of slew_icrs, but they are of ECDFS.

In [8]:

results.to_table()

Out[8]:

Table length=3

dataproduct_type	dataproduct_subtype	facility_name	calib_level	target_name	obs_id	obs_collection	obs_publisher_did	access_url	access_format	s_resolution	s_xel1	s_xel2	t_xel	t_min	t_max	t_exptime	t_resolution	em_xel	em_min	em_max	em_res_power	em_filter_name	o_ucd	pol_xel	instrument_name	lsst_visit	lsst_detector	lsst_tract	lsst_patch	lsst_band	lsst_filter	obs_title	s_ra	s_dec	s_fov	s_region
										arcsec				d	d	s	s		m	m													deg	deg	deg
object	object	object	int32	object	object	object	object	object	object	float64	int64	int64	int64	float64	float64	float64	float64	int64	float64	float64	float64	object	object	int64	object	int64	int64	int64	int64	object	object	object	float64	float64	float64	object
image	lsst.raw	Rubin:Simonyi	1	slew_icrs	CC_O_20241108_000266	LSST.DP1	ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=c9b72b48-0087-57a6-bd0c-2120836aedd4	https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3Dc9b72b48-0087-57a6-bd0c-2120836aedd4	application/x-votable+xml;content=datalink	--	4608	4096	--	60623.27366114595	60623.27401335648	30.0	--	--	4.026e-07	5.483e-07	--	g	phot.count	--	LSSTComCam	2024110800266	0	--	--	g	g_01	raw - g - CC_O_20241108_000266-R22_S00 2024-11-09T06:34:04.323010Z	52.99197077057104	-28.11438844236187	0.3179712897409338	POLYGON ICRS 52.898024 -27.978909 52.839871 -28.199792 53.086238 -28.249758 53.143786 -28.028839
image	lsst.raw	Rubin:Simonyi	1	slew_icrs	CC_O_20241108_000270	LSST.DP1	ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=041a0fb4-7358-5457-bcde-e453ef64d554	https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D041a0fb4-7358-5457-bcde-e453ef64d554	application/x-votable+xml;content=datalink	--	4608	4096	--	60623.28061678261	60623.280969085645	30.0	--	--	4.026e-07	5.483e-07	--	g	phot.count	--	LSSTComCam	2024110800270	0	--	--	g	g_01	raw - g - CC_O_20241108_000270-R22_S00 2024-11-09T06:44:05.290018Z	52.909077390196416	-28.180548573758742	0.3179909938756991	POLYGON ICRS 52.815074 -28.045065 52.756875 -28.265958 53.003408 -28.315918 53.060993 -28.094997
image	lsst.raw	Rubin:Simonyi	1	slew_icrs	CC_O_20241108_000274	LSST.DP1	ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=f1c33d01-37b6-5ba3-83e7-78b403bd7985	https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3Df1c33d01-37b6-5ba3-83e7-78b403bd7985	application/x-votable+xml;content=datalink	--	4608	4096	--	60623.28380763904	60623.28415983796	30.0	--	--	4.026e-07	5.483e-07	--	g	phot.count	--	LSSTComCam	2024110800274	0	--	--	g	g_01	raw - g - CC_O_20241108_000274-R22_S00 2024-11-09T06:48:40.980013Z	52.90653061685272	-28.044048182264174	0.3179661367732995	POLYGON ICRS 52.812634 -27.908564 52.754536 -28.129455 53.000743 -28.179420 53.058246 -27.958502

Keep circle to reuse below, but delete the results.

In [9]:

del results

Option to pass the coordinates to the SIA2 service as an astropy SkyCoord.

In [10]:

# coord = SkyCoord(target_ra, target_dec, unit="deg")
# coord_circle = (coord, 0.05)
# results = sia_service.search(pos=coord_circle, maxrec=use_maxrec)
# results.to_table()

In [11]:

# del coord, coord_circle, results

Option to define search boundaries as maximum and minimum coordinates.

In [12]:

# bounds = (target_ra-0.05, target_ra+0.05, target_dec-0.05, target_dec+0.05)
# results = sia_service.search(pos=bounds, maxrec=use_maxrec)
# results.to_table()

In [13]:

# del bounds, results

Option to define a polygon with four points.

In [14]:

# polygon = (target_ra-0.05, target_dec-0.05,
#            target_ra+0.05, target_dec-0.05,
#            target_ra+0.05, target_dec+0.05,
#            target_ra-0.05, target_dec+0.05)
# results = sia_service.search(pos=polygon, maxrec=use_maxrec)
# results.to_table()

In [15]:

# del polygon, results

3.2. Image type¶

Query by image type, which in the SIA2 service is specified by the calibration level (calib_level) and by the data product subtype (dpsubtype).

Calibration levels and data product subtypes names are:

• 1 : raw
• 2 : visit_image
• 3 : deep_coadd, template_coadd, difference_image

3.2.1. Visit images¶

Query for visit_images by setting calib_level=2 OR dpsubtype='lsst.visit_image', since there is only one dataproduct subtype (visit_image) for calibration level 2.

Return the results as Astropy tables using .to_table().

Set calib_level.

In [16]:

results = sia_service.search(pos=circle, calib_level=2).to_table()
print(len(results))

1656

Option to redo the query by setting dpsubtype to confirm they return the same number of visit images.

In [17]:

# results = sia_service.search(pos=circle, dpsubtype='lsst.visit_image').to_table()
# print(len(results))

Option to show the results table, which will be displayed as truncated.

In [18]:

# results

In [19]:

del results

3.2.2. Deep coadd images¶

Query for deep_coadd images by setting calib_level=3 and dpsubtype='lsst.deep_coadd'.

Do not pass use_maxrec here, retrieve all the deep_coadd images that match the search constraints.

In [20]:

results = sia_service.search(pos=circle, calib_level=3, dpsubtype='lsst.deep_coadd').to_table()
print(len(results))

12

Option to show the full table of results.

In [21]:

# results

Print the number of unique values of lsst_band (filters).

In [22]:

values, counts = np.unique(results['lsst_band'],
                           return_counts=True)
for value, count in zip(values, counts):
    print(value, count)

g 2
i 2
r 2
u 2
y 2
z 2

There are 2 overlapping patches for the search coordinates, and thus 2 deep_coadd images for each of the six filters.

In [23]:

del results, values, counts

3.2.3. Template coadd images¶

Query for template_coadd images by setting calib_level=3 and dpsubtype='lsst.template_coadd'.

Since the template_coadd images are also stored by patch, there are also 12 that match the query constraints.

In [24]:

results = sia_service.search(pos=circle, calib_level=3, dpsubtype='lsst.template_coadd')
print(len(results))

12

In [25]:

del results

3.2.4. Difference images¶

Query for difference_images by setting calib_level=3 and dpsubtype='lsst.difference_image'.

There will be as many difference images returned as there are visit_images.

In [26]:

results = sia_service.search(pos=circle, calib_level=3, dpsubtype='lsst.difference_image')
print(len(results))

1656

In [27]:

del results

3.3. Band¶

Query by band.

With the SIA2 Service, band is not the letter name of a filter, but a wavelength range (in meters).

For LSST, the filter transmission curves barely overlap and these values can be used as the band "edges".

In [28]:

band_edges = {'u': (3.000e-07, 3.932e-07),
              'g': (4.026e-07, 5.483e-07),
              'r': (5.510e-07, 6.891e-07),
              'i': (6.936e-07, 8.188e-07),
              'z': (8.192e-07, 9.201e-07),
              'y': (9.278e-07, 1.100e-06)}

Define the band to pass to the SIA2 service. Use the lower end of the LSST r-band, and the upper end of the LSST i-band, to retrieve both r- and i-band images.

In [29]:

band = (band_edges['r'][0], band_edges['i'][1])
print(band)

(5.51e-07, 8.188e-07)

Query the SIA2 service, and convert the results to an astropy table.

In [30]:

results = sia_service.search(pos=circle, band=band, calib_level=2)
table = results.to_table()

Print the number of r- and i-band images returned.

In [31]:

values, counts = np.unique(table['lsst_band'], return_counts=True)
for value, count in zip(values, counts):
    print(value, count)
del values, counts

i 321
r 466

In [32]:

del results, table

Alternatively, define a single wavelength (5e-07 m is in the g-band) and the SIA2 service will return all images with a bandpass that overlaps that wavelength.

In [33]:

results = sia_service.search(pos=circle, band=(5e-07), calib_level=2)
print(len(results))

414

In [34]:

del band, results

3.4. Time¶

The date and time of image acquisition can only be constrained for raw and visit_images (calibration levels 1 and 2).

Define time1 and time2 using astropy Time.

Use International Atomic Time (Temps Atomique International; TAI) to define astropy times, as the Rubin data image acquisition times are stored as TAI times. Times can be defined as calendar dates or MJD, as in time1 and time2 below.

In [35]:

time1 = Time("2024-11-09T00:00:00.0", format="isot", scale="tai")
time2 = Time(60624.0, format="mjd", scale="tai")

Query for visit_images that overlap circle and were obtained between time1 and time2.

The length of the results table should be 131.

In [36]:

results = sia_service.search(pos=circle, calib_level=2,
                             time=(time1, time2))
print(len(results))

131

Option to display the results table.

In [37]:

# results.to_table()

Delete the query constraints but keep the results table to use in Section 4.

In [38]:

del time1, time2

4. Retrieve and display an image¶

Use the results from the query in Section 3.4.

Extract the access URL from the first row of the results as datalink_url.

Retrieve the datalink VOTable document as dl_result. At this point the service requires authorization credentials, which are passed with the utility function get_pyvo_auth() from the lsst.rsp package.

Get the URL for the specific image as image_url.

In [39]:

datalink_url = results[0].access_url
dl_result = DatalinkResults.from_result_url(datalink_url, session=get_pyvo_auth())
image_url = dl_result.getrecord(0).get('access_url')

Option to display the URLs.

In [40]:

# datalink_url

In [41]:

# dl_result

In [42]:

# image_url

Read the entire image into an ExposureF object. The ExposureF format is designed to work with seamlessly with afw.display and so is recommended for all Rubin images.

In [43]:

visit_image = ExposureF(image_url)

Display the exposure in the Firefly tab.

In [44]:

afw_display.mtv(visit_image)

Set the mask transparency to 100% transparent.

In [45]:

afw_display.setMaskTransparency(100)

5. Exercises for the learner¶

The approximate central coordinates of the Rubin_SV_95_-25 field for Data Preview 1 are RA, Dec = $95, -25$. Show that the number of visit images obtained on MJD 60634 that overlap these coordinates was 43.

In [ ]: