103.2. Simple Image Access (SIA)#
103.2. Simple Image Access (SIA) service¶
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.
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.
target_ra = 53.076
target_dec = -28.110
Set afwDisplay to use Firefly, and define afw_display to show images in frame 1.
afwDisplay.setDefaultBackend("firefly")
afw_display = afwDisplay.Display(frame=1)
2. Instantiate the SIA2 service¶
Instantiate the SIA2 service and assert that it exists.
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.
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.
circle = (target_ra, target_dec, 0.05)
Run the SIA2 query.
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.
results.to_table()
| 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.
del results
Option to pass the coordinates to the SIA2 service as an astropy SkyCoord.
# 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()
# del coord, coord_circle, results
Option to define search boundaries as maximum and minimum coordinates.
# 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()
# del bounds, results
Option to define a polygon with four points.
# 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()
# 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.
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.
# 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.
# results
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.
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.
# results
Print the number of unique values of lsst_band (filters).
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.
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.
results = sia_service.search(pos=circle, calib_level=3, dpsubtype='lsst.template_coadd')
print(len(results))
12
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.
results = sia_service.search(pos=circle, calib_level=3, dpsubtype='lsst.difference_image')
print(len(results))
1656
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".
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.
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.
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.
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
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.
results = sia_service.search(pos=circle, band=(5e-07), calib_level=2)
print(len(results))
414
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.
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.
results = sia_service.search(pos=circle, calib_level=2,
time=(time1, time2))
print(len(results))
131
Option to display the results table.
# results.to_table()
Delete the query constraints but keep the results table to use in Section 4.
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.
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.
# datalink_url
# dl_result
# 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.
visit_image = ExposureF(image_url)
Display the exposure in the Firefly tab.
afw_display.mtv(visit_image)
Set the mask transparency to 100% transparent.
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.