103.12. Reading and writing images#
103.12. Reading and writing images¶
For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: r30.0.11
Last verified to run: 2026-10-27
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20
Learning objective: How to write a deep_coadd image to a FITS file and read it back, using both the LSST Science Pipelines tools and astropy.
LSST data products: deep_coadd
Packages: lsst.daf.butler, lsst.images, astropy.io.fits
Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.
Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.
1. Introduction¶
A deep_coadd is a CellCoadd object with two ways to save it to a FITS file: .image.write(), which saves only the image plane, and .write(), which saves the full object, including the mask, variance, point spread function (PSF), background models, and provenance.
Saving images to disk is useful for sharing a cutout, working offline, or using tools outside the Rubin Science Platform (RSP). There are also two ways to read a saved FITS file back: with the LSST Science Pipelines' read_archive function, which reconstructs the original, full-featured object, or with astropy, which returns the raw pixel data and FITS header only.
The same .write(), .image.write(), and read_archive() methods work the same way for visit_image and difference_image, which are also lsst.images objects.
Related tutorials: Tutorial 202.1 demonstrates deep_coadd data access and attributes in more detail.
1.1. Import packages¶
Import os and getpass to set up a scratch directory for saved files.
From the LSST Science Pipelines import the Butler, and read_archive to read a saved FITS file back as a native LSST object.
Import astropy.io.fits to read a saved FITS file with astropy, and numpy to compare pixel arrays.
import os
import getpass
from lsst.daf.butler import Butler
from lsst.images.serialization import read_archive
from astropy.io import fits
import numpy as np
1.2. Define parameters and functions¶
Instantiate the Butler, and set up a scratch directory to save files to: RSP users should save files to the shared scratch directory, /deleted-sundays/, in a sub-directory named with their username.
butler = Butler('dp2', collections=['dp2'])
username = getpass.getuser()
userdir = os.path.join('/deleted-sundays', username)
os.makedirs(userdir, exist_ok=True)
2. Get a deep_coadd¶
Query for a deep_coadd image overlapping coordinates RA, Dec near the center of the ECDFS field, obtained with the $r$-band filter, and retrieve it.
See tutorial 202.1 for the full range of methods to query for and retrieve deep_coadd images.
ra = 53.076
dec = -28.110
band = 'r'
dataset_refs = butler.query_datasets('deep_coadd',
where="band.name = :band AND \
patch.region OVERLAPS POINT(:ra, :dec)",
bind={'band': band, 'ra': ra, 'dec': dec})
deep_coadd = butler.get(dataset_refs[0])
deep_coadd
CellCoadd
├── skymap: lsst_cells_v2
├── tract: 5063
├── patch: id=15, index=(y=1, x=5), inner_bbox=[y=3000:6000, x=15000:18000], cells=([y=2850:6150,
│ x=14850:18150], cell_shape=(y=150, x=150))
├── band: r
├── bbox: [y=2850:6150, x=14850:18150]
├── mask_fractions: rejected
├── noise_realizations: 1 image (dtype float32)
├── aperture_corrections: 64 fields
├── image: (dtype float32)
├── mask (Mask)
│ └── schema (MaskSchema)
│ ├── dtype: uint8
│ └── Mask planes
│ ┏━━━━━┳━━━━━━━┳━━━━━━┳━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┓
│ ┃ Bit ┃ Index ┃ Mask ┃ Name ┃ Description ┃ Set pixels ┃
│ ┡━━━━━╇━━━━━━━╇━━━━━━╇━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━┩
│ │ 0 │ 0 │ 0x1 │ NO_DATA │ No data was available for this │ 25321 │
│ │ │ │ │ │ pixel. │ │
│ │ 1 │ 0 │ 0x2 │ INTERPOLATED │ Pixel value is the result of │ 32543 │
│ │ │ │ │ │ interpolating nearby good pixels. │ │
│ │ 2 │ 0 │ 0x4 │ COSMIC_RAY │ A cosmic ray affected this pixel │ 7222 │
│ │ │ │ │ │ on at least one input image (and │ │
│ │ │ │ │ │ was interpolated). │ │
│ │ 3 │ 0 │ 0x8 │ SATURATED │ More than 10% of the potential │ 55499 │
│ │ │ │ │ │ input visits had a saturated │ │
│ │ │ │ │ │ pixel at this location │ │
│ │ │ │ │ │ ('potential' because saturated │ │
│ │ │ │ │ │ pixel values are not actually │ │
│ │ │ │ │ │ propagated to the coadd). │ │
│ │ │ │ │ │ SATURATED always implies │ │
│ │ │ │ │ │ REJECTED, and is often a reason │ │
│ │ │ │ │ │ for NO_DATA. │ │
│ │ 4 │ 0 │ 0x10 │ DETECTION_EDGE │ Pixel was too close to the edge │ 92204 │
│ │ │ │ │ │ of the patch to be considered for │ │
│ │ │ │ │ │ detection, due to the finite size │ │
│ │ │ │ │ │ of the detection kernel. │ │
│ │ 5 │ 0 │ 0x20 │ CLIPPED │ Region was identified as a │ 817 │
│ │ │ │ │ │ probable artifact when comparing │ │
│ │ │ │ │ │ multiple single-visit warps. │ │
│ │ │ │ │ │ CLIPPED always implies REJECTED. │ │
│ │ 6 │ 0 │ 0x40 │ REJECTED │ At least one input visit was left │ 492810 │
│ │ │ │ │ │ out of the coadd for this pixel │ │
│ │ │ │ │ │ due to masking. REJECTED always │ │
│ │ │ │ │ │ implies INEXACT_PSF. │ │
│ │ 7 │ 0 │ 0x80 │ DETECTED │ Pixel was part of a detected │ 2108469 │
│ │ │ │ │ │ source. │ │
│ │ 8 │ 1 │ 0x1 │ INEXACT_PSF │ The set of visits contributing to │ 493568 │
│ │ │ │ │ │ this pixel differs from the set │ │
│ │ │ │ │ │ of visits contributing to the PSF │ │
│ │ │ │ │ │ model for its cell. │ │
│ └─────┴───────┴──────┴────────────────┴───────────────────────────────────┴────────────┘
├── variance: (dtype float32)
├── sky_projection (ICRS coordinates)
│ ├── center pixel: (x=16499.5, y=4499.5) → 03h31m58.3s -28d06m14s (RA 52.993059°, Dec
│ │ -28.103914°)
│ ├── Image extent: 11 x 11 arcmin @ 0 deg E of N
│ ├── Nominal pixel scale: 0.20 arcsec
│ ├── fits_wcs: available
│ └── Corners
│ ┏━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━━┓
│ ┃ x ┃ y ┃ RA ┃ Dec ┃ RA (°) ┃ Dec (°) ┃
│ ┡━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━━┩
│ │ 14849.5 │ 2849.5 │ 03h32m23.3s │ -28d11m44s │ 53.096981 │ -28.195602 │
│ │ 14849.5 │ 6149.5 │ 03h32m23.3s │ -28d00m44s │ 53.096965 │ -28.012288 │
│ │ 18149.5 │ 6149.5 │ 03h31m33.4s │ -28d00m44s │ 52.889312 │ -28.012146 │
│ │ 18149.5 │ 2849.5 │ 03h31m33.4s │ -28d11m44s │ 52.888979 │ -28.195459 │
│ └─────────┴────────┴─────────────┴────────────┴───────────┴────────────┘
├── psf (CellPointSpreadFunction)
│ ├── bounds: [y=2850:6150, x=14850:18150] in grid ([y=2850:6150, x=14850:18150],
│ │ cell_shape=(y=150, x=150))
│ └── kernel_bbox: [y=-17:18, x=-17:18]
├── provenance (CoaddProvenance)
│ ├── instrument: LSSTCam
│ ├── physical_filter: r_57
│ ├── input images: 98 from 31 visits
│ ├── day_obs: 20250728 - 20260105
│ ├── cells: 484 of 484 with contributions
│ └── per cell: 13 - 28 input images (median 22)
└── backgrounds: pretty, object (subtracted)
image_filename = os.path.join(userdir, 'deep_coadd_image.fits')
deep_coadd.image.write(image_filename)
print('File size:', os.path.getsize(image_filename), 'bytes')
File size: 39245760 bytes
3.2. Write the full object¶
Use .write() to save the full CellCoadd object, including the mask, variance, PSF, background models, and provenance, as a FITS file.
full_filename = os.path.join(userdir, 'deep_coadd_full.fits')
deep_coadd.write(full_filename)
print('File size:', os.path.getsize(full_filename), 'bytes')
File size: 110111040 bytes
The full object is substantially larger. Use astropy's fits.info to compare the FITS extensions in each file.
fits.info(image_filename)
print('=' * 55)
fits.info(full_filename)
Filename: /deleted-sundays/gmartin/deep_coadd_image.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 11 () 1 IMAGE 1 CompImageHDU 35 (3300, 3300) float32 2 JSON 1 BinTableHDU 12 1R x 1C [PB(11316)] 3 INDEX 1 BinTableHDU 24 3R x 7C [5A, J, 8A, L, K, K, K]
======================================================= Filename: /deleted-sundays/gmartin/deep_coadd_full.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 14 () 1 IMAGE 1 CompImageHDU 35 (3300, 3300) float32 2 MASK 1 CompImageHDU 61 (3300, 3300) int32 3 VARIANCE 1 CompImageHDU 35 (3300, 3300) float32 4 MASK_FRACTIONS/REJECTED 1 CompImageHDU 17 (3300, 3300) float32 5 NOISE_REALIZATIONS/0 1 CompImageHDU 17 (3300, 3300) float32 6 PSF 1 CompImageHDU 11 (35, 35, 22, 22) float64 7 APERTURE_CORRECTIONS/TABLE 1 BinTableHDU 142 484R x 66C [K, K, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D] 8 PROVENANCE/INPUTS 1 BinTableHDU 27 98R x 6C [K, I, J, B, 7D, 7D] 9 PROVENANCE/CONTRIBUTIONS 1 BinTableHDU 41 10387R x 12C [I, I, K, I, L, D, D, D, D, D, D, L] 10 BACKGROUND/FIELDS/PRETTY/DATA 1 CompImageHDU 9 (52, 52) float32 11 BACKGROUND/FIELDS/OBJECT/DATA 1 CompImageHDU 9 (26, 26) float32 12 JSON 1 BinTableHDU 12 1R x 1C [PB(32241)] 13 INDEX 1 BinTableHDU 24 13R x 7C [29A, J, 8A, L, K, K, K]
The image-only file still has JSON and INDEX extensions in addition to IMAGE -- these store metadata needed by the LSST Science Pipelines tools to reconstruct the object. The full object file additionally has mask, variance, PSF, background, and provenance extensions.
3.3. Download a file locally¶
The shared scratch directory is not visible in the JupyterLab file browser, so a file saved there cannot be downloaded directly. Rather than saving a second copy into the home directory, create a temporary symbolic link to the file from the home directory; it will then appear in the file browser, where it can be downloaded via the right-click menu.
home_symlink = os.path.join(os.path.expanduser('~'), 'deep_coadd_image.fits')
os.symlink(image_filename, home_symlink)
print('Symlinked to', home_symlink)
Symlinked to /home/gmartin/deep_coadd_image.fits
4. Read deep_coadd back¶
4.1. Read with LSST Science Pipelines tools¶
Use read_archive to read a saved FITS file back as a native LSST object. This gives a richer result than reading with astropy: the object retains its mask, variance, PSF, and other components, and all of its methods.
For example, sky_projection still converts between pixel and sky coordinates, psf.compute_kernel_image() still computes a PSF image at a given position, and apply_background() still works, exactly as they would on the original deep_coadd. This is useful because the read-back object can be used directly in any further analysis that relies on these methods, without needing to re-fetch the original from the Butler. Reading the same file with astropy only returns a plain pixel array and a FITS header.
coadd_read = read_archive(full_filename)
print(type(coadd_read))
print('Image data is unchanged:', np.array_equal(coadd_read.image.array, deep_coadd.image.array))
<class 'lsst.images.cells._coadd.CellCoadd'> Image data is unchanged: True
4.2. Read with astropy¶
Use astropy.io.fits to read the saved FITS file. Only the raw pixel data and FITS header are returned, and the data are in the extension named IMAGE, not the (empty) primary extension.
For deep_coadd, the WCS in this header is exact. For visit_image and difference_image, the FITS WCS is only an approximation -- good enough for display or object finding, not for precision astrometry -- so use the LSST Science Pipelines object's own sky_projection methods instead (via read_archive in Section 4.1), as described in the DP1 known issues documentation.
with fits.open(image_filename) as hdul:
astropy_image = hdul['IMAGE'].data
print(type(astropy_image))
print('Image data is unchanged:', np.array_equal(astropy_image, deep_coadd.image.array))
<class 'numpy.ndarray'> Image data is unchanged: True
Clean up, including removing images from the scratch space and symbolic links from the home directory.
os.remove(home_symlink)
os.remove(image_filename)
os.remove(full_filename)
del dataset_refs, deep_coadd, coadd_read, astropy_image
del username, userdir, image_filename, full_filename, home_symlink