1. Introduction
1.1. The square kilometre array
Radio Astronomy is undergoing a paradigm shift with the planning for a number of next-generation instruments, such as the Square Kilometre Array (SKA, Dewdney et al. Reference Dewdney, Hall, Schilizzi and Lazio2009), the next-generation Very Large Array (ngVLA, Selina et al. Reference Selina and Murphy2018) and the next-generation Event Horizon Telescope (ngEHT, Doeleman et al. Reference Doeleman2023). All of these provide an order of magnitude increase in bandwidth and a few orders of magnitude in collecting area (and sensitivity) relative to current radio telescopes. This enhancement will provide us with the opportunities to survey the radio sky in exquisite detail, detect the signal from the epoch of reionisation when the first stars were born and measure the spectral signal from millions of galaxies. To grasp this potential massive improvement in our understanding, the radio astronomy community must manage and process unprecedented volumes of data. We are focusing on the SKA, as Australia is a founding member of the collaboration. The SKA will be built in Australia for frequencies spanning 50 to 350 MHz (SKA-Low, Labate et al. Reference Labate2022) and in South Africa for frequencies from 350 MHz to 15 GHz (SKA-Mid, Swart et al. Reference Swart, Dewdney and Cremonini2022). Phase-1 will have about 500 40 m aperture array elements in SKA-Low and about 200 15 m parabolic dishes in SKA-Mid. Construction has commenced and preliminary observations will be made from 2025. The data rates out of the correlators will be about 1 TB/s, which will need to be captured in a local buffer and then processed on the day - as storage will be limited (see discussions in McMullin et al. Reference McMullin, Marshall, Spyromilio and Usuda2020).
1.2. Radio astronomy data
These next-generation instruments are all interferometers, as this is the most cost effective way to build a high resolution instrument. Multiple (significantly cheaper) small dishes can be distributed spatially, to provide the high spatial frequencies in the image. Interferometry consists of the measurements of the coherence function, at a spatial separation. In astronomy these correspond to the correlation of two signals, taken by antennas separated by a ‘baseline’ vector and, as the signal is travelling from astronomical distances, represents the Fourier term of the transform of the brightness distribution of the observable sky, at the spatial frequency given by the baseline (e.g. Van Citter 1934, discussed in Thompson et al. Reference Thompson, Moran and Swenson2017). These measurements of the spatial coherence over the interferometric array are the visibilities, which are the output of the correlator that connects all of the stations in the array.
Radio data presents a unique challenge. Much of the data is noise as the signal from any spatially concentrated region of high radio emission is spread between all of the visibility samples. Furthermore, the radio sky is largely empty and flux from the sources (expressed as T
$_\mathrm{sky}$
, the measured intensity in Kelvin) is much weaker than the visibility signals, which are dominated by the receiver temperature, T
$_\mathrm{rec}$
.
This data analysis challenge is combined with a data volume challenge. Radio astronomy data volumes from current generation telescopes are of PB-scale. This data is also often stored as a Measurement-Set (Kemball & Wieringa Reference Kemball and Wieringa2000), a format in which visibility and single-dish data are stored to accommodate synthesis. Although this format has been historically very useful, it does not scale particularly well and often the science process requires non-optimal access, giving rise to additional I/O load. This data challenge will only increase once next-generation telescopes become operational. The SKA strategy for multi-epoch (i.e. ‘deep’) projects is to combine fully processed ‘daily’ images. The data-reduction for image formation allows these to be stored for every observing session and then co-added (or image-stacked) to achieve the final sensitivity. The downside of this approach is that the deconvolution in the imaging, which is to correct for the imperfect instrumental sampling, is only partially applied and the final product is a combination of ‘cleaned’ and ‘uncleaned’ emission. The traditional processing operates simultaneously with the individual datasets co-added to a common grid, which is imaged, fully deconvolved at the final sensitivity and model subtracted in major cycles. This obviously both requires all the dataset to be available simultaneously, which is well beyond scope for SKA and for the Australian SKA Pathfinder (ASKAP). We are investigating an alternative strategy that allows for full deconvolution yet does not require the retention of all the raw datasets. In this case we store a ‘daily’ grid which is the residual visibility grid formed in the final major cycle of daily imaging, along with the response function and extracted model for that observation. This approach we call grid stacking, as we will co-add the uv-grids and do a final (minor) cycle of deconvolution at the final sensitivity. In principle, the uv-grids are not significantly larger than the image; there are the same number of voxels, with the former formed from complex numbers and the later from real numbers. In the companion paper (Rhee et al. Reference Rhee2026) we present a detailed comparative cube analysis between these methods. However, one also needs to store the Point Spread Function (PSF) and the Pre-Conditioning Function (PCF), which are also stored as complex numbers. The PSF grid represents the size, location and weighting of the complex-valued convolutional kernels applied during the gridding process. Both are used in the final imaging to apply weighting to individual visibility cells. The PCF grid represents the baseline sampling of the system, as the Wiener filtering requires the count of the samples per cell without kernel weighting. In this case the imaginary axis is used as an independent storage location for the average kernel size for that cell. The sum of these products can approach the same data volume as the original dataset. Thus, for uv-gridding to be useful we need to compress these three grids; fortuitously these grids are largely made up of zeros, as the baseline sampling in a radio interferometer is very sparse and they compress very efficiently.
An alternative solution to the data volume problem is to apply baseline dependent averaging (BDA) to the data before processing (e.g. Wijnholds et al. Reference Wijnholds, Willis and Salvini2018). This is effective because short baselines have longer (slower) response times to the changes in the sky structure (and the atmosphere), because of their physical length. However, this means the data is averaged before correcting for any non-linear processes in Radio Interferometric imaging, for example some of the kernel corrections. Our approach is to apply these correction kernels to the raw data and then average to the uv-grid (which is the ultimate limit of BDA) and store these uv-grids in compressed format rather than the original visibility data. This is the approach we have implemented in the software for ASKAP and here we report on the impact of both the lossy and lossless compression options.
1.3. The Australian SKA pathfinder
The ASKAP (Johnston et al. Reference Johnston2007; Hotan et al. Reference Hotan2021) radio telescope is one of the SKA precursors and is opening up a new window for large extragalactic
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
surveys beyond the local Universe due to its wide spectral bandwidth and large instantaneous field-of-view (FoV). ASKAP consists of 36 dishes, each of diameter 12 m and equipped with phased array feeds (PAFs) forming multiple (up to 36) receiving beams electronically (DeBoer et al. Reference DeBoer2009; Hampson et al. Reference Hampson2012). The baseline lengths of the full array are from 22 m to 6.3 km. The phased array feed technology allows ASKAP to have a large
$30\,\mathrm{deg}^{2}$
FoV (Bunton & Hay Reference Bunton and Hay2010) and a wide bandwidth of
$288\,\mathrm{MHz}$
with a channel resolution of 18.5-0.58 kHz in the observing frequency between 0.7 and
$1.8\,\mathrm{GHz}$
, which makes it an optimal survey instrument, enabling it to conduct both wide and deep surveys in a comparatively short period of time (Hotan et al. Reference Hotan2021).
1.4. Deep Investigation of Neutral Gas Origins (DINGO)
Deep Investigation of Neutral Gas Origins (DINGO) (Meyer Reference Meyer2009; Rhee et al. Reference Rhee2023) is an ASKAP deep
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
survey project aiming to provide a cosmologically representative dataset for
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
emission, enabling studies of the
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
gas content of galaxies over the past 4 billion years. The sky coverage of the DINGO survey is wider than deep
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
surveys previously conducted and the ongoing deep
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
surveys being carried out with other telescopes such as the JVLAFootnote
a
, MeerKATFootnote
b
and FASTFootnote
c
(CHILESFootnote
d
, LADUMAFootnote
e
/MIGHTEE-HIFootnote
f
, and FUDSFootnote
g
respectively). Due to its large volume coverage (60 deg
$^2$
compared to 20 deg
${^2}$
for MIGHTEE (Maddox et al. Reference Maddox2021) with a similar sensitivity and redshift range), the DINGO survey will reduce cosmic variance on
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
measurements, thereby providing a unique legacy
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
dataset.
DINGO observations are focused on the Galaxy and Mass Assembly (GAMA) (Driver et al. Reference Driver2022) 23 h (G23) field, centred at
$\alpha, \delta$
(J2000) =
$22^\mathrm{h}59^\mathrm{m} 00^\mathrm{s}00, -32^{\circ}18^{\prime}00^{\prime\prime} 0$
. These observations used the full array of ASKAP’s 36 antennas with the 288 MHz bandwidth (15,552 channels) in the observing frequency ranges of
$859.5-1147.5\,\mathrm{ MHz}$
(band 1) and
$1151.5-1439.5\,\mathrm{MHz}$
(band 2). The channel width is
$18.5\,\mathrm{kHz}$
, equivalent to a velocity resolution of
$\sim 3.9(1+z)\,\mathrm{km\,s^{-1}}$
in cosmologically nearby galaxies. The DINGO full survey allocation is 3,200 hr, split into two 1,600 hr allocations for the two frequency ranges of band 1 (
$z \sim$
0.2 to 0.4) and band 2 (z of 0 to 0.08) and alternating between two tiles, each with two dithering-footprints, A and B, stepped in a half beam-width, to ensure more uniform sky sensitivity. We would expect the final data volume will be from 100 days per tile, each of 8 hours and with 36 beams and 2 bands. As each MS file is 57 GB/hour in size the storage requirement for the full dataset will be 6.6 PB, which is beyond feasibility.
The default ASKAP and SKA strategy to address this is to image each day and stack these reduced data products to reach the results. Numerous investigations have suggested that this approach would be sub-optimal (for example, Dodson et al. Reference Dodson2013, Reference Dodson2016; Rozgonyi Reference Rozgonyi2021). This has now demonstrated in Rhee et al. Reference Rhee2026), which is the companion data-analysis paper for this data-methods discussion, where the results from grid-stacking and image-stacking are compared to the traditional approach. What Rhee et al. Reference Rhee2026) shows is that neither of the stacking approaches fully recover the results from traditional imaging, because the deconvolution in the imaging process is non-linear. This leads to a 1% shortfall in the recovered flux in the grid-stacking compared to traditional, but a 8% shortfall for image-stacking. The latter also suffers from significant image quality degradation (see Figures 5 and 3 in Rhee et al. Reference Rhee2026), respectively). The reason for this is that no further deconvolution is possible in image-stacking after the daily deconvolution, whereas in grid-stacking further minor-cycle deconvolution is possible but major cycles is not. This means that the weaker sources in the image-stacked cubes, in particular the low-column density HI, is a mixture of deconvolved and non-deconvolved emission, whereas this problem is much reduced in grid-stacking (and traditional) processing.
There are further approaches that could be developed further. One could perform a further minor cycle on the image-stacked residual cubes (as discussed in Rhee et al. Reference Rhee2026), but one would now need to store the daily residual image, the models and the image-domain PSF (at four times the image size). These would now be larger than the original data-products and would not be highly compressible as there are no longer empty cells. The stacked uv-grids could be extended to support major cycles, by extending the dimensionality, either by separating the data by w-plane or by baseline. As an uncompressed product these would be much larger than the original data, but as the number of occupied cells would not increase significantly (Dodson et al. Reference Dodson2013) they should compress to more or less the same storage footprint. Neither of these ideas have been explored in this work nor are currently supported by ASKAPsoft.
So far a total of
$\sim$
200 hr observations have been conducted in the higher frequency band, between 2019–2024 and with both footprints for one tile. These are formed of 25 individual
$\sim$
8 hr long observations, consisting of 36 beams and two footprints, with 15,552 channels spanning 288 MHz, centred on 1295.5 MHz, of which only the upper half (above the RFI from the GPS satellites) was processed. See Table 1 in Rhee et al. Reference Rhee2026) for full details. We are using these early observations to test the compression options for the proposed alternative strategy of uv-grid stacking.
2. Methods
2.1. ADaptive I/O system version 2 (ADIOS2)
The Adaptable Input Output System version 2 (Godoy et al. Reference Godoy2020), is a software framework with a simple input/output abstraction and a self-describing data model centred around distributed data arrays, allowing multiple applications to publish and subscribe data at large levels of concurrency. It also introduces a larger organising concept, the ‘step’, for driving data production and consumption within applications. ADIOS2 is primarily focused on high-performance and parallel I/O, with its parallel storage performance, the file format, the memory management, and data aggregation algorithms being designed together to be highly scalable in every axis (many processes, many variables, large amounts of data, many output steps). Our I/O backend is implemented using the file-based storage I/O engine in the ADIOS2 library. The engine designates separate MPI ranks as aggregators, which collect data in a streaming fashion, buffering them in large chunks before writing to a disc or a solid-state drive to avoid small I/O requests. As such the data can be independently written by each MPI imaging process, which in the DINGO analysis is many threads that each manage the data reduction for a small range of channels (18 in the analysis presented here). Thus ADIOS2 provides a native parallel I/O library for high performance high speed streaming data transport and the ADIOS image class distributes data between multiple table cells, allowing us to fully utilise the high-performance compute platform. The self-describing data are stored in binary-packed (.bp) format to allow for rapid metadata extraction. Of particular relevance to the astronomical community is that it can be embedded in the CASACORE libraries as an alternative storage manager (if enabled at compile time), and is able to write tables for both images and visibilities. We note that the ADIOS Image Class is not, strictly, in CASACORE, but is ‘research code’ currently available via the CASAREST library. A pre-compiled DINGO-ASKAPSoft image is available on request and we expect this functionality to be eventually available in the default ASKAPSoft release. For the DINGO processing the raw input data column (typically ‘DATA’) could be in an ADIOS format, but it is not in this pipeline; however the image cube is, and written out as columns ‘map’ and ‘mask’. ADIOS2 recently developed a new mechanism to allow applications to use state-of-the-art lossless and lossy compression algorithms. This mechanism makes use of tight integration between I/O and reduction and allows applications to take full advantage of the self-describing formatting and lossy compression techniques.
2.2. MultiGrid adaptive reduction of data (MGARD)
MGARD (Gong et al. Reference Gong2023) offers error-controlled lossy compression rooted in multi-grid theories. It transforms floating-point scientific data into a multilevel representation, followed by quantisation, lossless encoding, and ultimately generating a self-describing compressed buffer. One of MGARD’s notable features is its array of error control options, including
$L^\infty$
,
$L^2$
, point-wise relative
$L^\infty$
, and options to define varied error bounds across regions or different frequency components. This flexibility is valuable for preserving Quantities-of-Interest (QoI) (Gong et al. Reference Gong2022) derived from the reconstructed data. For region-adaptive compression, MGARD accommodates Regions-of-Interest (RoI) specified through either bounding boxes or masks, with the latter especially useful for irregular shaped RoIs. In cases where RoI information is not provided, MGARD employs internal functions to identify regions rich in detail, leveraging data turbulence measured across multiple scales.
2.3. ASKAPsoft
ASKAPsoft is a package that contains the software necessary for processing data from the ASKAP telescope. Its primary purpose is for the full-scale processing of ASKAP data, from the observed visibility data to spectral-line and continuum images. ASKAPsoft provides the well-established routines to image a spectral line dataset. That is to: read the data, apply weighting kernels to set the image parameters (Field of View, sensitivity to low surface brightness or compact objects, etc), resample the data onto a regular sampled grid for inversion using the Fast Fourier Transform, iteratively deconvolve for the limited sampling of Fourier terms and finally, for deep images, stacking of multiple epochs of observing for the final image. The DINGO pipeline reorganises these tasks so that grids are preserved for stacking, rather than the images. This reduces the number of inversions required and improves the quality of the deconvolution.
Imaging radio interferometric data at SKA-scales is expected to be I/O bound due to the massive size of the datasets. The computational costs are dominated by the gridding and the inversion steps, and these two are expected to have approximately similar requirements. Thus any reduction in the size of the datasets would have a significant impact on the total processing time and any reduction in the gridding or inversion would have a significant (albeit smaller) impact on the compute costs.
ASKAPsoft contains a wide variety of programmes and scripts that are useful in the analysis, manipulation, and processing of radio astronomy data. We will primarily use two applications present in ASKAPsoft: imager & cdeconvolver. imager creates spectral-line image cubes, which use frequency as an analogue for distance, allowing for a three-dimensional view of the sky. Thus this single programme includes the reading, weighting, gridding, deconvolution and inversion steps. For the investigations here we used identical imaging parameters parameters as in Rhee et al. Reference Rhee2026) for imager: 30
$^{\prime\prime}$
resolution and 6
$^{\prime\prime}$
cell size, 3 major cycles each terminating when the 5
$\sigma$
noise-level is reached, except for regions where flux has already been cleaned, in which case the threshold is 0.5
$\sigma$
. In our case, imager is used to produce visibility grids from each schedule block, as an intermediate product, which we then combine to produce the final image. That is, the final inversion step is not performed and the normally intermediate data products (i.e. the residual grids and the model components subtracted) are saved for later processing. Due to the sparse nature of these grids, compression is particularly efficient. Furthermore, in principle, at this point the image processing parameters can be changed as needed to adjust weights and even field of view. The deep imaging mentioned in section 1.4 requires stacking these grids over 3,200 hours of observed data to improve sensitivity to sources within the final image. For this task we use cdeconvolver, a bespoke application designed specifically to perform this stacking and complete the imaging process. Further information on these applications can be found at https://www.atnf.csiro.au/computing/software/askapsoft/sdp/docs/current/index.html
2.4. DINGO pipeline
The DINGO workflow is designed using the DALiuGE framework (Wu et al. Reference Wu2017), which is a scientific workflow development and execution framework to process large astronomical data sets at a scale required by the Square Kilometre Array, with a particular focus on efficient data I/O management. The executed workflow is represented by a Directed Acyclic Graph (DAG), where the vertices are an alternating sequence of stateless computational tasks and data nodes, edges connect the output of one task with one input of a data node which in turn connects to the input of the subsequent task. The graphs can be developed using the graphical editor EAGLEFootnote h , allowing the intuitive representation of complex data reduction workflows, consisting of both datasets and algorithmic components. EAGLE also allows for the run-time implementation and monitoring of the execution of the pipelines on distributed resources. By mapping the logical view of a pipeline to its physical realisation, DALiuGE separates the concerns of multiple stakeholders, allowing them to collectively optimise large-scale data processing solutions in a coherent manner. Figure 1 shows the pipeline as constructed in the EAGLE GUI for DALiuGE, where we built the pipeline from components autoloaded from the codebase (using, for example, Doxygen/Sphinx formatted component documentation for the derivation of the input and output parameters).
The graphical representation in EAGLE of the DINGO gridding pipeline. Specific parameters are provided, split among the scattered instances, merged with the global template and the individual ASKAPsoft imager applications are run to generate the schedule block residual visibility grids and the model components subtracted. Following this, in another similar graph, the gridded visibilities are combined, imaged and the model components are added back to produce the final stacked image.
We start with a MeasurementSet data format that contains the calibrated, continuum-subtracted visibilities (the continuum here refers to the components of the data that are not limited in frequency span). This is the input to imager, described above, which produces a residual visibility grid, a PSF grid, a PCF grid and the spectral model components subtracted to produce the residual visibilities. The visibility grid is a grid representation of the 3-dimensional visibilities projected onto a 2-dimensional grid for each frequency channel (producing a 3-dimensional grid). These grids are passed to the cdeconvolver application which, if more than one observation is provided, will sum the grids together as they are read in. These summed grids are then imaged with minor cycles for analysis.
2.4.1. Performance measurements
In Rhee et al. Reference Rhee2026, Table 3) we have made detailed comparisons of the scientific outcomes of grid-stacking, in comparison to both traditional processing of all the daily epochs simultaneously, that is, without consideration of storage limits, and the alternative storage-constrained approach of image-stacking. The focus of this paper is the performance of ADIOS2 and lossless and lossy gridding, which obviously only is relevant for grid-stacking. The use or otherwise of ADIOS has minimal impact on the required Memory, compute and read I/O, which are dominated by the management of the raw visibilities. For the processing of a daily image-cube with grid stacking we found that processing 7776 channels required 80min on 217 cores, with a maximum memory footprint of 30GB, with a required read I/O of 40GB. The write I/O without using ADIOS was about 30GB. The fractional impact of ADIOS on the latter figure is investigated in Section 3. In addition to the spectral line imaging step measured here, the grid-stacking requires the models to be restored, and all methods require subsequent continuum subtraction in the image domain to remove residual continuum emission, and two stages of linear mosaicking, but these are not significant additional computational cost.
2.5. Compression
The focus of the compression here is that of the grids. Due to the limited range of the values in the PCF grid (i.e. limited to half integers), lossless compression provided a compression ratio of
$\sim$
100, with lossy compression providing ratios similar to that of the visibility and PSF grids. For this reason, we are only comparing the compression of visibility and PSF grids. We use lossless and lossy compression with error bounds of
$10^{-3}$
,
$10^{-4}$
,
$10^{-5}$
, and
$10^{-6}$
both as relative and absolute error bounds, although these metrics are somewhat overlapping. The cdeconvolver application failed to complete the imaging with a relative error of
$10^{-3}$
, due to a failure to converge. We believe that this is due to the errors accumulating in the PSF reconstruction which, as it is applied iteratively, is most sensitive to numerical error. To investigate this further we would have used separate compression scales for the PSF as opposed to the visibilities, but this was out of scope of our investigations. The relative error bound is derived against the range of the data being compressed (i.e. against half of the minimum to maximum data values) whereas the absolute error bound is applied without scaling. The typical range of values for the input datafiles are about
$\pm$
100, so the absolute error bounds are two orders of magnitude less than the relative values. For the application to gridded data we quickly discovered that it was vital to ensure that zero values were not changed to values near zero but within the error bound, which is the default behaviour of multi-grid methods. This was solved by compressing both the raw data and a mask which records the position of the zeros. This is vital because the sum of visibilities in a grid cell is normalised against the weight of samples included in that sum, as stored in the PSF/PCF grid. Small near-zero values causes the imaging to diverge. The lossless compression uses the zstd algorithm to compress the grids, although bzip2 also provides a similar level of compression. This provided a consistent compression ratio and (as expected) did not alter the decompressed data in any way.
2.6. Parallel I/O
The integration of ADIOS2 into the I/O stage of ASKAPsoft leveraged the use of CASACORE and the implementation of the existing ADIOS2 storage manager (Wang et al. Reference Wang, Harris and Wicenec2016). An ADIOS2 Image module was developed to bridge the gap between imager and cdeconvolver’s use of image inputs and outputs, and the storage manager’s table interface. These applications use MPI to implement parallel processing which, when passing the communicators to the storage manager (and by extension ADIOS2), enables the I/O to occur in parallel. These capacities have been or are in the process of being merged back into the ASKAPsoft package.
3. Results
3.1. Compression comparison
Figure 2 shows the compression ratios of the visibility and PSF grids as the error bound are increased from zero (lossless) to
$10^{-3}$
for both absolute and relative error measurements. Lossless compression is consistent at a value of 7.5, whereas the lossy compression ratio increases from that limit with increasing error bound, and the compression ratios of the real and imaginary parts of the visibility grid are consistent. The absolute error bounds for the visibilities are effectively more than an order finer than the numerically equivalent relative error bound, and this can be seen in higher compression for the relative results.
The compression ratio for increasing specified error bounds. Shown are the visibility and PSF grid compression ratios for lossless compression, lossy compression with absolute error bounds, and that with relative error bounds.
The absolute residuals between the images produced after lossy compression and the uncompressed equivalent. The colour scales are not normalised as the range of values for each image is significantly different. Ignoring the image corners, the smallest error bounds on the right-hand side lead to maximum residuals of
$\sim10^{-5}$
.
Of the eight tests of the compression, seven completed the imaging and produced, qualitatively, decent images. The reason for the failure of relative-
$10^{-3}$
result is likely caused by the numerical degradation of the PSF, and thus errors introduced in the convolution, and a divergence in the image deconvolution.
The images produced from the compressed data are visually close to identical; it is only in the image differences that the impact of compression can be seen. Figure 3, plotting the difference between the images made with lossy compression and those without, shows that majority of the value deviation is at the corners of the image. Correcting for primary beam weighting across the FoV would exclude these points. We suspect this problem has a common origin with the complete failure of the image-convergence at relative error bounds greater that
$10^{-3}$
.
For the more aggressive compression, that is, the
$10^{-4}$
relative and
$10^{-3}$
absolute error bounds, increased residuals near the centre of the image appear, which would have an impact on the science outcomes. At higher precision error bounds the errors introduced across image are uniform.
An alternate distortion measure is the 2-point correlation of the residuals. This can be produced by performing a Discrete Fourier Transform (DFT) on the residuals, coherently binning the pixels in radius and calculating the square of the absolute value of these pixels. This is mainly to detect ripples in the residuals from poor deconvolution and/or calibration errors. The result is shown in Figure 4 and describes the prominence of patterns of certain scales within the residuals. The higher error bounds all show a sharp cutoff at
$\ell$
=289, which corresponds to the largest radial distance of a non-zero cell on the visibility grid, that is equivalent to the resolution – 30ʹʹ, whereas the residuals for highly compressed data with
$10^{-3}$
absolute error and
$10^{-4}$
relative error show a negative logarithmic trend between
$\ell$
289 and 410. An excess of values above this spatial frequency shows a leakage of values within the PSF grid, adding small scale residuals below the smallest resolvable scale, and is found only for PSFs which have been over-compressed.
The unnormalised 2-point correlation of the residuals of the images produced with varying error bounds. The x-axis represents the spatial frequency bins of the image, where
$\ell$
=0 is equivalent to structures over the whole field of view and
$\ell$
=700 is equivalent to structures of 12ʹʹ, double the pixel size. The colours are as for Figure 3 and same first two error bounds show power at angular scales greater than the resolution, indicating reconstruction errors.
3.2 Parallel I/O comparison
Figure 5 shows the comparison between the standard CASACORE tile manager to write the output grids, the FitsParallel I/O libraries as used for writing large files in ASKAPsoft, and using MPI-based parallel-enabled ADIOS2 tile manager in CASACORE for writing. These are selected using the arguments in the imager parset file and we have taken the average of five runs to reduce timing noise on the multi-user system. In comparison, the single threaded CASA storage format could not compete with these parallel writing modes. However, we note that the most recent version of ASKAPsoft has improved the performance of this mode, but we have not yet assessed this option.
The real runtime (normalised by No. of channels per core) of the imager application CASA here refers to the standard image I/O manager and ADIOS2 refers to the ADIOS2 I/O manager run in parallel, both under CASACORE. Fits refers to the parallel fits writer natively coded in ASKAPsoft. ADIOS2 and FITS are both run in parallel via MPI whereas the CASA Tiled-Storage Manager is run in serial (due to the limitations of this storage manager). The completion time in the top panel is the real time to completion whereas the completion time in the bottom panel is normalised by the number of cores per channel. The CASA data point at 7776 channels is a lower limit (due to the processing timing out at 5 hours).
Altering the input format of the raw data to lossy compressed format could also have a significant impact on the processing speed (Dodson et al. Reference Dodson2025) but this was not investigated, as it was out of scope.
3.3. Comparison of compression on processing time
Figure 6 shows this comparison for the 100-channel case where the different storage managers provide comparable results, but with MGARD compression enabled during the writing stage of the processing. Lossless compression appears to perform significantly better than the same processing with lossy compression. In the current implementation, MGARD is run completely on CPUs, whereas the best performance comes from using the GPUs for compression in parallel. We plan to implement and test this compression using GPUs in the future, and expect that this will be more time-efficient. However, if the processing is I/O bound the computational costs will not be highly significant.
The fractional processing time difference between using the default CASA table manager and parallel threaded ADIOS2 with different error bounds for both absolute and relative error for a 100-channel dataset, where the impact from multithreaded I/O is minimal. The median of 5 runs (points) are shown with solid lines. Along the x-axis we plot the error bound for different compression levels, implemented using ADIOS2 and MGARD to compress the data during the writing stage of the processing. The fractional difference in processing time from that of the CASA manager (270 sec) is plotted on the y-axis. For the lossless compression (error bound of zero) the impact is negative, meaning that the ADIOS writing is 10% faster than the standard tilemanager, even though there is the additional compression compute required. When lossy compression is used the impact is positive, meaning that the additional complication of detrending the lossy data before encoding slows the processing by about 30%. The expectation is that using the GPU code would remove or reduced this penalty.
3.4 Compression ratios vs distortion
The histograms of the image residuals (Figure 7) indicate the quality of the reconstructed data after imaging. MGARD produces a consistent residual distribution in terms of distribution shape and the maximum residuals are consistent with the specified error bounds during compression. The residuals produced from the images compressed with relative error bounds surpass one standard deviation for the
$10^{-4}$
and
$10^{-5}$
cases, and those for the absolute error bounds surpass
$1\sigma$
for the
$10^{-3}$
case. That is, the distortions introduced by the compression, will start to be detectable against the noise levels in the images. The images shown in Figure 3 show that majority of these high residuals are situated in the corners of the image which are cut off or normalised out during regular imaging.
The distribution of residuals between the images produced with lossy compressed data and the ground truth image. The colours are as for Figure 3. The residuals are normalised by the image RMS of each 100-channel cube. In the first two instances (absolute error bound of
$10^{-3}$
and relative error bound of
$10^{-4}$
) a significant number of pixel differences are greater than 1.0
$\sigma$
, the image RMS.
4. Discussion and conclusions
4.1. Impact on data volumes
The compression of the, largely empty, visibility and PSF grids is significant. The incremental compression by introducing lossy, compared to lossless, compression would allow for even smaller storage requirements. Compression levels of up to 20-fold do not introduce detectable errors in the reconstructed images. However, as there are associated risks to a multi-year observing project in using lossy compression we have decided for the DINGO project to accept the 7-fold compression we achieve with lossless compression. The main driver for this project was to ensure we could compress the data sufficiently to be able to complete the survey within our Pawsey storage allocation. Lossless uv-gridding achieved this goal, and has zero risk of degradation of the data products. Thus going-forward this is the option we have adopted.
4.2. Impact on I/O speed
We demonstrated the sensitivity of the processing to parallel writing in the comparison of the serial casa writer and the parallel fits or adios writer. We note that the casa storage manager has recently been improved in ASKAPsoft, however we have not measured the performance of this new manager. We still favour the ADIOS library-based solution rather than a bespoke solution, firstly for the ease of code maintenance and also for the inclusion of the advanced features. We have found a direct impact on the total processing time with the compressed datasets; a 7-fold compression leads to a 7-fold reduction in the writing time. These results will depend on the degree to which the tasks are I/O bound, the hardware of the compute, and also to degree of compression applied. However the reduction in processing time allows for exploration of more imaging parameters and thus a better science result.
4.3. Impact on compute
The computational overhead for the lossless compression is small, and is actually compensated for by adios I/O efficiencies. Lossy compression using the CPU imposes a small processing time overhead; the experience is that for GPU implementation these are even smaller (Gong et al. Reference Gong2023). With the resources available we were not able to explore that aspect, and we are not using the lossy compression mode in DINGO, so this was not investigated further.
4.4. Impact on image quality
We confirmed that, as expected, data that has had lossless compression applied reproduced the non-compressed results perfectly. We have explored the image of the lossy compression for DINGO and find, for moderate compression error bounds, the impact is minimal. At larger compression error bounds (e.g. a relative error bound of
$10^{-3}$
) the image reconstruction starts to diverge. This we believe because of the degradation of the accuracy in the image-domain PSF when it is used for the CLEAN minor cycles. This is shown by the growth of the pixel error distribution shown in Figure 4. These errors build up and dominate the reconstruction. At lower error bounds the PSF is sufficiently accurate to allow for good deconvolution and image reconstruction.
4.5. Consequences for SKA data management
The SKA will, and have, invested a significant amount of capital into the storage, transmission and computing infrastructure for both nodes of the SKA. The tests described in the preceding sections have been specifically designed around the bespoke nature of deep imaging with ASKAP. Preliminary tests have shown that applying MGARD to raw SKA data (simulated) will yield similar factors of compression, the inclusion of which would simplify the implementation of storage and transmission solutions considerably for the SKA (Dodson et al. Reference Dodson2025). The processing done for SKA and ASKAP to produce images is almost identical, meaning that implementing ADIOS2 parallel I/O in the SKA pipelines would only improve the efficiency and speed of processing. We would recommend that this implementation is done in memory, that is, that the xarray storage is extended to include MGARD compression in-situ, so that the requirements on both the physical memory and the memory bandwidth can be reduced.
4.6. Conclusion
We have successfully demonstrated the application of grid stacking observations for deep
${\mathrm{H}\,{\scriptstyle \mathrm{I}}}$
observations using MGARD compression, via the ADIOS2/CASACORE framework, within standard radio astronomy software. The integration with ADIOS2 provides a parallel read/writer with performance that matches the parallel FITS writer. The combined impact will be on improved I/O performance, and thus pipeline run times, due to the reduced dataset sizes. Comparing the images made with various compression approaches allows us to quantify the impact. We show that we can apply lossy compression to the data files and achieve compression ratios up to 20, using well-defined error bounds, without impacting the results. Nevertheless, out of an abundance of caution, we have settled on the safer lossless compression option as that fulfils the Pawsey storage requirements for DINGO. In addition, the parallel reading and writing provided by ADIOS offers an additional improvement in I/O, and this is readily integrated with most common Radio Astronomy applications via the CASACORE libraries. Compression of the visibilities offers an attractive solution to the SKA I/O challenge, and we have demonstrated that the MGARD approach to compression can guarantee that the data is not degraded.
Acknowledgements
This scientific work uses data obtained from Inyarrimanha Ilgari Bundara/the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. CSIRO’s ASKAP radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42). Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy. ASKAP uses the resources of the Pawsey Supercomputing Research Centre. Establishment of ASKAP, Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory and the Pawsey Supercomputing Research Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund.
Part of this research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013, and in part by the Pawsey PACER project, and by the U.S. Department of Energy (DOE) RAPIDS-3 SciDAC project under contract number DE-AC05-00OR22725.
Some of this work was supported by the Australian SKA Regional Centre (AusSRC), Australia’s portion of the international SKA Regional Centre Network (SRCNet), funded by the Australian Government through the Department of Industry, Science, and Resources (DISR; grant SKARC000001). AusSRC is an equal collaboration between CSIRO, Curtin University, the Pawsey Supercomputing Research Centre, and the University of Western Australia.
This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), Matplotlib (Hunter Reference Hunter2007), and Numpy (Harris et al. Reference Harris2020).
Data availability statement
There is no proprietary period for ASKAP Survey Science data products, image cubes are available from the CASDA archive but visibilities are not retained long term. The visibility grids for these data products are available on request to MM, and will eventually be accessible on CASDA. The scripts for processing and data-analysis underlying this article will be shared upon reasonable request to the corresponding author, RD.
∼10−5
ℓ
ℓ
10−3
10−4
σ
Open access
