308.5. Image stamps#
308.5. Image stamps¶
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: 2026-06-16
Repository: github.com/lsst/tutorial-notebooks
Learning objective: Generate many (bulk) image stamps for Solar System objects using the Rubin image cutout service.
LSST data products: DiaObject, SSObject, MPCORB, Visit, visit_image, template_coadd, difference_image
Packages: matplotlib, numpy, lsst.rsp,lsst.afw, lsst.geom, lsst.resources, lsst.utils.plotting, pyvo.dal.adhoc, reproject
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¶
This notebook demonstrates how to produce bulk image stamps for Solar System objects using the Rubin cutout service. This notebook differs from the 103-series notebooks on image access and display by demonstrating how to generate cutouts for objects that move to different ra and dec positions across each detection.
Related tutorials: See also tutorial notebook 103.9 on image stamps, which demonstrates how to perform bulk cutouts using the cutout service that returns only the image extension to minimize the data transferred and tutorial notebook 103.4 that demonstrates the return of exposure cutouts as Rubin LSST image type ExposureF and MaskedImageF. Both of these notebooks retrieve image stamps for a single ra and dec.
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).
From the lsst package, import modules for accessing the Table Access Protocol (TAP) service, image display functions, and the Simple Image Access (SIA) service from the LSST Science Pipelines (pipelines.lsst.io).
Functions from pyvo are imported to perform and parse results from Datalink and SODA queries. This is used to obtain images from the SIA service.
The reproject package is an astropy-adjacent package that provides access to coordinate transformations between images with world coordinate systems (WCS), and will be used to align visit images with different WCS.
import math
import pandas as pd
import numpy as np
import io
import getpass
import glob
import os
import shutil
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from lsst.rsp import get_tap_service
from lsst.rsp.service import get_siav2_service
from lsst.rsp.utils import get_pyvo_auth
import lsst.afw.display as afwDisplay
from lsst.afw.image import ExposureF, ImageF
from lsst.afw.fits import MemFileManager
from pyvo.dal.adhoc import DatalinkResults, SodaQuery
from astropy import units as u
from astropy.wcs import WCS
from astropy.io import fits
from reproject import reproject_interp
from IPython.display import Image as dimg
from PIL import Image
1.2. Define parameters and functions¶
Create an instance of the TAP service, and assert that it exists.
service = get_tap_service("tap")
assert service is not None
Instantiate the Simple Image Access service (SIA).
sia_service = get_siav2_service("dp1")
assert sia_service is not None
Set the afwDisplay package backend to matplotlib.
afwDisplay.setDefaultBackend('matplotlib')
1.2.1. Set up a temporary directory¶
The Rubin cutout service allows the user to save many image stamps as fits files locally on the Rubin Science Platform (RSP).
RSP users should save temporary files in the shared scratch directory, /deleted-sundays/, in sub-directories named with their RSP username.
Get the username and ensure an appropriate /deleted-sundays/ sub-directory exists.
username = getpass.getuser()
print('username: ', username)
userdir = '/deleted-sundays/' + username
if not os.path.exists(userdir):
os.makedirs(userdir)
print('Created ', userdir)
else:
print('Directory already existed: ', userdir)
username: sarahgreenstreet Directory already existed: /deleted-sundays/sarahgreenstreet
Within that folder, create a sub-folder named dp1_308_5_temp in which to save files created by this tutorial.
At the end of the notebook, the last cell will clear the contents and remove this temporary folder.
tempdir = userdir + '/dp1_308_5_temp'
if not os.path.exists(tempdir):
os.makedirs(tempdir)
print('Created ', tempdir)
else:
print('Directory already existed: ', tempdir)
Directory already existed: /deleted-sundays/sarahgreenstreet/dp1_308_5_temp
Delete the username and userdir, but keep tempdir to be used.
del username, userdir
1.2.2. Generate a gif from pngs¶
The function below will generate a gif to visualize all of the image stamp files in png format produced in the notebook.
def make_gif(folder):
"""
Generate a GIF for all *.png images contained in a folder
Parameters
----------
folder : 'string'
string containing the path to the folder
default filename is animation.gif
Returns
-------
"""
frames = [Image.open(img) for img in sorted(glob.glob(f"{folder}/*.png"))]
frame_one = frames[0]
frame_one.save(folder+"/animation.gif", format="GIF",
append_images=frames, save_all=True, duration=500, loop=1)
1.2.3. Generate cutouts¶
Define a function to get an image cutout from the cutout service.
def get_cutout(dl_result, ra, dec, session, fov):
"""
Wrapper function to generate a cutout using the cutout tool.
Parameters
----------
dl_result : 'pyvo.dal.DatalinkResults' The Datalink result object
containing the 'cutout-sync' ad-hoc service descriptor.
ra, dec : 'float'
the ra and dec of the cutout center
session : 'requests.Session'
An authenticated session object required to
access proprietary or restricted Rubin data services.
fov : 'float'
the edge size of the cutout in decimal degrees
Returns
-------
ImageF(mem) : 'lsst.afw.image.ImageF'
The cutout image data instantiated directly in memory via MemFileManager.
hdul : 'astropy.io.fits.HDUList'
The identical cutout data parsed into
an Astropy FITS object for standard header and extension manipulation.
"""
sq = SodaQuery.from_resource(dl_result,
dl_result.get_adhocservice_by_id("cutout-sync"),
session=session)
sq.circle = (ra * u.deg, dec * u.deg, fov / 2. * u.deg)
cutout_bytes = sq.execute_stream().read()
sq.raise_if_error()
mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
hdul = fits.open(io.BytesIO(cutout_bytes))
return ImageF(mem), hdul
1.2.4. Set up subplots¶
Define a function to organize the display with matplotlib of bulk image stamps.
def make_subplot_grid(n_subplots, figsize_per_plot=(4, 3)):
"""
Create an optimal grid of matplotlib subplots.
Parameters
----------
n_subplots : int
The total number of subplots required.
figsize_per_plot : tuple of float or int, optional
The (width, height) in inches for each individual subplot.
The total figure size is scaled automatically based on this
and the calculated grid size. Default is (4, 3).
Returns
-------
fig : matplotlib.figure.Figure
The generated matplotlib figure object.
axes : numpy.ndarray of matplotlib.axes.Axes
A flattened, 1D array of exactly `n_subplots` axes objects,
ready for plotting.
"""
n_cols = math.ceil(math.sqrt(n_subplots))
n_rows = math.ceil(n_subplots / n_cols)
figsize = (figsize_per_plot[0] * n_cols,
figsize_per_plot[1] * n_rows)
fig, axes = plt.subplots(
n_rows,
n_cols,
figsize=figsize)
axes = axes.flatten()
for ax in axes[n_subplots:]:
ax.remove()
return fig, axes[:n_subplots]
2. Define target¶
2.1. Find target Solar System object¶
Query the MPCORB table for the object designation, orbital parameters, and ssObjectId for all DP1 Solar System objects.
query = "SELECT mpc.mpcDesignation, "\
"mpc.ssObjectId, "\
"mpc.q, mpc.e, mpc.incl, mpc.mpcH, "\
"sso.discoverySubmissionDate, "\
"sso.numObs "\
"FROM dp1.MPCORB as mpc "\
"INNER JOIN dp1.SSObject as sso "\
"ON mpc.ssObjectId = sso.ssObjectId "\
"ORDER BY mpc.mpcDesignation"
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
Fetch the job results and assign to an astropy result table. There are 431 Solar System objects in DP1.
assert job.phase == 'COMPLETED'
result = job.fetch_result()
print(len(result))
431
Convert the result astropy table to a pandas dataframe.
result_df = pd.DataFrame(result)
Print the results, sorting the table based on the parameters of interest. For example, find a new DP1 discovery with ~30 observations in total.
df_discoveries = result_df.loc[result_df['discoverySubmissionDate'] > 0]
df_discoveries.sort_values(by=["numObs"]).tail(10)
| discoverySubmissionDate | e | incl | mpcDesignation | mpcH | numObs | q | ssObjectId | |
|---|---|---|---|---|---|---|---|---|
| 350 | 60853.579024 | 0.106758 | 6.845302 | 2024 WE104 | 17.799000 | 29 | 2.738098 | 23133931615302197 |
| 353 | 60853.579024 | 0.111326 | 5.762680 | 2024 WF103 | 19.474001 | 29 | 1.971333 | 23133931615301683 |
| 389 | 60853.579024 | 0.044314 | 7.499122 | 2024 WP103 | 17.933001 | 29 | 3.027090 | 23133931615301936 |
| 357 | 60853.579024 | 0.151044 | 5.355632 | 2024 WG103 | 18.221001 | 29 | 2.347606 | 23133931615301684 |
| 413 | 60853.579024 | 0.296358 | 6.391025 | 2024 WV103 | 20.608000 | 30 | 1.796838 | 23133931615301942 |
| 417 | 60853.579024 | 0.189241 | 12.220148 | 2024 WW103 | 18.768000 | 30 | 2.189507 | 23133931615301943 |
| 373 | 60853.579024 | 0.228471 | 4.587451 | 2024 WL103 | 20.250000 | 32 | 2.024287 | 23133931615301687 |
| 349 | 60853.579024 | 0.070071 | 24.063976 | 2024 WE103 | 20.363001 | 63 | 1.782912 | 23133931615301682 |
| 365 | 60853.579024 | 0.058592 | 24.350102 | 2024 WJ103 | 19.969999 | 65 | 1.813405 | 23133931615301681 |
| 361 | 60853.579024 | 0.109862 | 22.242227 | 2024 WH103 | 20.798000 | 82 | 1.602490 | 23133931615301680 |
An object that satisfies our parameters of interest is 2024 WP103, which is a Rubin DP1 discovery with 29 observations.
2.2. Find Solar System object visits¶
Query the DiaSource catalog for the ra, dec, band, detector, and visit for each observation of 2024 WP103, inputting the ssObjectId (23133931615301936) into the query below. Order the query results by the observation midpoint.
query = "SELECT mpc.mpcDesignation, "\
"mpc.q, mpc.e, mpc.incl, mpc.mpcH, "\
"mpc.ssObjectId, "\
"diasource.ssObjectId, "\
"diasource.ra, diasource.dec, "\
"diasource.band, diasource.psfFlux, "\
"diasource.detector, diasource.visit, "\
"diasource.midpointMjdTai "\
"FROM dp1.DiaSource as diasource "\
"INNER JOIN dp1.MPCORB as mpc "\
"ON diasource.ssObjectId = mpc.ssObjectId "\
"WHERE diasource.ssObjectId = 23133931615301936 "\
"ORDER BY diasource.midpointMjdTai"
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
Fetch the job results and assign to an astropy result table. There should be 29 observations for Rubin discovery 2024 WP103 (ssObjectId = 23133931615301936) in DP1.
assert job.phase == 'COMPLETED'
result = job.fetch_result()
print(len(result))
29
Convert the result to an astropy table.
tab = result.to_table()
Option to print results.
tab
| mpcDesignation | q | e | incl | mpcH | ssObjectId | ssObjectId2 | ra | dec | band | psfFlux | detector | visit | midpointMjdTai |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AU | deg | mag | deg | deg | nJy | d | |||||||
| str13 | float64 | float64 | float64 | float32 | int64 | int64 | float64 | float64 | str1 | float32 | int64 | int64 | float64 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.258240927869565 | 6.7311106548718795 | r | 3552.95 | 1 | 2024112300221 | 60638.125741319556 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.25778791619301 | 6.731071010520855 | r | 4183.06 | 2 | 2024112300225 | 60638.128363188574 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.256403779969375 | 6.7309971979448555 | i | 5919.87 | 4 | 2024112300232 | 60638.13631882514 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.25597632564128 | 6.731023054495113 | i | 5639.46 | 1 | 2024112300236 | 60638.1387492129 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.25555901706945 | 6.731041699580604 | i | 5753.86 | 2 | 2024112300240 | 60638.14114406817 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.25512832059666 | 6.730973536295063 | i | 5987.14 | 4 | 2024112300244 | 60638.14361382516 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 38.254710852792726 | 6.730955279313136 | i | 6088.61 | 4 | 2024112300248 | 60638.145979554436 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.78750030181114 | 6.720030654357857 | r | 3692.69 | 3 | 2024112600104 | 60641.0479109837 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.78705042613154 | 6.720030968896434 | r | 3488.95 | 3 | 2024112600108 | 60641.05036715285 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.63205639071729 | 6.720145850719419 | i | 5875.05 | 5 | 2024112700100 | 60642.06069765047 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47941042896497 | 6.722362632977377 | r | 3757.06 | 4 | 2024112800165 | 60643.08383392948 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.479062569886985 | 6.722376736275439 | r | 4002.98 | 4 | 2024112800169 | 60643.086150856456 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47864682376877 | 6.722364010508033 | g | 2184.4 | 4 | 2024112800173 | 60643.088763668995 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47829356337315 | 6.722365547656145 | g | 2269.94 | 4 | 2024112800177 | 60643.09108065391 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.4779665648256 | 6.722378524568007 | g | 2299.75 | 4 | 2024112800181 | 60643.093257951274 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47753195167145 | 6.722411922683821 | i | 4742.57 | 4 | 2024112800185 | 60643.09600285883 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.477217735767994 | 6.722407035931727 | i | 5355.75 | 4 | 2024112800189 | 60643.098207766234 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.476864816801466 | 6.722404539178317 | i | 5354.24 | 4 | 2024112800193 | 60643.100482164285 |
Restrict the detections to those in a single band. Make a new table with just the g-band observations (9 in total). Store the provisional designation (mpcDesignation= 2024 WP103) for labeling image stamp displays below.
tab_g = tab[tab['band'] == 'g']
mpcDesignation = tab_g['mpcDesignation'][0]
Option to print results.
tab_g
| mpcDesignation | q | e | incl | mpcH | ssObjectId | ssObjectId2 | ra | dec | band | psfFlux | detector | visit | midpointMjdTai |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AU | deg | mag | deg | deg | nJy | d | |||||||
| str13 | float64 | float64 | float64 | float32 | int64 | int64 | float64 | float64 | str1 | float32 | int64 | int64 | float64 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.78657276684927 | 6.720035379049524 | g | 2243.77 | 3 | 2024112600112 | 60641.05333741903 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.785457578599555 | 6.720055682159962 | g | 1933.17 | 3 | 2024112600115 | 60641.06044383095 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.785152892884824 | 6.720062593376525 | g | 1868.57 | 0 | 2024112600118 | 60641.06222083327 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.63352968623554 | 6.720153708835776 | g | 2543.87 | 5 | 2024112700084 | 60642.0513119733 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.63316030920984 | 6.720143839512422 | g | 2559.71 | 5 | 2024112700088 | 60642.0535819272 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.63286554978612 | 6.720175714240858 | g | 2447.1 | 5 | 2024112700092 | 60642.055765798636 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47864682376877 | 6.722364010508033 | g | 2184.4 | 4 | 2024112800173 | 60643.088763668995 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.47829356337315 | 6.722365547656145 | g | 2269.94 | 4 | 2024112800177 | 60643.09108065391 |
| 2024 WP103 | 3.02709038379305 | 0.044314099044914 | 7.4991222397328 | 17.933 | 23133931615301936 | 23133931615301936 | 37.4779665648256 | 6.722378524568007 | g | 2299.75 | 4 | 2024112800181 | 60643.093257951274 |
3. Generate cutouts¶
3.1. Query ObsCore table for first visit¶
To generate a cutout for the first g-band detection for the Solar System target, query the ObsCore table for the first visit (science) image for the target.
For stationary targets, this can be done using the ra and dec coordinates of the target. Solar System targets, however, move to different ra and dec positions between detections, and many images likely cover those (ra, dec) coordinates at times when the moving object was not at that location. To be more specific in the image query to only return images when a Solar System object was at a given location, use the visit and detector parameters, because they correspond to the images taken at a given sky location when the moving Solar System object was at that position.
To generate a cutout, first query the International Virtual Observatory Alliance (IVOA) ObsCore table to retrieve the access url that points to the location of the image on the remote server. Use the visit and detector of the first detection in the tab_g table above.
query = "SELECT dataproduct_type, dataproduct_subtype, calib_level, lsst_band, em_min, em_max, lsst_tract, lsst_patch, "\
"lsst_filter, lsst_visit, lsst_detector, t_exptime, t_min, t_max, s_ra, s_dec, s_fov, obs_id, "\
"obs_collection, o_ucd, facility_name, instrument_name, obs_title, s_region, access_url, access_format "\
"FROM ivoa.ObsCore WHERE lsst_visit = "+str(tab_g['visit'][0])+" AND lsst_detector = "+str(tab_g['detector'][0])+" AND "\
"obs_collection = 'LSST.DP1' AND calib_level = 2 AND dataproduct_type = 'image' AND "\
"instrument_name = 'LSSTComCam' AND "\
"dataproduct_subtype = 'lsst.visit_image' "\
"ORDER by t_min 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()
first_visit = job.fetch_result()
first_visit.to_table()
Job phase is COMPLETED
| dataproduct_type | dataproduct_subtype | calib_level | lsst_band | em_min | em_max | lsst_tract | lsst_patch | lsst_filter | lsst_visit | lsst_detector | t_exptime | t_min | t_max | s_ra | s_dec | s_fov | obs_id | obs_collection | o_ucd | facility_name | instrument_name | obs_title | s_region | access_url | access_format |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| m | m | s | d | d | deg | deg | deg | ||||||||||||||||||
| object | object | int32 | object | float64 | float64 | int64 | int64 | object | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 | object | object | object | object | object | object | object | object | object |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600112 | 3 | 30.0 | 60641.052733032506 | 60641.05308532407 | 37.795749255533444 | 6.654049228050561 | 0.3176848500023962 | CC_O_20241126_000112 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241126_000112-R22_S10 2024-11-27T01:15:56.134008Z | POLYGON ICRS 37.951751 6.688922 37.757816 6.808360 37.639760 6.619084 37.833691 6.499764 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1977-c18e-7eec-8837-ee2e142137f1 | application/x-votable+xml;content=datalink |
Retrieve the access url from the query (or datalink) which points to the image location on the remote server. Authorization to access the image is required and provided using get_pyvo_auth().
dl_result = DatalinkResults.from_result_url(first_visit['access_url'][0],
session=get_pyvo_auth())
Define the image cutout to be 0.004 degrees (~14 arcseconds) on a side.
fov = 0.004
Use the IVOA Server-side Operations for Data Access (SODA) protocol to make and return the cutout centered on the ra and dec coordinates of the first visit in the tab_g table above.
sq = SodaQuery.from_resource(dl_result,
dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
session=get_pyvo_auth())
sq.circle = (tab_g['ra'][0] * u.deg, tab_g['dec'][0] * u.deg, fov * u.deg)
cutout_bytes = sq.execute_stream().read()
sq.raise_if_error()
Receive the cutout "bytes" into memory, and convert it to an ExposureF image type.
mem = MemFileManager(len(cutout_bytes))
mem.setData(cutout_bytes, len(cutout_bytes))
exposure = ExposureF(mem)
Display the science image extension of the ExposureF cutout using the Rubin/LSST image display package, afwDisplay.
display = afwDisplay.Display()
display.scale('linear', 'zscale')
display.image(exposure.image)
plt.show()
Figure 1: Science image extension of the
ExposureFcutout for the first g-band observation of Rubin discovery 2024 WP103.
3.2. Generate cutout in different format¶
Now, return the cutout in ImageF format instead of ExposureF. This is the same procedure, except one calls get_adhocservice_by_id("cutout-sync") instead of get_adhocservice_by_id("cutout-sync-exposure"), and the return format will be like a fits file. Use astropy and matplotlib tools to plot instead of LSST pipeline tools.
sci, scihdul = get_cutout(dl_result, tab_g['ra'][0], tab_g['dec'][0], get_pyvo_auth(), fov)
sci_header = scihdul[1].header
scidata = scihdul[1].data
plt.subplot(1, 2, 1, projection=WCS(scihdul[1].header))
plt.imshow(scidata, origin='lower', vmin=np.nanpercentile(scidata, 1),
vmax=np.nanpercentile(scidata, 99))
<matplotlib.image.AxesImage at 0x7c53e8392db0>
Figure 2:
ImageFcutout format for the first g-band observtion of Rubin discovery 2024 WP103.
4. Create single-epoch triplet stamps¶
4.1. Query ObsCore table for all visits¶
To create a triplet stamp (i.e., the science, template, and difference images) for the first g-band detection of the Solar System target, query the ObsCore table for the visit (science), template_coadd, and difference images for the target. Rather than restricting the queries to the first detections, return all visit, template, and difference images for the target. This will avoid the need to run the queries again in Section 5 below when bulk image stamps are generated.
As described in Section 3.1 above, use the visit and detector parameters in the image query to return images containing the Solar System target.
Set the visit and detector conditions that cutout queries must satisfy to match those for each g-band detection of the Solar System object in the tab_g table. Returned visit images must match both the visit and detector for each detection.
pointings = []
for row in tab_g:
pointings.append((row['visit'], row['detector']))
pointing_conditions = []
for visit, detector in pointings:
cond = f"(lsst_visit = {visit} AND lsst_detector = {detector})"
pointing_conditions.append(cond)
pointings_sql = " OR ".join(pointing_conditions)
Build the query for the full list of visit images for the Solar System object, matching both the visit and detector for each detection. Sort the returned list of visit images by time.
query = f"""
SELECT dataproduct_type, dataproduct_subtype, calib_level, lsst_band, em_min, em_max,
lsst_tract, lsst_patch, lsst_filter, lsst_visit, lsst_detector, t_exptime, t_min, t_max,
s_ra, s_dec, s_fov, obs_id, obs_collection, o_ucd, facility_name, instrument_name,
obs_title, s_region, access_url, access_format
FROM ivoa.ObsCore
WHERE obs_collection = 'LSST.DP1'
AND calib_level = 2
AND dataproduct_type = 'image'
AND instrument_name = 'LSSTComCam'
AND dataproduct_subtype = 'lsst.visit_image'
AND ({pointings_sql})
ORDER BY t_min ASC
"""
Execute the query. Store the results in a table called scitab.
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
if job.phase == 'ERROR':
job.raise_if_error()
scitab = job.fetch_result()
scitab.to_table()
| dataproduct_type | dataproduct_subtype | calib_level | lsst_band | em_min | em_max | lsst_tract | lsst_patch | lsst_filter | lsst_visit | lsst_detector | t_exptime | t_min | t_max | s_ra | s_dec | s_fov | obs_id | obs_collection | o_ucd | facility_name | instrument_name | obs_title | s_region | access_url | access_format |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| m | m | s | d | d | deg | deg | deg | ||||||||||||||||||
| object | object | int32 | object | float64 | float64 | int64 | int64 | object | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 | object | object | object | object | object | object | object | object | object |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600112 | 3 | 30.0 | 60641.052733032506 | 60641.05308532407 | 37.795749255533444 | 6.654049228050561 | 0.3176848500023962 | CC_O_20241126_000112 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241126_000112-R22_S10 2024-11-27T01:15:56.134008Z | POLYGON ICRS 37.951751 6.688922 37.757816 6.808360 37.639760 6.619084 37.833691 6.499764 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1977-c18e-7eec-8837-ee2e142137f1 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600115 | 3 | 30.0 | 60641.059839455884 | 60641.06019172454 | 37.83524180001764 | 6.673969161299994 | 0.3176664242054425 | CC_O_20241126_000115 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241126_000115-R22_S10 2024-11-27T01:26:10.128988Z | POLYGON ICRS 37.991246 6.708828 37.797303 6.828267 37.679250 6.639019 37.873191 6.519696 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1977-ca2a-7599-832a-c39703e525ad | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600118 | 0 | 30.0 | 60641.06161649293 | 60641.06196869213 | 37.931582332958406 | 6.758994226880448 | 0.3178948841948642 | CC_O_20241126_000118 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241126_000118-R22_S00 2024-11-27T01:28:43.664989Z | POLYGON ICRS 38.087651 6.793757 37.893718 6.913432 37.775560 6.724115 37.969400 6.604585 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1945-bb8d-76b6-8d34-b34edda6b334 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700084 | 5 | 30.0 | 60642.05070755771 | 60642.051059907404 | 37.60090178481414 | 6.769413090260296 | 0.3177300950256694 | CC_O_20241127_000084 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241127_000084-R22_S12 2024-11-28T01:13:01.132986Z | POLYGON ICRS 37.755665 6.809679 37.557543 6.922309 37.446171 6.729138 37.644205 6.616487 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-2821-7daf-b9c0-09e3c882ec63 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700088 | 5 | 30.0 | 60642.05297753495 | 60642.05332983796 | 37.640969646357455 | 6.789995360031279 | 0.317715650572084 | CC_O_20241127_000088 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241127_000088-R22_S12 2024-11-28T01:16:17.259019Z | POLYGON ICRS 37.795733 6.830255 37.597607 6.942892 37.486239 6.749721 37.684278 6.637072 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-1eec-7543-8eea-b5321f4b97b7 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700092 | 5 | 30.0 | 60642.055161446806 | 60642.05551366898 | 37.61705785136126 | 6.674054573851342 | 0.31771511237595107 | CC_O_20241127_000092 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241127_000092-R22_S12 2024-11-28T01:19:25.949004Z | POLYGON ICRS 37.771785 6.714312 37.573703 6.826949 37.462366 6.633779 37.660357 6.521135 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-501f-7775-9518-2e47da7c58c6 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800173 | 4 | 30.0 | 60643.08815931716 | 60643.088511539354 | 37.42815925163586 | 6.769329377688846 | 0.3174988483113268 | CC_O_20241128_000173 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241128_000173-R22_S11 2024-11-29T02:06:56.965002Z | POLYGON ICRS 37.566873 6.848262 37.346184 6.905623 37.289493 6.690356 37.510086 6.633022 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-e986-7e19-aa65-552e1ed22ad2 | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800177 | 4 | 30.0 | 60643.0904762731 | 60643.09082855324 | 37.461253560207915 | 6.799537713185662 | 0.31749993202137033 | CC_O_20241128_000177 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241128_000177-R22_S11 2024-11-29T02:10:17.149996Z | POLYGON ICRS 37.599973 6.878472 37.379276 6.935835 37.322578 6.720562 37.543188 6.663230 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-df7b-7d70-b6eb-d4988e35681c | application/x-votable+xml;content=datalink |
| image | lsst.visit_image | 2 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800181 | 4 | 30.0 | 60643.09265355301 | 60643.09300586805 | 37.46843549716876 | 6.68145654847747 | 0.31751231026779375 | CC_O_20241128_000181 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | visit_image - g - CC_O_20241128_000181-R22_S11 2024-11-29T02:13:25.266980Z | POLYGON ICRS 37.607125 6.760401 37.386474 6.817755 37.329792 6.602476 37.550347 6.545141 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-ef6e-7020-874d-d441f0d73cba | application/x-votable+xml;content=datalink |
Use the SIA service to search for template_coadd (calibration level 3) images that overlap a search circle centered on the ra and dec coordinates for each detection of the Solar System object (in tab_g). Filter the results for g-band template_coadd images to match the band of the visit images. Store the results in a table called reftab.
Assert that the number of returned template_coadd images equals the number of visit images. If the number of returned template images does not match the number of visits, the visit images cannot be rotated to align with the template images in Section 4.3 below.
circle = (tab_g['ra'][0], tab_g['dec'][0], 0.001)
results = sia_service.search(pos=circle, calib_level=3)
lvl3_table = results.to_table()
num_visits = len(scitab['t_max'])
print('number of visits:', num_visits)
for i in range(num_visits-1):
circle = (tab_g['ra'][i+1], tab_g['dec'][i+1], 0.001)
results = sia_service.search(pos=circle, calib_level=3)
tmp_table = results.to_table()
lvl3_table = np.append(lvl3_table, tmp_table, axis=0)
tx = np.where((lvl3_table['dataproduct_subtype'] == 'lsst.template_coadd')
& (lvl3_table['lsst_band'] == 'g'))[0]
reftab = pd.DataFrame(lvl3_table[tx])
assert len(reftab) == num_visits
reftab
number of visits: 9
| dataproduct_type | dataproduct_subtype | facility_name | calib_level | target_name | obs_id | obs_collection | obs_publisher_did | access_url | access_format | ... | lsst_detector | lsst_tract | lsst_patch | lsst_band | lsst_filter | obs_title | s_ra | s_dec | s_fov | s_region | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-52 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 52 | g | template_coadd - g - tract=10463 patch=52 | 37.763960 | 6.777396 | 0.267120 | POLYGON ICRS 37.668794 6.683025 37.858964 6.68... | ||
| 1 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-52 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 52 | g | template_coadd - g - tract=10463 patch=52 | 37.763960 | 6.777396 | 0.267120 | POLYGON ICRS 37.668794 6.683025 37.858964 6.68... | ||
| 2 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-52 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 52 | g | template_coadd - g - tract=10463 patch=52 | 37.763960 | 6.777396 | 0.267120 | POLYGON ICRS 37.668794 6.683025 37.858964 6.68... | ||
| 3 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-53 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 53 | g | template_coadd - g - tract=10463 patch=53 | 37.596127 | 6.777511 | 0.267126 | POLYGON ICRS 37.500990 6.683107 37.691167 6.68... | ||
| 4 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-53 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 53 | g | template_coadd - g - tract=10463 patch=53 | 37.596127 | 6.777511 | 0.267126 | POLYGON ICRS 37.500990 6.683107 37.691167 6.68... | ||
| 5 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-53 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 53 | g | template_coadd - g - tract=10463 patch=53 | 37.596127 | 6.777511 | 0.267126 | POLYGON ICRS 37.500990 6.683107 37.691167 6.68... | ||
| 6 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-54 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 54 | g | template_coadd - g - tract=10463 patch=54 | 37.428290 | 6.777568 | 0.267128 | POLYGON ICRS 37.333183 6.683131 37.523364 6.68... | ||
| 7 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-54 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 54 | g | template_coadd - g - tract=10463 patch=54 | 37.428290 | 6.777568 | 0.267128 | POLYGON ICRS 37.333183 6.683131 37.523364 6.68... | ||
| 8 | image | lsst.template_coadd | Rubin:Simonyi | 3 | lsst_cells_v1-10463-54 | LSST.DP1 | ivo://org.rubinobs/usdac/lsst-dp1?repo=dp1&id=... | https://data.lsst.cloud/api/datalink/links?ID=... | application/x-votable+xml;content=datalink | ... | 0 | 10463 | 54 | g | template_coadd - g - tract=10463 patch=54 | 37.428290 | 6.777568 | 0.267128 | POLYGON ICRS 37.333183 6.683131 37.523364 6.68... |
9 rows × 37 columns
Similar to the visit_images above, build the query to search for (calibration level 3) difference_images that overlap the full list of visit images for the Solar System object, matching the visit and detector conditions for each detection as defined in pointings_sql above. Sort by time to match the visit images (scitab) above.
query = f"""
SELECT dataproduct_type, dataproduct_subtype, calib_level, lsst_band, em_min, em_max,
lsst_tract, lsst_patch, lsst_filter, lsst_visit, lsst_detector, t_exptime, t_min, t_max,
s_ra, s_dec, s_fov, obs_id, obs_collection, o_ucd, facility_name, instrument_name,
obs_title, s_region, access_url, access_format
FROM ivoa.ObsCore
WHERE obs_collection = 'LSST.DP1'
AND calib_level = 3
AND dataproduct_type = 'image'
AND instrument_name = 'LSSTComCam'
AND dataproduct_subtype = 'lsst.difference_image'
AND ({pointings_sql})
ORDER BY t_min ASC
"""
Execute the query. Store the results in a table called difftab.
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
if job.phase == 'ERROR':
job.raise_if_error()
difftab = job.fetch_result()
difftab.to_table()
| dataproduct_type | dataproduct_subtype | calib_level | lsst_band | em_min | em_max | lsst_tract | lsst_patch | lsst_filter | lsst_visit | lsst_detector | t_exptime | t_min | t_max | s_ra | s_dec | s_fov | obs_id | obs_collection | o_ucd | facility_name | instrument_name | obs_title | s_region | access_url | access_format |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| m | m | s | d | d | deg | deg | deg | ||||||||||||||||||
| object | object | int32 | object | float64 | float64 | int64 | int64 | object | int64 | int64 | float64 | float64 | float64 | float64 | float64 | float64 | object | object | object | object | object | object | object | object | object |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600112 | 3 | 30.0 | 60641.052733032506 | 60641.05308532407 | 37.795749255533444 | 6.654049228050561 | 0.3176848500023962 | CC_O_20241126_000112 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241126_000112-R22_S10 2024-11-27T01:15:56.134008Z | POLYGON ICRS 37.951751 6.688922 37.757816 6.808360 37.639760 6.619084 37.833691 6.499764 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1977-cb94-76a4-ac6c-430f0df4b35e | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600115 | 3 | 30.0 | 60641.059839455884 | 60641.06019172454 | 37.83524180001764 | 6.673969161299994 | 0.3176664242054425 | CC_O_20241126_000115 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241126_000115-R22_S10 2024-11-27T01:26:10.128988Z | POLYGON ICRS 37.991246 6.708828 37.797303 6.828267 37.679250 6.639019 37.873191 6.519696 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1977-d42f-77f6-abf6-bf720ce0a5fe | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112600118 | 0 | 30.0 | 60641.06161649293 | 60641.06196869213 | 37.931582332958406 | 6.758994226880448 | 0.3178948841948642 | CC_O_20241126_000118 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241126_000118-R22_S00 2024-11-27T01:28:43.664989Z | POLYGON ICRS 38.087651 6.793757 37.893718 6.913432 37.775560 6.724115 37.969400 6.604585 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1945-c5cd-73b4-8551-a65e768c111f | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700084 | 5 | 30.0 | 60642.05070755771 | 60642.051059907404 | 37.60090178481414 | 6.769413090260296 | 0.3177300950256694 | CC_O_20241127_000084 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241127_000084-R22_S12 2024-11-28T01:13:01.132986Z | POLYGON ICRS 37.755665 6.809679 37.557543 6.922309 37.446171 6.729138 37.644205 6.616487 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-325c-762c-81f0-49def7cf9459 | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700088 | 5 | 30.0 | 60642.05297753495 | 60642.05332983796 | 37.640969646357455 | 6.789995360031279 | 0.317715650572084 | CC_O_20241127_000088 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241127_000088-R22_S12 2024-11-28T01:16:17.259019Z | POLYGON ICRS 37.795733 6.830255 37.597607 6.942892 37.486239 6.749721 37.684278 6.637072 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-292c-73b5-9c93-11cc7c2e0b09 | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112700092 | 5 | 30.0 | 60642.055161446806 | 60642.05551366898 | 37.61705785136126 | 6.674054573851342 | 0.31771511237595107 | CC_O_20241127_000092 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241127_000092-R22_S12 2024-11-28T01:19:25.949004Z | POLYGON ICRS 37.771785 6.714312 37.573703 6.826949 37.462366 6.633779 37.660357 6.521135 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1999-5a34-7475-8443-e97a968c6c41 | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800173 | 4 | 30.0 | 60643.08815931716 | 60643.088511539354 | 37.42815925163586 | 6.769329377688846 | 0.3174988483113268 | CC_O_20241128_000173 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241128_000173-R22_S11 2024-11-29T02:06:56.965002Z | POLYGON ICRS 37.566873 6.848262 37.346184 6.905623 37.289493 6.690356 37.510086 6.633022 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-f3a5-7fac-b994-5b9c5d3947b7 | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800177 | 4 | 30.0 | 60643.0904762731 | 60643.09082855324 | 37.461253560207915 | 6.799537713185662 | 0.31749993202137033 | CC_O_20241128_000177 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241128_000177-R22_S11 2024-11-29T02:10:17.149996Z | POLYGON ICRS 37.599973 6.878472 37.379276 6.935835 37.322578 6.720562 37.543188 6.663230 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-e9a0-70bd-8b5a-857b2d644c9e | application/x-votable+xml;content=datalink |
| image | lsst.difference_image | 3 | g | 4.026e-07 | 5.483e-07 | -- | -- | g_01 | 2024112800181 | 4 | 30.0 | 60643.09265355301 | 60643.09300586805 | 37.46843549716876 | 6.68145654847747 | 0.31751231026779375 | CC_O_20241128_000181 | LSST.DP1 | phot.flux.density | Rubin:Simonyi | LSSTComCam | difference_image - g - CC_O_20241128_000181-R22_S11 2024-11-29T02:13:25.266980Z | POLYGON ICRS 37.607125 6.760401 37.386474 6.817755 37.329792 6.602476 37.550347 6.545141 | https://data.lsst.cloud/api/datalink/links?ID=ivo%3A%2F%2Forg.rubinobs%2Flsst-dp1%3Frepo%3Ddp1%26id%3D019b1989-f97d-7d2b-ba37-d6bd6b9679dd | application/x-votable+xml;content=datalink |
Print the date of the earliest (g-band) image for which bulk image stamps will be retrieved.
print('begin date:', scitab['t_max'][0], difftab['t_max'][0])
print('end date:', scitab['t_max'][len(scitab['t_max'])-1])
print('number of visits:', len(scitab['t_max']))
begin date: 60641.05308532407 60641.05308532407 end date: 60643.09300586805 number of visits: 9
4.2. Retrieve stamps¶
Retrieve the access_url from the query (or datalink) which points to the image location on the remote server for the first entries in each of the visit (scitab), template (reftab), and difference (difftab) image table results. The datalinks stored as dl_results below are in a format that can be used by the IVOA tools to generate the cutouts. The image products being retrieved for the science and difference images correspond to the earliest images.
dl_result_sci = DatalinkResults.from_result_url(scitab['access_url'][0], session=get_pyvo_auth())
dl_result_ref = DatalinkResults.from_result_url(reftab['access_url'][0], session=get_pyvo_auth())
dl_result_diff = DatalinkResults.from_result_url(difftab['access_url'][0], session=get_pyvo_auth())
Retrieve bulk image stamps with an edge size of fov (defined to be 0.004 degrees above). The get_cutout function returns both the image array (first parameter) and the corresponding Header Data Unit List (HDUL), a Python object used by astropy to contain one or more Header Data Units (HDUs) that contain astronomical image data, metadata, or tables.
sci, scihdul = get_cutout(dl_result_sci, tab_g['ra'][0], tab_g['dec'][0], get_pyvo_auth(), fov)
ref, refhdul = get_cutout(dl_result_ref, tab_g['ra'][0], tab_g['dec'][0], get_pyvo_auth(), fov)
diff, diffhdul = get_cutout(dl_result_diff, tab_g['ra'][0], tab_g['dec'][0], get_pyvo_auth(), fov)
First, visualize the cutouts of the science, template, and difference images for the earliest g-band images by displaying them side by side.
fig, ax = plt.subplots(1, 3, figsize=(8, 3))
plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(sci)
plt.title('science')
plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig)
display2.scale('linear', 'zscale')
display2.mtv(ref)
plt.title('template')
plt.sca(ax[2])
display3 = afwDisplay.Display(frame=fig)
display3.scale('linear', 'zscale')
display3.mtv(diff)
ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()
plt.title('difference')
plt.tight_layout()
fig.suptitle(str(mpcDesignation)+', MJD='+str(np.round(scitab['t_max'][0], 2)))
plt.show()
Figure 3: The cutout triplet for the first g-band detection of this Solar System object, Rubin discovery 2024 WP103.
4.3. Rotate visit and difference images¶
As Figure 3 shows, visit images are not aligned with the template image (North up, East left) since they are taken at many different telescope orientations. This section demonstrates how to align the visit image to match the orientation of the template and difference cutouts using the reproject python package. reproject requires the world coordinate system (WCS) of the existing and desired coordinate systems, and the image pixels of the image to reproject, in order to perform the calculation.
The cell below extracts these from the HUDL. The structure of the ImageF object returned by the cutout service has two HDUL extensions. The extension 0 contains a primary header of telescope metadata and no image pixels, while the extension 1 contains the science header and the science image pixels.
sci_primary_header = scihdul[0].header
sci_hdu = scihdul[1]
sci_header = scihdul[1].header
scidata = scihdul[1].data
sciwcs = WCS(scihdul[1].header)
ref_primary_header = refhdul[0].header
ref_hdu = refhdul[1]
ref_header = refhdul[1].header
refdata = refhdul[1].data
refwcs = WCS(refhdul[1].header)
diff_primary_header = diffhdul[0].header
diff_hdu = diffhdul[1]
diff_header = diffhdul[1].header
diffdata = diffhdul[1].data
diffwcs = WCS(diffhdul[1].header)
Use the reproject_interp function to put the science visit image onto the WCS of the template image.
reproj_scidata, footprint = reproject_interp((scidata, sciwcs), refwcs)
reproj_diffdata, footprint = reproject_interp((diffdata, diffwcs), refwcs)
Plot the reprojected science visit image and reprojected difference image together with the template image to compare.
fig = plt.figure(figsize=(20, 12))
ax1 = fig.add_subplot(1, 3, 1, projection=WCS(ref_hdu.header))
ax1.imshow(reproj_scidata, origin='lower', vmin=np.nanpercentile(reproj_scidata, 1),
vmax=np.nanpercentile(reproj_scidata, 99))
ax1.coords['ra'].set_axislabel('Right Ascension')
ax1.coords['dec'].set_axislabel('Declination')
ax1.set_title('Reprojected Visit Image')
ax2 = fig.add_subplot(1, 3, 2, projection=WCS(ref_hdu.header))
ax2.imshow(refdata, origin='lower', vmin=np.nanpercentile(refdata, 1),
vmax=np.nanpercentile(refdata, 99))
ax2.coords['ra'].set_axislabel('Right Ascension')
ax2.coords['dec'].set_axislabel('Declination')
ax2.coords['dec'].set_axislabel_position('r')
ax2.coords['dec'].set_ticklabel_position('r')
ax2.set_title('Template Image')
ax3 = fig.add_subplot(1, 3, 3, projection=WCS(ref_hdu.header))
ax3.imshow(reproj_diffdata, origin='lower', vmin=np.nanpercentile(reproj_diffdata, 1),
vmax=np.nanpercentile(reproj_diffdata, 99))
ax3.coords['ra'].set_axislabel('Right Ascension')
ax3.coords['dec'].set_axislabel('Declination')
ax3.coords['dec'].set_axislabel_position('r')
ax3.coords['dec'].set_ticklabel_position('r')
ax3.set_title('Reprojected Difference Image')
plt.subplots_adjust(wspace=1.0)
plt.show()
Figure 4: The science visit image (left) and difference image (right) after being reprojected onto the same WCS as the template image cutout (center).
5. Bulk image stamps¶
Bulk image stamps are a key visualization tool for Solar System objects. Below, retrieve many science visit image stamps and their corresponding template images at matching ra and dec coordinates to the target Solar System object. Save the image stamps as png files.
If >20 visits span the time during which the target was a DIA source, limit the number to visualize only the first 20 (to save time). The number of stamps can be limited by replacing the first line with num_visits = 20.
num_visits = len(scitab['t_max'])
print('number of visits:', num_visits)
for i in range(num_visits):
plt.title(f"MJD {np.round(scitab['t_max'][i], 4)}")
dl_result_sci = DatalinkResults.from_result_url(scitab['access_url'][i],
session=get_pyvo_auth())
sci, scihdul = get_cutout(dl_result_sci, tab_g['ra'][i], tab_g['dec'][i], get_pyvo_auth(), fov)
sci_hdu = scihdul[1]
sci_header = scihdul[1].header
scidata = scihdul[1].data
sciwcs = WCS(scihdul[1].header)
dl_result_ref = DatalinkResults.from_result_url(reftab['access_url'][i],
session=get_pyvo_auth())
ref, refhdul = get_cutout(dl_result_ref, tab_g['ra'][i], tab_g['dec'][i], get_pyvo_auth(), fov)
ref_hdu = refhdul[1]
ref_header = refhdul[1].header
refdata = refhdul[1].data
refwcs = WCS(refhdul[1].header)
reproj_scidata, footprint = reproject_interp((scidata, sciwcs), refwcs)
plt.imshow(reproj_scidata, origin='lower', vmin=np.nanpercentile(reproj_scidata, 1),
vmax=np.nanpercentile(reproj_scidata, 99))
figname = os.path.join(tempdir, 'stamp_' + str(i) + '.png')
if os.path.isfile(figname):
os.remove(figname)
plt.savefig(figname)
number of visits: 9
Figure 5: Final bulk science visit image saved to disk.
Visualize the bulk image stamps that were saved to disk.
fig, axes = make_subplot_grid(num_visits)
for i, ax in enumerate(axes):
j = i
figname = os.path.join(tempdir, 'stamp_' + str(j) + '.png')
img = mpimg.imread(figname)
ax.imshow(img)
plt.tight_layout()
fig.suptitle(str(mpcDesignation))
plt.subplots_adjust(top=0.95)
plt.show()
Figure 6: Image stamps from the science visit images of the Solar System object, Rubin discovery 2024 WP103 (
ssObjectId= 23133931615301936) in DP1, in chronological order, after being reprojected onto the same WCS as the template image.
6. Visualize cutouts with gif¶
Below, the cell will use the make_gif function defined in Section 1.2.3 to make a single gif movie out of all the pngs in the folder from the visit images in chronological order.
print(tempdir)
make_gif(tempdir)
dimg(data=open(tempdir+'/animation.gif', 'rb').read())
/deleted-sundays/sarahgreenstreet/dp1_308_5_temp
<IPython.core.display.Image object>
Figure 7: A gif animation generated from the reprojected visit image stamps in Figure 6.
Delete tempdir created to store cutouts and gif animation and all of its contents. NOTE: This permanently deletes the directory and all of its contents.
try:
shutil.rmtree(tempdir)
print(f"Directory '{tempdir}' and all its contents removed successfully.")
except FileNotFoundError:
print("The directory does not exist.")
Directory '/deleted-sundays/sarahgreenstreet/dp1_308_5_temp' and all its contents removed successfully.
7. Exercise for the learner¶
Generate r-band cutouts for 2009 AO58 (q=2.673640 au, e=0.148635, inc=7.449680 deg, and H=16.806999).