103. 3I/ATLAS image stamps#

103_3I_ATLAS_image_stamps

103. Image stamps for 3I/ATLAS (C/2025 N1)¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Chandler et al. (2026)
Container Size: Large
LSST Science Pipelines version: v29.2.0
Last verified to run: 2026-03-26
Repository: github.com/lsst/tutorial-notebooks\

Learning objective: To open and display the LSSTCam image stamps for observations of interstellar object 3I/ATLAS that were used in Chandler et al. (2026).

LSST data products: FITS files.

Packages: astropy

Credit: Cite Chandler et al. (2026) and DOI zenodo.org/records/19244429 if the FITS files are used in a publication. Tutorial developed by the Rubin Community Science team with input from Colin Orion Chandler. Please also reference the DOI for Rubin tutorials, 10.11578/rubin/dc.20250909.20, if this notebook is used for the preparation of, e.g., journal articles, software releases.

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¶

The third known interstellar object, 3I/ATLAS, was coincidentally observed by Rubin Observatory with LSSTCam during commissioning.

An analysis of commissioning images that include 3I/ATLAS is presented in Chandler et al. (2026), but these images have not yet been part of an official Rubin data release. Thus, a dedicated DOI for just this subset has been assigned: zenodo.org/records/19244429.

The purpose of this tutorial is to make image stamps containing 3I/ATLAS, plus their metadata and measurements, easily accessible via the Rubin Science Platform on a temporary basis, until these images (or at least those that pass image quality control and science validation) are released as part of Rubin's Data Preview 2.

Caveats:

Commissioning images may have (very) degraded image quality.
The pixel value units are in counts, not calibrated flux (nJy) as is standard for Rubin images.
Not all stamps are centered on 3I/ATLAS.
The photometry reported in Chandler et al. (2026) was performed on difference images, which are not available.

Related tutorials: See also the Solar System 200- and 300-level tutorials for Data Preview 1 (DP1).

1.1. Import packages¶

Import standard python packages for array manipulation (numpy) and plotting (matplotlib).

Import astropy modules to work with FITS files: fits, WCS, and SkyCoord (astropy documentation).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colormaps
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord

1.2. Define parameters¶

A copy of the 3I/ATLAS image stamps has been saved in a shared, read-only directory for easy access by all RSP users.

Set the path to the data directory.

In [2]:

path = "/rubin/cst_repos/tutorial-notebooks-data/data/3I_ATLAS_images_Chandler+26/"

2. Read in metadata files¶

Two files accompany the image stamp FITS files.

Read in each file, assert they contain data for 97 observations, and show the first row of each as an example.

2.1. Table 1 from Chandler et al. (2026)¶

Chandler+26_Table1.dat: Table 1 from Chander et al. (2026) containing observations and measured properties of 3I/ATLAS.

Row: integer row number
Date, Midpoint: UTC TAI time for the midpoint of the observation
VisitID, Chip: integer visit identifier and LSSTCam detector number
RA, Dec: measured coordinates of 3I/ATLAS (deg)
eRA, eDec: error in measured coordinates (arcsec)
Band: LSST filter (u, g, r, i, z, or y)
Mag, eMag: apparent magnitude and uncertainty
Phase: solar angle phase (deg)
r: heliocentric distance (au)

In [3]:

fn_table1 = "Chandler+26_Table1.dat"
table1 = np.genfromtxt(path + fn_table1, delimiter='', dtype=None,
                       names=True, filling_values=np.nan)
assert len(table1) == 97
print(table1[0])

(1, '2025-06-21', '08:11:32', 2025062000620, '104', 276.5133, -18.755701, 0.043, 0.05, 'z', nan, nan, 1.63, 4.837)

2.2. Image stamp metadata¶

image_stamp_metadata.dat: Metadata for the image stamp FITS files.

filename: name of the FITS file for the image stamp
centered: True if the stamp is centered on the NASA JPL Horizons coordinates
visit, detector: integer visit identifier and LSSTCam detector number
band: LSST filter (u, g, r, i, z, or y)
date_mid_tai: UTC TAI time for the midpoint of the observation
mjd_mid_tai: MJD time for the midpoint of the observation
exp_time: exposure time (s)
rot_pa: rotation position angle (deg)
ra_horizons, dec_horizons: NASA JPL Horizons coordinates (deg)

