EuclidRust – Rust FITS reader for Euclid, 10x faster than astropy

High-performance FITS I/O for Python, written in Rust and exposed via PyO3 as a normal, pip install-able Python extension module. Built as a testbed for data-processing tools for the ESA Euclid mission, but useful for any FITS-based workflow. The Python package it contains, and the name you import, is euclidrust.

import euclidrust

stats = euclidrust.fits_stats("tile.fits")
cutout = euclidrust.extract_cutout("tile.fits", ra_deg=83.82, dec_deg=-5.39, width=50, height=50)
euclidrust.export_cutout_parquet("tile.fits", 83.82, -5.39, 50, 50, "cutout.parquet")

About this project. The code, benchmarks and documentation here were written by Claude Code (Anthropic's coding agent). The author decided what to build, in what order, how to measure it and when to stop, and steered the agent throughout. Treat the project as experimental: read the code, run the benchmark yourself, and judge the numbers rather than the author.

Contents: What files it works with · How Python and Rust fit together · Install · Use it · Many images at once · NVIDIA GPU · Tests · Run the benchmark · Troubleshooting · Motivation · Benchmark results · Limitations · Roadmap

Short answer: 2-D FITS images, raw or tile-compressed (Rice / GZIP), whether the image sits in the primary HDU or in an extension.

FITS (Flexible Image Transport System) is the standard file format of astronomy. A FITS file is a sequence of blocks called HDUs ("header/data units"): each has a text header and a payload, which can be an image or a table. euclidrust reads images: a rectangular grid of pixels (a sky mosaic, a detector frame...), plus the header keywords that say how pixels map to the sky (the WCS).

Where the image sits depends on the product. A Euclid MER mosaic keeps its image in the first HDU (the "primary HDU"), which is what euclidrust reads by default. Calibrated frames (VIS, NISP) and raw frames have an empty primary HDU and one image extension per detector, with names such as DET11.SCI or 1-1.E. For those you pick the HDU with the hdu= argument (an index or a name), and list_hdus(path) shows what a file contains.

Supported:

Requirement	Details
An image HDU	2-D ( NAXIS = 2 ). The primary HDU by default, or any image extension throughhdu= (0-based index orEXTNAME )
Compression	None, or FITS tile compression with RICE_1 ,GZIP_1 orGZIP_2 (the.fits.fz format written byfpack and by astropy'sCompImageHDU ), including quantized floating-point images (NO_DITHER ,SUBTRACTIVE_DITHER_1 ,SUBTRACTIVE_DITHER_2 ). A compressed image always lives in an extension, so passhdu=1 (or its name)
Pixel type ( BITPIX )	8 ,16 ,32 ,64 (integers),-32 ,-64 (floating point).BZERO /BSCALE are applied, so unsigned 16-bit raw data (BZERO = 32768 ) works. 64-bit integers, signed or unsigned (FITS stores unsigned 64-bit as signed withBZERO = 2^63 ), are handledexactly , in raw and in gzip-compressed images (cfitsio compresses 64-bit integers with gzip only), including a whole-numberBSCALE /BZERO : the exact values are inmin_int /max_int (fits_stats ) anddata_int (extract_cutout ), as Python ints, and the Parquetvalue column isInt64 /UInt64 ; the ordinary float fields hold the nearestfloat64
WCS, only needed for cutouts	CTYPE1/2 =RA---xxx /DEC--xxx withxxx one ofTAN ,SIN ,ARC ,STG ,ZEA ,ZPN (the zenithal projections),TAN-SIP (tangent plus the SIP polynomial distortion) orTPV (tangent plus thePV1_m /PV2_m polynomials used by the Euclid VIS pipeline), together withCRVAL1/2 ,CRPIX1/2 and a linear matrix written asCD (CD1_1 ,CD1_2 ,CD2_1 ,CD2_2 ), asCDELT1/2 with an optionalPC matrix, or asCDELT1/2 with the legacyCROTA2 .LONPOLE is honoured.fits_stats needs no WCS

Not supported. Anything else is rejected with a ValueError that says why, rather than being read wrongly (these cases are covered by tests/test_input_validation.py):

Kind of file	Error you get
Tile compression with another algorithm ( HCOMPRESS_1 ,PLIO_1 )	HDU 1 (SCI): tile-compressed image uses ZCMPTYPE = HCOMPRESS_1 which is not supported yet (supported: RICE_1, GZIP_1, GZIP_2)
An extension that is a table (catalogue, spectrum...)	HDU 1 is a BINTABLE table, not an image
3-D cubes ( NAXIS = 3 )	NAXIS = 3 is not supported ...: only 2-D images are handled
Not a FITS file	truncated file: no END card found
A projection that is not implemented (cylindrical or all-sky ones such as CAR ,AIT ,MOL ,HPX ...), a mismatchedCTYPE pair, or a missingCTYPE	Cutouts raise unsupported WCS projection: found CTYPE1/CTYPE2 = ...; supported are RA---TAN/DEC--TAN (also with -SIP distortion), RA---TPV/DEC--TPV and the zenithal projections SIN, ARC, STG, ZEA and ZPN ... (fits_stats still works)
Distortion stored in lookup tables (HST CPDIS /D2IMDIS keywords)	Cutouts raise this header declares distortion corrections stored in lookup tables ... which euclidrust does not apply (fits_stats still works)
No WCS at all (no CRVAL /CRPIX / matrix keywords)	Cutouts raise this file has no WCS header (CRVAL, CRPIX and either a CD matrix or CDELT/PC) (fits_stats still works)

Null pixels. NaN pixels of floating-point images are left out of min, max and mean and counted in n_nan (infinities are values and stay in). The BLANK null value of integer images is not honoured: those pixels are counted as ordinary values.

If you forget hdu= on a file whose primary HDU is empty, the error lists the image HDUs it found and tells you what to pass, for example [1] DET11.SCI (2040x2040, BITPIX -32) ... Pass hdu=<index or EXTNAME>.

What this means for Euclid Q1 data. I read the headers of a handful of files in each product family of the public Q1 release (this is a sample, not a full survey of the archive) and ran the code on real files from four of them, comparing with astropy:

Product	Where the image is	Status
MER mosaics ( EUC_MER_BGSUB-MOSAIC-... , uncompressed)	primary HDU	Works (stats and cutouts). All benchmark numbers come from the NISP Y-band tile. The other bands and theRMS /FLAG mosaics have the same layout but were not run
NISP calibrated frames ( EUC_NIR_W-CAL-IMAGE-BKG... )	extensions DETnn.SCI ,.RMS ,.DQ (TAN WCS)	Works withhdu="DET11.SCI" etc. Checked on one detector: stats and a cutout identical to astropy
VIS calibrated frames ( EUC_VIS_SWL-BKG... )	extensions 1-1.E ... (TPV WCS)	Works withhdu="1-1.E" etc. Checked on one detector: stats identical to astropy; over 300 random sky positions the cutout center agrees with astropy/wcslib to 3 × 10^-8 pixel, and all cutout pixels are identical
RAW frames, NISP ( EUC_LE1_NISP... )	extensions DETnn.SCI (unsigned 16-bit viaBZERO )	Stats work (checked on one detector). Cutouts are refused: the WCS is in detector-mosaic coordinates (Y_MOSA ), not sky coordinates
RAW frames, VIS ( EUC_LE1_VIS... )	extensions 1-1.E ...	Same layout as the other frames, not run
MER segmentation maps ( EUC_MER_FINAL-SEGMAP... )	primary HDU, BITPIX = 64	Stats work (checked on the real 2.9 GB file: min and max identical to astropy, exact inmin_int /max_int , mean agreeing to ~13 digits). NoCD matrix in the header, so no cutouts
Catalogues, PSF files, other tables	table extensions	Not images, out of scope

None of the Q1 image files I looked at were tile-compressed, so there was no compressed Euclid file to test on. Compressed-image support was validated on copies of real Euclid data that I compressed with astropy/cfitsio (see Compressed images).

To check a file yourself before running anything on it, list its HDUs:

import euclidrust
for hdu in euclidrust.list_hdus("file.fits"):
    print(hdu)

or look at the header keywords with astropy:

from astropy.io import fits

h = fits.getheader("file.fits", 1)          # 0 = primary HDU, 1 = first extension...
print("NAXIS", h["NAXIS"], "| BITPIX", h["BITPIX"])
print("CTYPE", h.get("CTYPE1"), h.get("CTYPE2"))
print("linear matrix keywords:", [k for k in ("CD1_1", "CDELT1", "PC1_1", "CROTA2") if k in h])

You want NAXIS 2, a BITPIX from the table above and, for cutouts, a CTYPE pair from the projection list above and some linear-matrix keywords (CD1_1, or CDELT1).

If you are an astronomer who knows Python but has never touched Rust, this is the one thing to understand: you never write or run Rust yourself. You install the package and import it like numpy or astropy.

Under the hood there are three pieces:

Piece	Where	What it does
The Rust code	euclidrust/src/lib.rs	Does the actual work: memory-maps the FITS file, runs the parallel computation, converts sky coordinates to pixels
PyO3	a dependency in euclidrust/Cargo.toml	A Rust library that turns Rust functions and structs into Python functions and objects, and converts arguments between the two (for example, hdu= accepts a Pythonint orstr , which PyO3 maps to a Rustenum )
maturin	build-backend ineuclidrust/pyproject.toml	Compiles the Rust into a shared library and installs it in your Python environment ( pip calls it for you)

PyO3 shows up in lib.rs as a handful of annotations. Here is the whole mapping, taken from the real code:

#[derive(FromPyObject)]                         // -> accepts a Python int OR str
enum HduSelector { Index(usize), Name(String) }

#[pyfunction]                                   // -> a Python function
fn fits_stats(path: &str, hdu: Option<HduSelector>) -> PyResult<FitsStats> { /* ... */ }

#[pyclass]                                      // -> a Python class
struct FitsStats {
    #[pyo3(get)] min: f64,                      // -> readable attribute: stats.min
    #[pyo3(get)] max: f64,
    #[pyo3(get)] mean: f64,
    /* ... */
}

#[pymodule]                                     // -> `import euclidrust`
fn euclidrust(_py: Python<'_>, m: &Bound<'_, PyModule>) -> PyResult<()> {
    m.add_function(wrap_pyfunction!(fits_stats, m)?)?;
    /* ... */
    Ok(())
}

and on the Python side it is just:

import euclidrust                       # loads the compiled Rust library
s = euclidrust.fits_stats("tile.fits")  # calls the Rust function
print(s.min, s.max, s.mean)              # reads the fields of the Rust struct

So the build chain is: Rust source -> (maturin, run by pip) -> compiled library such as euclidrust.cpython-313-x86_64-linux-gnu.so (.pyd on Windows) -> installed into your virtual environment -> import euclidrust. The scripts in scripts/ contain no Rust and no PyO3 code; they only import euclidrust. Errors raised in Rust arrive in Python as ordinary ValueError exceptions.

You need three things on your machine:

• Python 3.8 or newer (developed and tested with 3.13).
• A Rust toolchain, version 1.88 or newer (the Arrow / Parquet dependencies require it; developed with 1.100). Install it withrustup :
  • Linux / macOS: curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
  • Windows: download rustup-init.exe fromhttps://rustup.rs/ . It will ask you to install the Visual Studio C++ Build Tools; accept.
  • Then open a new terminal and check withrustc --version .
• Linux / macOS:
• A C linker , which is normally already there (build-essential on Debian/Ubuntu, Xcode command-line tools on macOS, the Build Tools on Windows).

Everything else (maturin, PyO3 and the other Rust libraries) is downloaded automatically during the build, so the first install needs internet access.

Tested on Debian 13 (x86_64) with Python 3.13 and Rust 1.100, on an Intel Core i7-8559U (4 cores / 8 threads) with 31 GiB of RAM. macOS and Windows have not been tested; if you try them, please open an issue with what happened.

git clone https://github.com/MaxIKweeger/EuclidRust.git
cd EuclidRust

python3 -m venv .venv                 # a private Python environment
source .venv/bin/activate             # Windows PowerShell: .venv\Scripts\Activate.ps1

pip install ./euclidrust             # compiles the Rust code (a few minutes); GPU support is included

The first compile takes a few minutes because the release profile enables full link-time optimisation; later rebuilds are faster. When it finishes, check that it worked:

python -c "import euclidrust; print('OK, threads available:', euclidrust.detected_cpus())"

If that prints OK, threads available: <a number>, you are done.

Developing the Rust code? Use pip install maturin, then run maturin develop --release inside euclidrust/. That installs the module in editable mode and rebuilds it after you change src/lib.rs.

You need a FITS file of the kind described in What kind of files does it work with?. Euclid MER mosaics are of this kind; see the benchmark section for where to get one.

import numpy as np
import euclidrust

path = "tile.fits"

for hdu in euclidrust.list_hdus(path):
    print(hdu)

s = euclidrust.fits_stats(path)
print(s.min, s.max, s.mean)          # physical values (BZERO/BSCALE applied)
print(s.naxis1, s.naxis2, s.bitpix)  # image shape and data type
print(s.threads_used, s.elapsed_ms)  # how it ran

c = euclidrust.extract_cutout(path, ra_deg=273.0716597, dec_deg=68.0,
                               width=50, height=50)
image = np.array(c.data).reshape(c.height, c.width)   # row-major -> 2-D array
print(c.pixel_x_center, c.pixel_y_center)             # FITS 1-indexed pixel position

r = euclidrust.export_cutout_parquet(path, 273.0716597, 68.0, 50, 50, "cutout.parquet")
print(r.rows, r.out_path)

import pandas as pd                  # pip install pandas pyarrow
print(pd.read_parquet("cutout.parquet").head())       # columns: x, y, value

Files whose images are in extensions (Euclid calibrated and raw frames) need the hdu argument, a name or a 0-based index, on all three functions:

frame = "EUC_NIR_W-CAL-IMAGE-BKG_H-2681-0_20240930T175032.222263Z.fits"

s = euclidrust.fits_stats(frame, hdu="DET11.SCI")              # by EXTNAME
s = euclidrust.fits_stats(frame, hdu=1)                        # by index (same HDU)
c = euclidrust.extract_cutout(frame, 266.09, 67.32, 50, 50, hdu="DET11.SCI")

(names are matched case-insensitively; if several HDUs share a name, the first one is used).

Compressed files (.fits.fz, as written by fpack) work the same way: the image is in extension 1, so pass hdu=1.

stats = euclidrust.fits_stats("image.fits.fz", hdu=1)
cutout = euclidrust.extract_cutout("image.fits.fz", 150.0, 2.2, 50, 50, hdu=1)   # decodes only the tiles it needs

Lossy (quantized) floating-point images decode to exactly the values astropy/cfitsio would give, which are not the original pre-compression values; that loss belongs to the format, not to euclidrust.

Reference:

Function	Returns	Fields
fits_stats(path, hdu=None)	FitsStats	min ,max ,mean ,naxis1 ,naxis2 ,bitpix ,threads_used ,elapsed_ms , andmin_int ,max_int (exact integers for 64-bit integer images, otherwiseNone ),n_nan (number of NaN pixels, whichmin ,max andmean leave out)
extract_cutout(path, ra_deg, dec_deg, width, height, hdu=None)	Cutout	data (flat list of floats, row-major),data_int (the same pixels as exact Python ints for 64-bit integer images, otherwiseNone ),width ,height ,pixel_x_center ,pixel_y_center ,elapsed_ms
export_cutout_parquet(path, ra_deg, dec_deg, width, height, out_path, hdu=None)	ParquetExport	out_path ,rows ,elapsed_ms
fits_stats_many(paths, hdu=None, skip_errors=False)	list of FitsStats (orNone )	one per file, in the same order; see Many images at once
fits_stats_all_hdus(path, name_contains=None)	list of HduStats	per image HDU: index ,name ,naxis1 ,naxis2 ,bitpix ,min ,max ,mean ,min_int ,max_int ,n_nan ,elapsed_ms
extract_cutouts(path, ra_deg, dec_deg, width, height, hdu=None, skip_errors=False)	list of Cutout (orNone )	one per position, same fields as extract_cutout
list_hdus(path)	list of HduInfo	per HDU: index ,kind (PRIMARY ,IMAGE ,BINTABLE ...),name ,bitpix ,shape (NAXIS1, NAXIS2...),readable ,note (why it is not readable)
detected_cpus()	int	number of logical cores found on this machine

Notes: hdu=None (the default) reads the primary HDU, exactly like before. Cutout.data is a plain Python list, so convert it with np.array(...) as above (fine for cutouts; it is not meant for multi-megapixel regions). The Parquet file has one row per pixel with columns x, y (absolute pixel coordinates in the source image) and value (an exact Int64 / UInt64 column for 64-bit integer images, Float64 otherwise). All errors (file not found, wrong kind of HDU, cutout outside the image, unsupported WCS...) are raised as ValueError.

When you have many images, three functions do the looping in Rust, in parallel on all the cores, with the GIL released:

results = euclidrust.fits_stats_many(paths)
results = euclidrust.fits_stats_many(paths, hdu=1)               # hdu applies to every file
results = euclidrust.fits_stats_many(paths, skip_errors=True)    # None for a file that cannot be read

for r in euclidrust.fits_stats_all_hdus("frame.fits", name_contains="SCI"):
    print(r.index, r.name, r.min, r.max, r.mean)

cutouts = euclidrust.extract_cutouts("tile.fits", ras, decs, 50, 50)
images = [np.array(c.data).reshape(50, 50) for c in cutouts]

Details:

• Results come back in the order of the inputs. fits_stats_all_hdus returns the HDUs in file order, skipping the ones that are not readable images (the empty primary HDU, tables...);name_contains keeps only the HDUs whoseEXTNAME contains that text, case-insensitively (for example"SCI" to leave out theRMS andDQ images).
• By default the first problem raises a ValueError that names the file (or the position, for cutouts). Withskip_errors=True , the failed items areNone instead, so you can process a whole list and look at the failures afterwards.
• extract_cutouts opens the file and reads its header once, then extracts the positions in parallel.
• The GIL is released during all the calls (single or batch), so other Python threads keep running while Rust works. For example, a counter in the main thread reached 731,463 while one 155 msfits_stats call on the 1.47 GB mosaic ran in another thread.

What this buys you, measured on real Euclid data (median of 3 runs, warm cache; the 40-file case uses 40 copies of one detector, so everything is in the OS cache):

Task	astropy	euclidrust : a loop of single calls	euclidrust : one batch call
Stats of the 48 image HDUs of a NISP frame (799 MB)	484 ms	102 ms	72 ms
Stats of the 144 image HDUs of a VIS frame (2.4 GB)	1309 ms	693 ms	299 ms
Stats of 40 files (17 MB each, uncompressed)	274 ms	68 ms	55 ms
Stats of 40 Rice-compressed files	3317 ms	583 ms	572 ms
1000 cutouts 50×50 from the 1.47 GB mosaic	108 ms (WCS built once) / 2835 ms (WCS and file reopened each time)	216 ms	39 ms

Every result is identical to astropy's and to the single calls'. How to read the table: the batch call helps most where each item is small (many cutouts, frames with many small HDUs). On big or compressed files a single call already uses all the cores, so a batch call is barely faster than a loop; the large gain over astropy there comes from reading and decoding in parallel, not from batching.

euclidrust is one binary for every machine. When a computation could use an NVIDIA GPU, it looks for one; if there is none (or the driver is missing, or the card is too old) it uses the CPU and gives the same answers. There is no separate "GPU build" and nothing to install besides the NVIDIA driver you already have to display graphics: the CUDA toolkit is not needed, because the GPU programs are built into the library.

euclidrust.gpu_info()

euclidrust.fits_stats(path)                   # device="auto" is the default: the library chooses
euclidrust.fits_stats(path, device="cpu")     # always the CPU
euclidrust.fits_stats(path, device="gpu")     # always the GPU; a ValueError says why if it cannot
s = euclidrust.fits_stats(path); print(s.device)   # "cpu" or the name of the GPU that computed it

img = euclidrust.GpuImage(path)               # copy an image to the card once...
img.stats(); img.stats()                      # ...then compute on it with no copy at all

The environment variable EUCLIDRUST_DEVICE (cpu, gpu or auto) sets the default for a whole program, for example EUCLIDRUST_DEVICE=cpu to keep everything on the CPU.

Image	auto uses	Why
Rice-compressed ( RICE_1 ), at least 256 tiles and 64 MB	the GPU	Only the compressed data (five times smaller) crosses the PCIe bus, and decoding thousands of tiles in parallel is what a GPU does well: one GPU thread decodes one tile and accumulates its statistics, without ever storing the decompressed image
Raw (uncompressed)	the CPU	The 1.47 GB tile has to cross the PCIe bus and the host memory twice, which costs more than the whole CPU computation
Anything else (gzip, lossless floats, 64-bit integers, small files)	the CPU	Not implemented on the GPU, or too small to repay starting CUDA (about 0.1 s, once per process)

Measured on the Euclid Q1 tile used throughout this README (milliseconds, median, lower is better; the full table, with the machines, is in benchmarks/COMPARISON.md):

	astropy CPU i7	astropy CPU i9	EuclidRust CPU i7 NUC	EuclidRust CPU i9	Euclid GPU RTX 4070 Ti
Raw tile 1.47 GB: min / max / mean	478	1,187	132	119	194
same, image already on the GPU					10.1
50×50 cutout (always on the CPU)	2.70	3.33	0.22	0.23	(CPU)
Rice-compressed tile 0.31 GB: min / max / mean	3,603	3,849	1,211	384	67
50×50 cutout (always on the CPU)	17.1	19.3	3.55	2.51	(CPU)

Three results look odd at first sight. Each one was investigated and measured, and none is a measurement error: the same file (1,474,565,760 bytes), the same Python and astropy versions (numpy 2.5.4, astropy 8.0.2), the data already in the operating system's cache, the median of several runs, and the same minimum, maximum and mean to the last digit on every machine.

1. astropy is slower on the i9 (1,187 ms) than on the i7 (478 ms), although the i9 is the more powerful processor. astropy and numpy compute onone core , so the 10 cores of the i9 do not help them. FITS stores numbers big-endian (>f4 ), and numpy has to swap the bytes of every value before it can compute. I measured numpy'smin() on 64 MB of float32: on ordinary (little-endian) data the two machines are comparable (20.3 GB/s on the PC, 23.8 GB/s on the NUC), but on big-endian data the PC with Windows reaches 6.9 GB/s and the NUC with Linux 15.0 GB/s. So that column mostly measures numpy's big-endian path on Windows, not the processor. Reading the file adds about 0.5 s on the PC. I did not investigate why numpy is slower there (compiler or build options are the likely suspects).
2. EuclidRust gains little on the i9 (119 ms) over the i7 (132 ms) for the raw tile, despite more than twice the cores and a third more clock speed. The work is limited by how fast the memory delivers 1.47 GB, not by the processor. The PC's memory runs at 2933 MT/s (measured at about 37 GB/s with many threads, after an earlier XMP misconfiguration was fixed — it ran at 2133 MT/s, ~28 GB/s, before), the NUC's at 2400 MT/s. Adding cores cannot help a memory-bound computation much. On the Rice-compressed tile, decoding each tile leaves more work per byte fetched, so the i9 pulls further ahead of the i7 there (384 ms vs 1,211 ms).
3. The GPU is slower than the CPU on the raw image (194 ms against 119 ms). The image has to cross the PCIe 3.0 bus (12.9 GB/s measured, so at least 114 ms for 1.47 GB) after a first copy into page-locked memory, whereas the CPU reads it once. The computation on the card takes 10 ms: the time is spent moving the data. The GPU wins where the transfer is small or paid once: the Rice-compressed file (0.31 GB, 67 ms) and an image already in GPU memory (10 ms).

The speed-ups over astropy are real for the user, who runs astropy as it is, but they combine two things: using every core, and not byte-swapping. Compare EuclidRust with the other columns of the same machine to see the first one alone.

What to take from it:

• Compressed images: the GPU is the clear winner , 5.7 times faster than the same code on the 20-thread CPU, and 57 times faster than astropy, withidentical results (minimum, maximum and mean agree to the last digit, the GPU running the same arithmetic as the CPU code).

• Raw images, one computation: the GPU is slower than the CPU here (194 ms against 119 ms). This PC's CPU only offers PCIe 3.0 (12.9 GB/s measured, so 114 ms at best just to move the data) and the file has to be read into page-locked memory first. Once the image is on the card, the computation takes10 ms , 12 times faster than the CPU: that is whereGpuImage pays, when you compute several things on the same image.

• Cutouts, batches and everything else run on the CPU.

• Astropy has no GPU support at all (it is single-threaded C code); a Python user who wants a GPU today has to read the data with astropy and send it to another library.

• Hardware and driver: an NVIDIA GPU of compute capability 7.5 or higher (RTX 20 series and newer) and a driver recent enough for the built-in kernels (PTX version 9.3, from CUDA 13.3). With an older driver, the library recompiles the kernels itself if the CUDA toolkit is installed (its compiler, NVRTC), and otherwise falls back to the CPU.

• Only tested on one GPU: a GeForce RTX 4070 Ti under Windows 11 (driver 610.88). Linux with an NVIDIA card runs the same code but has not been tried. The no-GPU case was tested on Windows (a simulated machine without CUDA) and on the NUC (Debian, no NVIDIA card).

• Not on the GPU: gzip-compressed images, lossless floating-point tiles, 64-bit integers (they need exact integer arithmetic), and everything about WCS and cutouts. Withdevice="gpu" these raise an error that says so; withauto they use the CPU.

• One GPU, one call at a time: device 0 is used (CUDA_VISIBLE_DEVICES chooses which card that is), and GPU calls from several threads wait for each other.

• Memory: the library keeps 128 MiB of page-locked RAM for the transfers once the GPU has been used (two staging buffers of 64 MiB: I measured 16 MiB, 64 MiB and 256 MiB, and 64 is the fastest, because the limit is memory bandwidth, not quantity), plus the image itself in the card's memory forGpuImage .

• First GPU call of a process: starting CUDA costs about 0.1 to 0.2 s. For a script that makes a single call on a small file, that is more than the CPU takes:auto accounts for it with the 64 MB / 256 tiles thresholds above.

To reproduce the table: run python scripts/benchmark_vs_astropy_real.py FILE --gpu --label "i9 (PC)" --workload "raw tile" --json out.json on each machine, then python scripts/make_comparison_table.py *.json > benchmarks/COMPARISON.md. After editing the GPU programs (euclidrust/src/kernels.cu), regenerate the built-in PTX with python scripts/build_ptx.py (this one needs the CUDA toolkit) and rebuild.

tests/test_input_validation.py builds small FITS files with astropy and checks that supported files work (primary HDU, extensions, BITPIX 16 and 64 including exact unsigned 64-bit values and gzip-compressed 64-bit images, BZERO/ BSCALE), that unsupported ones are rejected with a clear error, that every supported projection, matrix form, the TPV distortion and the SIP distortion (including real SIP headers from astropy's own test data) agree with astropy/wcslib, and that tile-compressed images (RICE_1, GZIP_1, GZIP_2, quantization and dithering) decode exactly like astropy/cfitsio, and that the batch functions return what the single calls return. It needs no download. With euclidrust installed and pip install -r scripts/requirements.txt done:

python tests/test_input_validation.py     # prints "all checks passed"

(pytest tests/ works too.)

The script scripts/benchmark_vs_astropy_real.py times euclidrust against astropy on a file of your choice, and checks that both give the same answer.

1. Get a FITS file. Use any file that matches the requirements above. To reproduce the numbers in this README, use the Euclid Q1 tile they were measured on, from the public IRSA mirror (anonymous S3 bucket nasa-irsa-euclid-q1, no account needed, 1.47 GB):

Download it with your browser, curl or wget, and put it anywhere. To try a calibrated frame instead (images in extensions), these are public too:

• NISP, 799 MB (TAN WCS):https://nasa-irsa-euclid-q1.s3.amazonaws.com/q1/NIR/2681/EUC_NIR_W-CAL-IMAGE-BKG_H-2681-0_20240930T175032.222263Z.fits
• VIS, 2.4 GB (TPV WCS):https://nasa-irsa-euclid-q1.s3.amazonaws.com/q1/VIS/2681/EUC_VIS_SWL-BKG-002681-00-1-0000000__20241017T042839.727792Z.fits

2. Install the benchmark's two dependencies (inside the same virtual environment you installed euclidrust in):

pip install -r scripts/requirements.txt     # numpy + astropy

3. Run it, giving the path to the file:

python scripts/benchmark_vs_astropy_real.py /path/to/tile.fits
python scripts/benchmark_vs_astropy_real.py /path/to/frame.fits --hdu DET11.SCI   # image in an extension

Optional flags: --hdu (the HDU holding the image, an EXTNAME or a 0-based index; default: the primary HDU), --ra / --dec (cutout center in degrees; default: the sky position of the image center, taken from the file's own WCS), --size (cutout width and height in pixels, default 50) and --runs (repetitions per measurement, default 5). If the HDU is tile-compressed, astropy has to decompress the whole image for the global stats and reads the cutout through .section; the script handles both. If the file has no sky WCS (raw detector frames, for example), the cutout comparison is skipped with an explanation and the global-stats comparison still runs. If the file uses BZERO/ BSCALE, astropy cannot memory-map it and reads it into RAM; the script says so, and the comparison is then less favourable to astropy than usual for reasons unrelated to euclidrust.

4. Read the output. It prints, for each of the two operations, the median time of both tools, the speedup, and a correctness check:

--- 1. Global stats (min/max/mean over the whole file) ---
astropy (memmap, numpy)      :   476.49 ms
euclidrust (Rust/rayon)      :   130.05 ms
Speedup                      : ×3.66  (72.7% less time)
Result check: astropy min/max/mean = (-111.43283081054688, 15220.939453125, 1.0181564092636108)
              rust    min/max/mean = (-111.43283081054688, 15220.939453125, 1.0181568676826722)

--- 2. 50x50 cutout extraction by RA/Dec (273.071586, 68.000028) ---
astropy (WCS + memmap)       :    2.670 ms
euclidrust (Rust/mmap)       :    0.224 ms
Speedup                      : ×11.91  (91.6% less time)
...
              all 2500 pixels identical: True

(That is a real run on the development machine; your times will differ with your hardware, and a re-run on the same machine moves by a few percent.) min/ max and every cutout pixel should match exactly; mean differs in the 7th significant figure because the two tools sum floats in a different order. Run it twice: the first run may be slower while the operating system loads the file into its cache.

Symptom	Cause and fix
rustc: command not found or pip fails with "Cargo, the Rust package manager, is not installed"	Rust is not installed, or the terminal was opened before installing it. Install via rustup and open anew terminal.
Build error mentioning edition2024 orrust-version	Your Rust is older than 1.88. Run rustup update .
linker 'cc' not found (Linux) /link.exe not found (Windows)	Install a C toolchain: sudo apt install build-essential , or the Visual Studio C++ Build Tools.
ModuleNotFoundError: No module named 'euclidrust'	You are in a different Python than the one you installed into. Activate the virtual environment ( source .venv/bin/activate ) and retry; check withpip show euclidrust .
ValueError: the primary HDU holds no image (NAXIS = 0); it contains: ...	The image is in an extension (Euclid calibrated/raw frames). Pass hdu= with the name or index shown in the message.
ValueError: no HDU named '...' /no HDU with index ...	Wrong hdu value. The message lists the available HDUs;euclidrust.list_hdus(path) shows them all.
ValueError: ... tile-compressed image uses ZCMPTYPE = ... or... is a BINTABLE table, not an image	That HDU uses a compression algorithm that is not implemented (only RICE_1 ,GZIP_1 ,GZIP_2 are), or it is a table.
ValueError: NAXIS = 3 is not supported orno END card found	A 3-D cube, or a file that is not FITS. Only 2-D images are handled.
ValueError: unsupported WCS projection ...	The cutout functions implement the zenithal projections ( TAN ,SIN ,ARC ,STG ,ZEA ,ZPN ),TAN-SIP andTPV . Cylindrical and all-sky projections, and WCS in other coordinate systems (galactic, ecliptic), are not supported yet.
ValueError: this header declares distortion corrections stored in lookup tables ...	The header has HST-style CPDIS /D2IMDIS keywords. Their correction tables are not applied, so cutouts are refused.fits_stats still works.
ValueError: this file has no WCS header	The cutout functions need CRVAL1/2 ,CRPIX1/2 and a linear matrix (CD , orCDELT with an optionalPC ) in the chosen HDU.fits_stats still works without them.
ValueError: cutout outside the image	The RA/Dec you asked for falls outside the image. Check the coordinates, and for frames make sure hdu= is the detector you mean (each detector covers a different patch of sky).
The benchmark exits asking for numpy /astropy	Run pip install -r scripts/requirements.txt in the same environment.

Euclid will deliver petabytes of high-resolution imaging and catalogue data over its six-year survey. The standard astronomy stack for this kind of workload — astropy.io.fits on top of the C library cfitsio — has two structural limitations that become painful at this scale, especially for independent researchers or small labs without access to the official consortium's HPC clusters:

• The GIL. Even on a machine with 8 logical cores, a single Python thread holds the Global Interpreter Lock whilecfitsio does its work, so "read 1000 FITS files in parallel" rarely scales the way you'd hope from thread count alone.
• Unnecessary copies. Reading a single cutout out of a multi-gigabyte tile often pulls far more data into memory than the cutout itself requires.

euclidrust addresses both with a small Rust core: mmap-based zero-copy access to uncompressed FITS data, and a rayon thread pool that auto-sizes itself to the host's logical core count (the GIL is released while it works) — there is no hardcoded thread count, so the same binary scales itself from a laptop to a many-core workstation.

Three single-image operations (plus list_hdus to see what a file contains, and batch versions of them, see Many images at once), each independently benchmarked against astropy on identical input and compared for numerical correctness, not just speed:

Function	What it computes	Notes
fits_stats(path, hdu=None)	Global min /max /mean over one image HDU (the primary one by default)	mmap +rayon , chunked reduction tuned for this workload
extract_cutout(path, ra_deg, dec_deg, width, height, hdu=None)	Pixel-space cutout around a sky position	Zenithal WCS (TAN, SIN, ARC, STG, ZEA, ZPN), with TPV or SIP distortion; FITS 1-indexed convention
export_cutout_parquet(path, ra_deg, dec_deg, width, height, out_path, hdu=None)	Same cutout, written as Parquet	One row per pixel: columns x ,y (absolute pixel coords in the source image),value (physical value,BZERO /BSCALE applied)

All numbers are measured on a real Euclid Quick Data Release 1 (Q1) product from the public IRSA mirror (see Run the benchmark for the exact file):

• File :EUC_MER_BGSUB-MOSAIC-NIR-Y_TILE102160339-...fits — a NISP Y-band, background-subtracted MER mosaic tile near the Euclid Deep Field North. 19200×19200 px,float32 (BITPIX = -32 ), 1.47 GB.
• WCS : standardRA---TAN /DEC--TAN with aCD matrix and no SIP terms.

Measured on the development machine this project was built on: an Intel Core i7-8559U (4 cores / 8 threads, 2.7 GHz base, up to 4.5 GHz), 31 GiB of RAM, an NVMe SSD, no discrete GPU — a laptop-class machine, not a cluster. This matters: the point of these numbers is "what you get on hardware a working astronomer might actually have on their desk," not a best-case HPC result. Each number is the median of 5 runs after a warm-up pass, so the OS page cache is in a comparable state for both tools.

Operation	astropy	euclidrust	Speedup	Time saved
Global stats (min/max/mean, full 1.47 GB)	478.30 ms	130.35 ms	×3.67	−72.7%
Cutout extraction (50×50 px, by RA/Dec)	2.664 ms	0.219 ms	×12.16	−91.8%

Re-runs of the same script on the same machine land within a few percent of these (for example ×3.65 and ×12.36).

min and max matched bit-for-bit between the two implementations. mean matched to 7 significant figures (1.0181564... vs 1.0181568...) — the residual difference is floating-point summation order (NumPy's pairwise summation vs. this project's parallel tree reduction; floating-point addition is not strictly associative), not a computational error. The cutout pixels are identical.

Two things are worth noting about why the gap differs between the two operations:

• The cutout benchmark shows the larger relative gain, because it is dominated by per-call Python/WCS/NumPy object overhead — exactly the overhead this project is designed to remove. The benchmark builds theWCS object inside every call, like a naive script would; an astropy user who builds it once and reuses it will see a smaller gap.
• The global stats benchmark is ultimately bounded by memory bandwidth, not CPU: scanning a full tile moves more than a gigabyte through RAM regardless of which language does the summing, so the ceiling on the achievable speedup is lower and both tools get closer to it.

Tile-compressed images (RICE_1, GZIP_1, GZIP_2, with quantized or lossless floats) are decoded tile by tile, with the tiles spread over the threads; a cutout decodes only the tiles it touches. Results are checked against astropy (cfitsio): min/max and every cutout pixel are identical, and mean agrees to ~10^-9 relative, over 40 test configurations (algorithms, pixel types, tile shapes, dither methods and seeds, NaN pixels, unsigned 16-bit) and on the real Euclid data below.

I had no Euclid-produced compressed file (none of the Q1 images I checked were compressed), so the real-data runs use real Euclid images that I compressed myself with astropy's CompImageHDU (cfitsio underneath), with its defaults (quantization level 16, dithering) unless stated:

Real Euclid data, compressed with astropy	Global stats: astropy ->euclidrust	50×50 cutout: astropy (.section ) ->euclidrust
NISP detector, 2040², RICE_1	78.6 ms -> 14.5 ms (×5.4)	6.96 ms -> 0.54 ms (×12.9)
NISP detector, GZIP_2	219 ms -> 22.4 ms (×9.8)	10.5 ms -> 0.80 ms (×13.2)
NISP detector, GZIP_2 lossless	176 ms -> 12.8 ms (×13.7)	8.58 ms -> 0.48 ms (×17.8)
VIS detector, 2048×2066, RICE_1	72.7 ms -> 11.2 ms (×6.5)	7.37 ms -> 0.46 ms (×15.9)
MER mosaic, 19200², RICE_1 (1.47 GB -> 308 MB)	3719 ms -> 1216 ms (×3.06)	14.96 ms -> 3.71 ms (×4.0)

Median of 3 runs on the development machine, warm cache. For astropy, the global stats include decompressing the whole image (single-threaded), which is what an astropy user does; the gap is wider than for uncompressed files because euclidrust decompresses tiles on all the cores. On the big mosaic the tiles are whole image rows (19200 pixels wide), so a cutout still has to decode about fifty full rows, which is why the cutout gain is smaller there. These are decoding speeds only: the files were produced by cfitsio, and Euclid's own compressed products (if any) may use other settings.

The speedups above are the result of iterating against measurements, not a single clever trick:

1. Naive version — onerayon task per pixel viainto_par_iter() . Correct, but the per-item scheduling overhead dominates.
2. Chunked reduction — split the buffer into 65,536-pixel blocks, each reduced by a tight sequential loop (which LLVM auto-vectorizes atopt-level = 3 without any explicit SIMD intrinsics), withrayon only parallelizing across blocks. This cut the stats time by more than half (×2.24 on a ~1 GB file during development).
3. mmap vs. eager read — reading the whole file into aVec<u8> first was measured to be1.87× slower thanmmap , confirming that letting the OS page cache serve data on demand beats an eager bulk copy for this access pattern.
4. Stopping point — once (2) and (3) were in place, further micro-optimization (explicit SIMD intrinsics, finer chunk tuning) showed diminishing returns against the memory-bandwidth ceiling identified above. That's the deliberate stopping condition for this round of tuning, not an oversight.

Build settings: cargo release profile with lto = "fat", codegen-units = 1, opt-level = 3. Thread pool size is read at runtime from the host (rayon::current_num_threads() / num_cpus::get()) — never hardcoded, so the same wheel adapts to whatever machine runs it.

• 2-D images only. No tables or catalogues, no 3-D cubes, no random-groups files. Each call reads one HDU of one file (the batch functions loop over files, HDUs or positions for you).

• Compression: RICE_1, GZIP_1 and GZIP_2 only.HCOMPRESS_1 andPLIO_1 are not implemented and are refused with a clear error. The compressed files tested were produced by astropy/cfitsio, not by Euclid's pipeline.

• WCS coverage. Zenithal projections (TAN ,SIN ,ARC ,STG ,ZEA ,ZPN ), with TPV or SIP distortion. Not supported: cylindrical and all-sky projections (CAR ,AIT ,MOL ,HPX ...), coordinate systems other than RA/Dec, alternate WCS (A -Z : only the primary one is read), and distortion given by lookup tables (HSTCPDIS /D2IMDIS , which makes cutouts fail with an error). The sky frame is not read:RA andDec are taken in the file's own frame, andRADESYS /EQUINOX are neither read nor converted.

• Null pixels. NaN pixels are left out of the statistics; the integerBLANK null value is ignored (see above).

• A cutout is a crop, not a resampling. It is a rectangle of whole pixels whose corner isround(centre - size / 2) (halves round away from zero, whereas Python'sround rounds halves to even, so a hand-made comparison can differ by one pixel on an exact tie). The requested position can therefore be up to half a pixel from the middle of the crop; the exact fractional position is returned inpixel_x_center /pixel_y_center . There is no interpolation, rotation or sub-pixel shift, and sizes are in pixels, not in angle.

• Cutouts must lie fully inside the image. One that touches the edge or extends beyond it is refused (cutout outside the image ); there is no padding and no partial cutout.

• Cutout.data is a plain Python list of floats (anddata_int a list of ints): fine for thumbnails, not meant for multi-megapixel regions.

• Statistics.mean is the mean of the valid (non-NaN) pixels. It agrees with numpy to about 7 significant figures on float32 images, not to the last digit, because the sums are made in a different order. Images with NaN pixels take about three times longer than clean ones, because a second, NaN-aware pass is needed (still faster than numpy'snanmin +nanmax +nanmean ; measured: 70 ms against 164 ms on 256 MB with 7% NaN).

• 64-bit integers are exact only in the integer fields. For 64-bit integer images (signed, or unsigned viaBZERO = 2^63 ; raw or gzip-compressed) the exact values aremin_int /max_int anddata_int , and the Parquetvalue column. The float fields (min ,max ,mean ,data ) arefloat64 , exact only up to 2^53 (about 9 × 10^15), which is inherent to the type. The exact path needs whole-numberBSCALE (|BSCALE| ≤ 2^20) andBZERO (|BZERO| ≤ 2^64); with a fractionalBSCALE orBZERO the physical values are not integers anyway, so everything isfloat64 .

• The sky-to-pixel inverse of TPV and SIP is iterative. It was checked against wcslib (through astropy) on over 120 synthetic headers and on real ones: the real Euclid VIS frame (TPV, 3 × 10^-8 pixel) and the real SIP headers shipped with astropy's tests, including HST/ACS (4 × 10^-9 pixel at worst). It can pick the wrong root when a SIP distortion is so strong (hundreds of pixels) that the inverse is not unique, which does not happen with real instruments.

• Real-data validation is small. The benchmark numbers above come from a single Q1 NISP mosaic tile out of the ~35 TB Q1 release. Correctness (not speed) was also checked against astropy on one NISP frame, one VIS frame, one raw NISP detector and one segmentation map. That is enough to show the approach works on real archival products, not a representative sample across instruments, bands and tiles.

• Tested on two machines. (1) An Intel NUC running Debian 13 (Linux 6.12, x86_64): Intel Core i7-8559U (4 cores / 8 threads, 2.7 GHz, up to 4.5 GHz), 31 GiB of RAM, an NVMe SSD, Python 3.13.5, Rust 1.100, no NVIDIA GPU. (2) A desktop PC running Windows 11: Intel Core i9-10850K (10 cores / 20 threads), 64 GiB of RAM, an NVIDIA GeForce RTX 4070 Ti (12 GB), Python 3.12.10, Rust 1.95. macOS and ARM have not been tried, nor Linux with an NVIDIA GPU, nor machines with much less memory (the images are memory-mapped, so a laptop with 8 GB of RAM should cope with the 1.47 GB tile, but I have not measured it).

• No continuous integration, no pre-built packages. The tests are run by hand (there is no CI workflow). The package is not on PyPI: you build it yourself, which needs a Rust toolchain. The interface is version 0.1.0 and may change.

• The files are memory-mapped and must not change while they are read. If another process truncates or overwrites a file during a call, the operating system can kill the Python process (this is howmmap works; I have not tried to provoke it). Network or FUSE file systems are untested.

• On Windows, mapping a whole big file is expensive to undo. Closing a view of the 1.47 GB tile took 165 ms, more than the computation. The library therefore reads big images through many small mapped windows, one per thread, which spreads that cost (255 ms became 123 ms). Linux does not have this problem. Astropy'smemmap is also much slower on Windows (1.3 s against 0.48 s on the NUC for the same statistics).

• All the cores are used by default. The number of threads follows the machine; to limit it (a shared node, for example) set the environment variableRAYON_NUM_THREADS before starting Python, e.g.RAYON_NUM_THREADS=4 . The GIL is released during the calls.

• Memory with compressed images. Each worker thread decodes whole tiles intofloat64 (8 bytes per pixel), so an image compressed as a single huge tile needs about eight times its pixel count in bytes. There is no tile cache: many cutouts falling in the same tile decode it again each time.

• Batch results are held in memory.fits_stats_many andextract_cutouts return everything as Python lists at the end: there is no streaming and no progress callback, so a huge list of cutouts uses memory proportional to its size.hdu= is the same for every file of afits_stats_many call.

This is the first of four planned components for Euclid-scale data processing, chosen first specifically because it needed no GPU/cluster and is the foundation the others build on:

1. ✅ FITS streaming / cutout / Parquet export (this package). Next steps for this component, if someone needs them: the remaining compression algorithms (HCOMPRESS_1 ,PLIO_1 ) and the cylindrical / all-sky projections.
2. Catalogue cross-matching at the hundred-million-to-billion-object scale (HEALPix / spherical KD-tree indexing, memory-mapped for datasets larger than available RAM).
3. Morphology / weak-lensing profile fitting (PyO3-bound, parallel) — feasible at prototype scale on commodity hardware, not yet at billion-galaxy scale without a cluster.
4. Interactive WASM/WebGPU visualization of large mosaics and catalogues.

MIT licensed (see LICENSE) — free to use, copy, modify, merge, publish, distribute and sublicense, with no restriction beyond keeping the copyright notice. Experimental, actively developed.

source & further reading

github.com — original article Training and Inference on Huawei Ascend NPUs Cowork – Self-hosted task queue and dashboard for Claude, Codex, and AI agents Show HN: DeepSeek V3 in ~200 lines of pure Jax