In [4]:

fn_metadata = "image_stamp_metadata.dat"
metadata = np.genfromtxt(path + fn_metadata, delimiter='', dtype=None,
                         names=True, filling_values=np.nan)
assert len(metadata) == 97
print(metadata[0])

('cutout-2025062000620-104-z-ctr.fits', True, 2025062000620, 104, 'z', '2025-06-21T08:11:32.665', 60847.34135, 30.0, 30.09, 276.51285, -18.7557)

Plot the sky coordinates of 3I/ATLAS from NASA JPL Horizons for all Rubin Observations.

In [5]:

cmap = colormaps['viridis']
fig, ax = plt.subplots(1, 1, figsize=(7, 3))
im = ax.scatter(metadata['ra_horizons'], metadata['dec_horizons'],
                c=metadata['mjd_mid_tai'], cmap=cmap)
ax.set_xlabel('RA (deg)')
ax.set_ylabel('Dec (deg)')
fig.suptitle("JPL Horizons coordinates of 3I/ATLAS for Rubin observations.")
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])
fig.colorbar(im, cax=cbar_ax, label='MJD TAI')
fig.show()

No description has been provided for this image

Figure 1: A plot of the NASA JPL Horizons coordinates, RA vs. Dec, for the 97 Rubin observations. Points are colored by the MJD of the observation, as illustrated by the colorbar.

3. Image stamp FITS files¶

Choose an image to examine and display. For this tutorial the 6th image (index = 5) was selected as a demonstration.

In [6]:

img_i = 5

3.1. Read the pixel and header data¶

Read the image pixel data and header.

In [7]:

img_data, img_hdr = fits.getdata(path + metadata['filename'][img_i], header=True)

Option to display the header, which is limited to pixel data and World Coordinate System (WCS) information.

In [8]:

# img_hdr

Notice that pixel data is in counts (adu), not calibrated flux.

In [9]:

print(img_hdr.get('BUNIT'))

adu

3.2. Get the image WCS¶

Get the WCS from the header.

In [10]:

img_wcs = WCS(img_hdr)

Display the WCS parameters.

In [11]:

img_wcs

Out[11]:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 271.404313301785 -17.445035526301
CRPIX : 20177.165592 -9863.563196
CD1_1 CD1_2  : 2.38053181271285e-05 5.01925062832633e-05
CD2_1 CD2_2  : 5.01446130282366e-05 -2.37681232783881e-05
NAXIS : 600  600

Above, notice that the CRVAL coordinates (RA, Dec in degrees) do not correspond to pixels CRPIX 1, 1. This is because the WCS is inherited from the larger image from which this stamp was extracted.

Convert the Horizons coordinates to the astropy SkyCoord format, and use the footprint_contains method of the WCS to confirm that the coordinates are within the image stamp. Use the world_to_pixel method to convert the sky coordinates to pixel coordinates. Display the sky and pixel coordinates.

In [12]:

coord_horizons = SkyCoord(metadata['ra_horizons'][img_i], metadata['dec_horizons'][img_i],
                          frame="icrs", unit="deg")
assert img_wcs.footprint_contains(coord_horizons)
xy_horizons = img_wcs.world_to_pixel(coord_horizons)
print(coord_horizons)
print(xy_horizons)

<SkyCoord (ICRS): (ra, dec) in deg
    (271.44738, -18.68487)>
(array(300.52497177), array(375.06040619))

Repeat the above with the measured coordinates from Table 1 (if they exist).

In [13]:

if np.isfinite(table1['RA'][img_i]):
    coord_table1 = SkyCoord(table1['RA'][img_i], table1['Dec'][img_i], frame="icrs", unit="deg")
    assert img_wcs.footprint_contains(coord_table1)
    print(coord_table1)
    xy_table1 = img_wcs.world_to_pixel(coord_table1)
    print(xy_table1)
else:
    print('No measured coordinates for image index value', img_i)

<SkyCoord (ICRS): (ra, dec) in deg
    (271.447823, -18.684875)>
(array(303.67821129), array(381.92766861))

3.3. Display the image¶

Prepare to display the image.

Define the minimum and maximum pixel values to define black and white in the greyscale to represent pixel flux, and display the values.
Define the image extent in image pixel coordinates, and display the values.
Display the rotation of LSSTCam at the time the image was obtained.
Display whether the stamp is centered on 3I/ATLAS.

In [14]:

scale_min = np.mean(img_data) - 0.5*np.mean(img_data)
scale_max = np.mean(img_data) + 0.5*np.mean(img_data)
print('Greyscale min,max (adu): ', scale_min, scale_max)

Greyscale min,max (adu):  2593.078 7779.2334

In [15]:

img_extent = [img_hdr.get('CRPIX1'), img_hdr.get('CRPIX1') + img_hdr.get('NAXIS1'),
              img_hdr.get('CRPIX2'), img_hdr.get('CRPIX2') + img_hdr.get('NAXIS2')]
print('Image extent (x1, x2, y1, y2): ', img_extent)

Image extent (x1, x2, y1, y2):  [20177.165592, 20777.165592, -9863.563196, -9263.563196]

In [16]:

print('Rotation position angle (deg E of N): ', metadata['rot_pa'][img_i])
if (metadata['rot_pa'][img_i] > 90) & (metadata['rot_pa'][img_i] < 270):
    print(' > Note: the image will be displayed with "Dec" on the x-axis.')

Rotation position angle (deg E of N):  115.36
 > Note: the image will be displayed with "Dec" on the x-axis.

In [17]:

print('The stamp is centered on 3I/ATLAS? ', metadata['centered'][img_i])

The stamp is centered on 3I/ATLAS?  False

Show the image and mark the NASA JPL Horizons coordinates, the measured coordinates of 3I/ATLAS (if they exist), and the image stamp center.

In [18]:

fig, ax = plt.subplots(subplot_kw=dict(projection=img_wcs))
ax.imshow(img_data, cmap='gray', vmin=scale_min, vmax=scale_max,
          origin='lower', extent=img_extent)
ax.grid(color='white', ls='solid')
ax.plot(xy_horizons[0] + img_hdr.get('CRPIX1'),
        xy_horizons[1] + img_hdr.get('CRPIX2'),
        's', ms=8, mec='cyan', color="None")
if np.isfinite(table1['RA'][img_i]):
    ax.plot(xy_table1[0] + img_hdr.get('CRPIX1'),
            xy_table1[1] + img_hdr.get('CRPIX2'),
            'o', ms=8, mec='magenta', color="None")
ax.plot(301 + img_hdr.get('CRPIX1'),
        301 + img_hdr.get('CRPIX2'),
        '+', ms=8, mec='yellow', color="None")
fig.show()

Figure 2: The full image stamp, with white grid lines to show the lines of RA and Dec. Markers are displayed for the stamp center (yellow plus), the NASA JPL Horizons coordinates (cyan square), and the measured coordinates from Table 1 of Chandler et al. (2026; magenta circle).

Re-display a subset of the image centered on the NASA JPL Horizons coordinates.

In [19]:

fig, ax = plt.subplots(subplot_kw=dict(projection=img_wcs))
ax.imshow(img_data, cmap='gray', vmin=scale_min, vmax=scale_max,
          origin='lower', extent=img_extent)
ax.grid(color='white', ls='solid')
ax.plot(xy_horizons[0] + img_hdr.get('CRPIX1'),
        xy_horizons[1] + img_hdr.get('CRPIX2'),
        's', ms=8, mec='cyan', color="None")
if np.isfinite(table1['RA'][img_i]):
    ax.plot(xy_table1[0] + img_hdr.get('CRPIX1'),
            xy_table1[1] + img_hdr.get('CRPIX2'),
            'o', ms=8, mec='magenta', color="None")
ax.plot(301 + img_hdr.get('CRPIX1'),
        301 + img_hdr.get('CRPIX2'),
        '+', ms=8, mec='yellow', color="None")
halfsize = 50
ax.set_xlim([xy_horizons[0] + img_hdr.get('CRPIX1') - halfsize,
             xy_horizons[0] + img_hdr.get('CRPIX1') + halfsize])
ax.set_ylim([xy_horizons[1] + img_hdr.get('CRPIX2') - halfsize,
             xy_horizons[1] + img_hdr.get('CRPIX2') + halfsize])
fig.show()

Figure 3: A zoom-in on Figure 2 at the location of the NASA JPL Horizons coordinates (cyan square), showing also the measured coordinates from Table 1 of Chandler et al. (2026; magenta circle).