Super-resolution maps of icy moons: Application to the Cassini-VIMS observations of Enceladus with updated navigation and photometric correction

Abstract

The icy moons of the outer solar system have been found in the forefront of planetary exploration, due to their diverse geomorphology coupled with numerous discoveries hinting to possible habitable conditions. Of paramount importance for constraining the internal structure and properties of these bodies are high-precision cartographic products, allowing for linking surface morphologies and composition to internal activity. These cartographic products are constructed by reconciling numerous different observations, which, due to the complexities of an orbital tour, are recorded under a great range of different illumination and observation conditions, with varying spatial sampling. The aim of this study is to provide a complete pipeline for producing Super-Resolution (SR) maps for icy moons, by combining all the overlapping observations for a specific target. The development and testing of our pipeline are carried out for Enceladus, using the hyperspectral cubes recorded by the Visual and Infrared Mapping Spectrometer (VIMS) instrument, on board the Cassini Saturn orbiter. We first navigate the data, calculating accurate projections that can be combined into global mosaics of Enceladus. We then implement a global photometric correction model that compensates for the different acquisition conditions. Finally, we combine all the overlapping observations to calculate the new SR maps of Enceladus in the infrared. The new SR global map of Enceladus provides new insight into the moon’s geomorphology as revealed from its surface spectral heterogeneity at finer spatial resolution, demonstrating the capability of SR to augment our current understanding of icy moons and help prepare for future observations.

1 Introduction

Despite its relatively small size, Enceladus is a geologically active icy moon, with observations confirming the existence of active jets consisting mostly of water vapor and ice grains emanating from the so-called “tiger-stripes”, a complex of tectonically active fissures at the moon’s South Polar Terrain (SPT) (Porco et al., 2006; Spencer et al., 2006; J. H. Waite et al., 2006). The detection and analysis of the cryovolcanic activity on Enceladus hints at the existence of an internal ocean, enriched with volatiles and organics, in contact with the moon’s rocky core (Hansen et al., 2011; Peter et al., 2024; Postberg et al., 2011; J. Waite et al., 2009; J. Waite et al., 2017). All these characteristics make Enceladus a prime target for the search of habitable environments within our solar system, a status solidified by the report “Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032”, a road map for National Aeronautics and Space Administration (NASA) future missions, which identified an Enceladus mission as a priority both in the medium-class New Frontiers 5 and 6 calls (Enceladus Multiple Flybys mission) and as the second-highest priority in the large Flagship mission category (Enceladus Orbilander mission) (National Academies of Sciences & Medicine, 2022). On the European Space Agency (ESA) side, the “Voyage 2050 Final recommendations from the Voyage 2050 Senior Committee” report cites Enceladus as a prime candidate for a future astrobiological mission and underlines the potential for a dedicated mission within the Voyage 2050 time frame (Tacconi et al., 2021), which will be ESA’s fourth Large-class mission (L4).

The vast majority of our current knowledge for Enceladus comes from the Cassini-Huygens mission, the first dedicated mission to study Saturn’s system and one of the most ambitious missions in the history of robotic exploration of our solar system (Matson et al., 2002). Launched in 1997, the spacecraft was inserted into Saturn’s orbit in 2004. Planned for an initial duration of 4 years, the mission was finally extended until September 2017. The Cassini-Huygens mission was designed to explore the Saturn system, including the icy moons, and among them Enceladus. One of the scientific instruments employed for the study of the Saturnian system was the Visual and Infrared Mapping Spectrometer (VIMS) (Brown et al., 2004), which performed close observations of Enceladus during 23 targeted flybys throughout the 13 years of the Cassini-Huygens mission. The VIMS instrument acquires hyperspectral images using two independent detectors: one for the visible and one for the infrared (IR) part of the spectrum (with a total of 352 spectral channels ranging from 0.4 to 5.1 \(\mu \)m), providing not only information on the composition of the surface (e.g. Brossier et al., 2018; Clark et al., 2010; Griffith et al., 2019; McCord et al., 2008), but also the physical state of the surface materials (e.g. Brown et al., 2006; Filacchione et al., 2016; Jaumann et al., 2008; Newman et al., 2008; Robidel et al., 2020; Taffin et al., 2012). In this work we will only consider the VIMS-IR instrument, which provides a spectral range from 0.85 to 5.1 \(\mu \)m in 256 bands.

The precision of compositional mapping and its correlation with geological features are often limited by the spatial resolution of the captured images. Indeed, spatial resolution in hyperspectral images acquired by a spacecraft performing repeated flybys or in high-altitude orbits is highly variable and, in the case of the VIMS instrument, rarely better than a few kilometers, which strongly limits the geological interpretation. In the general case, high-resolution (HR) imaging data are possible during close flybys or at low-altitude orbits but will always be limited by the bandwidth, hence data volume (Xie et al., 2021). Having optimized to the best of our capacities the image acquisition from the hardware point of view, we are left with post-mission signal processing techniques in order to enhance the resolution of the already acquired images. One such solution is the implementation of Super-Resolution (SR) image reconstruction, defined as a broad set of techniques used to extract a HR image from a set of Low-Resolution (LR) observations, finding application in various fields such as satellite and aerial imaging (Shi et al., 2017; Zhang et al., 2014), medical image processing (Kennedy et al., 2006; Van Reeth et al., 2012), text image analysis (Capel & Zisserman, 2000), biometrics recognition (Alonso-Fernandez et al., 2018; Nguyen et al., 2018) and many more (Nasrollahi & Moeslund, 2014). SR methods output a HR image with a higher pixel density than the original image, that is not just an oversampled version of it, by retrieving the high frequency components expressed as non-redundant information due to the sub-pixel shifts between the LR images.

In the context of planetary observations, SR has been implemented in a variety of observations of Mars, either by applying the traditional multiple-image reconstruction model based on the Irani-Peleg equation (Tao & Muller, 2016), or by employing generative neural networks to perform single-image SR (Sun & Chen, 2022; Tao et al., 2021, 2022; Wang et al., 2021). At this point, it should be noted that Mars poses an exception in the exploration of the bodies in our solar system, in the sense that it has had a great number of dedicated missions resulting in a plethora of available data of high quality, due to the combination of orbital and landing components. However, this is not the case for the high astrobiological interest icy moons, with only a handful of dedicated orbital missions with relatively small data sets available. The SR techniques found in the literature are specifically tailored to leverage the advantages of their respective data sets and therefore present three main limitations when attempting to implement them in the case of icy moons: 1) they require a significant number of repeated passes of the same area (for multiple-image SR) or vast training sets (for single-image neural networks) that are not available for icy moons, 2) they rely on relatively homogeneous data sets whereas the available data sets for icy moons have been acquired during flybys, under a variety of observation geometries, illumination conditions and original sampling that need to be reconciled, and 3) they use convolution matrices (also known as masks, which assign a new brightness value to each pixel based on its brightness and those of its neighbors) that consist an image processing technique possibly introducing uncontrollable artifacts that do not correspond to real surface features.

From all the above, it becomes clear that we need to devise a completely new strategy for applying SR to icy moons. Inspired by previous work (Karkoschka et al., 2017; Karkoschka, 2015) and using as a starting point the latest VIMS Enceladus global mosaic included in the pivotal work of Robidel et al. (2020), the main objective of this study is to develop a general SR mapping method for icy moons, in order to produce surface maps with increased spatial resolution. At first we describe the data set used to develop and test our new SR methodology before moving into the detailed description of the three fundamental steps of our methodology: navigation, photometric correction and computation of the SR map. We present a new pipeline in Python for creating the navigation cubes necessary in merging individual cubes into global mosaics, increasing the accuracy of the covered area and maximizing the number of pixels used. We then determine a global surface photometric function which mitigates variations in the brightness due to the changing geometric conditions of the observations. Finally, we implement a novel SR algorithm, which combines the information from all the overlapping pixels over a given area to acquire products with increased spatial resolution. We finally conclude the presentation of our new tool by applying it to the mapping of Enceladus with the VIMS instrument and discussing the potential of using this new, SR global map to reveal geological and compositional units at finer spatial scale.

2 The VIMS data set for Enceladus

During its 13-year long mission, Cassini performed 23 targeted flybys of Enceladus (from E0 to E22), with flyby E9 not recording VIMS data. These flybys have recorded 3982 hyperspectral cubes, from which 3911 have been calibrated and included in the Nantes Cassini VIMS Data Portal at https://vims.univ-nantes.fr/ (Le Mouélic et al., 2019a), accounting for a total of 1 582 502 pixels. The remaining unavailable 71 targeted flyby cubes correspond to two cubes (1500059679_3 in flyby E2 and 1829237951_15 in flyby E22), for which the IR channel is not available and 69 noodles (single-sample cubes) in flybys E5 and E22, that require special treatment. Moreover, there were 124 non-targeted flybys of Enceladus, during which there were 7502 additional cubes recorded. In total, for Enceladus there are 147 flybys, containing a total of 11 502 calibrated cubes, available in the Nantes VIMS Data Portal.

The creation of a global mosaic for Enceladus can be achieved by a careful selection of the best pixels (most resolved and under the most favorable observation conditions). However, the implementation of the SR approach requires a delicate balance between competing incentives: excluding from the data set unsuitable pixels, but at the same time including enough overlapping observation to allow for the successful calculation of the SR map. This is achieved by filtering the available data down to the pixel level. Before that and as a first rough selection, we isolate suitable data at the cube level, before moving forwards with our pipeline and a more detailed selection at the pixel level. This step also guarantees that we will not be employing our computational resources for calculations on cubes that will be eventually rejected as a whole at the pixel level filtering. The Nantes VIMS Data Portal allows for the visual inspection of each cube separately, along with their metadata, which facilitates the swift identification of the cubes in the targeted and non-targeted flybys that meet the following criteria: 1) mean spatial sampling \(\leq \) 60 km/pixel, 2) minimum incidence angle \(\leq \) 80\(^\circ \), and 3) minimum emergence angle \(\leq \) 80\(^\circ \). The filtering of the cubes at the sampling limit of 60 km/pixel, is selected so that we have enough range for follow-up testing of gradually increasing maximum sampling, while remaining in the range of cubes that are resolved enough to be useful in the context of our maps. It should be noted that this is an indication of the mean spatial sampling of the whole cube and therefore at the level of individual pixel we can have sampling values greater than that limit. Moreover, as we present in Section 3.1, we implemented a different approach to calculate the pixel spatial sampling. The photometric angles are constrained following Robidel et al. (2020) in order to respect the limits of the photometric model that is implemented in the respective section that follows. At this first stage of the data selection process, we create a data set consisting of 1198 cubes from the targeted flybys accounting for 526 811 pixels and 364 cubes from 13 non-targeted flybys accounting for 147 784 pixels. The data set is complemented by five cubes targeting Enceladus not available in the Nantes VIMS Data Portal, but which are nonetheless used by Robidel et al. (2020): 1604167154_1, 1635842860_1, 1635842868_2, 1669819057_1 and 1694646566_1. These five cubes are calibrated manually using the third version of the Integrated Software for Imagers and Spectrometers (USGS-ISIS3) and contribute 11 792 additional pixels to the data set. Even with a first order selection at the cube level the majority of the pixels in the data set corresponds to the sky (349 768 pixels), while 53 003 pixels contain the target’s limb. The total amount of usable pixels, namely those that are fully contained within the target, is 271 824 pixels. This data set is ingested by the SR pipeline we have developed, which is presented in the following sections and is comprised of three main steps: 1) navigation, 2) photometric correction and 3) super-resolution (SR).

3 Super-resolution (SR) global mapping of Enceladus with VIMS

3.1 Step 1: Filtering and navigation of the VIMS data set

The first step in our pipeline for creating SR mosaics for icy moons is to create a module for calculating the navigation cube for each of the hyperspectral cubes. This transformation allows us to convert two-dimensional images at a specific wavelength (Fig. 1a) into coverage patches on the surface of the target (Figs. 1b and 1c). The navigation allows in particular for the calculation of the longitude, latitude and conditions of observations (illumination and observation angles, namely incidence, emergence and phase angles, referred hereafter as the photometric angles) of a given VIMS pixel in the planetographic reference frame of a given target. Even though navigation cubes already exist for the Cassini VIMS data set (https://vims.univ-nantes.fr/), we implemented our own approach in order to ensure continuity in the pipeline and increased navigational accuracy for the purposes of the registration step of SR methodology. The need for an accurate navigation of the hyperspectral cubes is dictated by the incentive to create the best mosaics possible for two reasons: 1) the more accurate the projection (within intrinsic ephemeris and spacecraft pointing error, the evaluation of which is discussed by Nicholson et al., 2020), the better the matching of the overlapping observations in the global mosaic and finally the higher the quality of the photometrically corrected mosaic with regards to visible seams; and 2) as we aim to combine overlapping pixels of the data set to create new pixels in our SR methodology, accurate projection of those overlapping pixels is of the utmost importance in order to make sure that we are combining spectral information from the same area.

The first step for the calculation of the navigation cube is to extract basic information from the calibrated cube (.cub file) header, namely the samples and lines dimensions, the name of the target, the start and stop time of the cube acquisition, the interline delay, the X- and Z-axis offsets, the sampling mode (“NORMAL” or “HI-RES”) and the exposure duration. Considering that the infrared channel of the VIMS instrument (VIMS-IR) operates in “whiskbroom” mode, hence recording the spectrum of only one spatial pixel per exposure, the essence of the navigation cubes lies in the calculation of the time at which each pixel has been recorded, by using the following cumulative time equation (as described by Brown and Soderblom (n.d.) and in the algorithm of Seignovert et al., 2023):

\begin{equation} \Delta t_{row}^{IR} = i \times \Delta t_{integ}, i \in [0, s_{row}] \tag {1a} \end{equation}

\begin{equation} \Delta t_{line}^{IR} = j \times (\Delta t_{delay}+s_{row} \times \Delta t_{integ}), j \in [0, s_{line}] \tag {1b} \end{equation}

\begin{equation} \Delta t_{pixel}^{IR} = \Delta t_{row}^{IR} + \Delta t_{line}^{IR}\tag {1c} \end{equation}

where \(\Delta t_{row}^{IR}\) the time elapsing between rows, \(i\) is the number of the row where the pixel is located, \(\Delta t_{integ}\) is the corrected integration time, \(s_{row}\) is the number of rows contained in the cube as read from the header, \(\Delta t_{line}^{IR}\) the time elapsing between lines, \(j\) is the number of the line where the pixel is located, \(\Delta t_{delay}\) is the delay between lines, \(s_{line}\) is the number of lines contained in the cube as read from the header, \(t_{pixel}^{IR}\) is the time elapsing between the start time of the cube and the acquisition of the specific pixel. The “ExposureDuration” and “InterlineDelayDuration” time values extracted from the cube’s header, are corrected by applying a correction coefficient of 1.01725 to compensate for the drift of the instrument’s clock (Brown & Soderblom, n.d.) (the start time of the acquisition does not require the same correction since it is measured by the spacecraft’s clock). Finally, the time of a single pixel acquisition (assuming instantaneous acquisition at half the integration time) is:

\begin{equation} t_{pixel} = t_{start} + \Delta t_{pixel}^{IR} + 0.5 \times \Delta t_{integ}.\tag {2} \end{equation}

By combining the time equation with the geometry of each pixel within the instrument’s Field Of View (FOV) and using this information as input for the SPICE (Spacecraft Planet Instrument C-matrix Events) Application Program Interfaces (APIs) (Acton, 1996), along with the relevant kernels containing information on the position and orientation of both spacecraft and target, we can calculate a “pinhole” projection of each point comprising the edges of each pixel and in this manner produce the navigation cube corresponding to each one of the calibrated cubes, which we can then visualize as 3D projections on the target represented in light blue (Fig. 1b) or by an equirectangular projection (Fig. 1c). The cube 1487302144_1 (flyby E0 and acquisition start time 2005-02-17T03:02:38.846) has been selected as an example to showcase how the lines of a VIMS-IR image do not necessarily align after the navigation of the cube, due to the relative movement of the spacecraft with respect to the surface of the target and the “whisk broom” scanning acquisition mode. As far as the selected number of points per pixel edge is concerned, we have deduced empirically that five points per edge provides a good balance between the accurate projection of the outline of a pixel and the computational time (less than 1 minute per cube).

Figure 1: (a) 2D image of cube 1487302144_1 at standard wavelength 1.80401 \(\mu \)m, (b) the corresponding projection in 3D of its pixels (in red) on the surface of Enceladus (in light blue) and (c) equirectangular projection, with the colorbar representing the I/F values at 1.80401 \(\mu \)m.

Regarding the shape of the target body, used in the calculation of the intercepts and photometric angles (specified in the SPICE APIS with the keyword method), there are two options available: 1) “ELLIPSOID” where the target body is modeled as a triaxial ellipsoid and 2) “DSK/UNPRIORITIZED” where a Digital Shape Kernel (DSK) is used, an integrated shape model subsystem for Digital Shape Models (DSMs). As far as Enceladus is concerned, a DSM has been constructed by Bland et al. (2020), from a dense photogrammetric control network using 625 Cassini Imaging Science Subsystem (ISS) images of Enceladus. The existence of a DSM is not a prerequisite for implementing our SR methodology but allows for a more realistic representation of the surface of the target and therefore more accurate projections of the pixels and more accurate values for the incidence and emission angles associated with each pixel (difference in phase angle is zero because the phase angle is independent of the surface normal vector that changes when transitioning from the triaxial ellipsoid to the DSM). From our investigation we deduce that the difference in the incidence and emission angles between the two models can be as high as around \(\pm 10{^\circ }\) (as it can be seen in more detail in the supporting material Appendix A).

We also need to address the calculation of the pixel sampling, which can not be defined explicitly anymore in the case of non-square pixels. The calculation, on first order, can be carried out using the following simplified formula \(s_p=\theta d\), where \(d\) is the distance between the center of the pixel and the spacecraft in km and \(\theta \) is the angular size of the pixel in rad/pixel. The VIMS-IR has two spatial resolution modes: “NORMAL” with 1\(\times \)2 spatial pixels binning resulting in Instantaneous Field of View (IFOV) of 0.5\(\times \)0.5 mrad and “HI-RES” where full detector resolution of 0.25\(\times \)0.5 mrad per IFOV (no binning) was applied (Mccord et al., 2004). Therefore, for mode=“NORMAL” the angular pixel size is \(\theta =0.5\) mrad/pixel and for mode=“HI-RES” it is \(\theta =0.25\) mrad/pixel (Clark et al., 2018). This approach in calculating the pixel spatial sampling assumes square pixels at nadir and tends to underestimate the pixel spatial sampling of “elongated” pixels (namely pixels close to the limb at HI-RES mode) or pixel-polygons, resulting as common overlapping areas in the SR part of the pipeline. The decreasing spatial sampling order at which the pixels are stacked does not influence the SR calculation, which consists of a linear combination of the overlapping pixels independently of their order, but affects the original and photometrically corrected mosaics, which present at the topmost layer the most resolved observations and the filtering of pixels based on their spatial sampling values. In order to take into account this limitation, we introduce a more flexible and accurate metric for the pixel spatial sampling, by calculating all the possible distances from each edge of the pixel projected polygon to every other edge (including the diagonals) and assigning as sampling of the pixel the largest of these distances as the most conservative approach. All distances are calculated in kilometers for the 3D case (flat polygons in 3D space) and therefore the spatial sampling values are independent of the projection employed. A comparison of the two different methods to calculate the pixel sampling (diagonal sampling versus distance sampling) can be seen in Figure 2, with the lower limit for the pixels acquired in “NORMAL” mode shown in red (square pixels with the ratio of the diagonal over the edge being equal to \(\sqrt {2}\)) and the “HI-RES” mode in orange (rectangular pixels with the large edge being twice the length of the small edge have a ratio of the diagonal over the small edge equal to \(\sqrt {5}\)). This figure clearly demonstrates that the commonly used calculation based on distance can largely underestimate the VIMS pixel sampling compared to our approach, which can assign spatial sampling values as high as 20 times larger.

Figure 2: Comparison for each target pixel in our selected data set of the sampling calculated with the simplified distance formula (distance sampling) versus the proposed evaluation of all the lengths of all the edges and the diagonals of each pixel-polygon (diagonal sampling). The red dashed line marks the lower limit for the pixels acquired in “NORMAL” mode and the orange dashed the lower limit for the pixels acquired in “HI-RES” mode.

Finally, we end up with navigation cubes for each of the VIMS cubes in our selected data set, containing the following geographic and geometric information for each pixel: longitude and latitude of the center and the corners, time of acquisition, photometric angles and pixel sampling. In Figure 3 we present the global mapping of the incidence (Fig. 3a), emission (Fig. 3b) and phase angles (Fig. 3c), as well as the pixel sampling (Fig. 3d) for each of the VIMS cubes used in this study. Even though global maps of Enceladus should follow the planetographic west-positive longitudes (0\(^\circ \) to 360\(^\circ \)W) system in accordance with International Astronomical Union (IAU) instructions (Archinal et al., 2018), in this work we present all the maps in the planetocentric 180\(^\circ \)W to 180\(^\circ \)E reference system to be consistent with other cartographic products employed for validation and reference (Bland et al., 2018; Crow-Willard & Pappalardo, 2015). The incidence and emission angles have been filtered to exclude values greater than \(80{^\circ }\) so as to avoid extreme observation and illumination conditions, as it will be explained in more detail in the next section. The sampling values are presented for pixels up to 50 km/pixel so as not to saturate the top layer of the map (pixels are sorted by increasing spatial sampling and plotted one on top of the other over the same areas, hence the top layer of the map includes the most spatially resolved observations), which shows the most resolved observation. In the global maps of the incidence and emission angles, we can observe the variation within the same cube, which is attributed to the topography stemming from the incorporation of the DSM of Enceladus.

Figure 3: Global maps of (a) incidence, (b) emission, and (c) phase angles for the selected data set, and (d) pixel spatial sampling up to 50 km/pixel.

Having calculated an accurate navigation of the pixels in the data set, we can now create global mosaics of Enceladus by assigning the relevant spectral values from the calibrated cubes to each pixel-polygon. However, this requires two additional processing steps of the data set in the spectral dimension of the cubes: 1) correction of wavelength shift and 2) removal of saturated parts of the spectra. To address the wavelength shift documented by Clark et al. (2018), we calculate a linear interpolation and resample it on a standard reference wavelength grid defined at the beginning of the mission, as already implemented on VIMS cubes of Titan by Le Mouélic et al. (2019a). To filter the saturated part of the data set, we implement a novel approach by taking into account the background noise recorded in each of the cube’s backplanes, which we subtract from the total capacity of the detector to define a dynamic saturation limit for each spectrum as a function of wavelength. The details for the wavelength shift and saturation corrections can be found in Appendix B.

The final product at the end of this first step of the methodology is presented in Figure 4 in the form of a global mosaic of the selected VIMS observations of Enceladus at standard wavelength 1.80401 \(\mu \)m (a characteristic reflectance peak in the water-ice spectrum) in (Fig. 4a) equirectangular projection, (Fig. 4b) orthographic projection centered at the north and (Fig. 4c) south pole respectively, once the entire data set has been corrected for the wavelength shift, signal saturation, and properly navigated.

Figure 4: Global mosaic of the VIMS observations of Enceladus at wavelength 1.80401 \(\mu \)m in (a) equirectangular projection and orthographic projection centered at the (b) north and (c) south pole, after the application of the global filtering at the pixel level and data processing described in the text (including dayside surface pixels selection, wavelength shift correction, saturated wavelengths removal and navigation).

3.2 Step 2: Production of photometrically corrected global map of Enceladus

The first step in producing global mosaics of the surface of Enceladus has been achieved by calculating the new navigation cubes and projecting each pixel of each observation from lower to higher pixel spatial sampling (in order to have the most resolved observations on top) at the surface of Enceladus. This process is independent of the wavelength depicted and therefore can be used for mapping any of the 256 VIMS-IR channels (or combination thereof). The global map produced in this way (Fig. 4) does not accurately represent the surface of Enceladus, since it does not appear consistently homogeneous, with striking artifacts such as visible seams at the areas where two different observations overlap. These differences are attributed to the fact that different observations of the same area happened at different stages of the orbital tour, presenting a variety of observation and illumination conditions. In order to mitigate the different illumination and geometric parameters of each observation we have to implement a photometric correction method as a second step, which normalizes all the observations to zero incidence, emission and phase angles. This allows us to merge, in the third and final step of our methodology, overlapping observations under the same observation conditions, thus with comparable and compatible radiance factor values.

The photometric correction of our data set requires the implementation of a radiative transfer model for surface regolith of atmosphereless bodies. Among numerous models that have been implemented, there exist two scattering theories for solid planetary surfaces that dominate the literature due to their accuracy and relative simplicity to implement: the semi-empirical Hapke model (integral) (Hapke, 1981, 1984, 1986, 2002, 2008, 2012; Hapke & Wells, 1981) and Kaasaleinen-Shkuratov model (statistical). Since the photometric correction of the data is an intermediate, but necessary, step in our methodology, we decided not to investigate in depth different models, but instead to rely on the seminal work of Filacchione et al. (2018b, 2018a) and Filacchione et al. (2022) and Robidel et al. (2020) that already demonstrated the efficiency of the Kaasaleinen-Shkuratov model for Enceladus, implementing a method proposed by Shkuratov et al. (2011), where the photometrically corrected radiance factor (\(r_{f_{model}}\)) is modeled in two parts: the disk function \(D(i,e,g)\) and the equigonal albedo \(A_{eq}(g,\lambda )\):

\begin{equation} r_{f_{model}}=\left ( \frac {I}{F} \right )_{model}(i,e,g,\lambda ) = D(i,e,g)A_{eq}(g,\lambda )\tag {3} \end{equation}

where \(i\), \(e\) and \(g\) are the local incidence, emission and phase angles respectively and \(\lambda \) the wavelength. Different disc functions have been evaluated by Robidel et al. (2020) and it was determined the best results, in terms of the quality of the mosaic produced, are achieved by implementing the model by Akimov (Kaidash et al., 2009):

\begin{equation} D(i,e,g)=D(g,\beta ,\gamma )=\cos \left ( \frac {g}{2} \right )\cos \left [\frac {\pi }{\pi -g}\left (\gamma -\frac {g}{2} \right )\right ]\frac {(\cos \beta )^{\frac {g}{\pi -g}}}{\cos \gamma }\tag {4} \end{equation}

where \(\gamma \) and \(\beta \) are the photometric longitude and latitude respectively:

\begin{equation} \gamma =\arctan \left ( \frac {\cos i-\cos e \times \cos g}{\cos e \times \sin g} \right ),\tag {5} \end{equation}

\begin{equation} \beta =\arccos \left ( \frac {\cos e}{\cos \gamma } \right ).\tag {6} \end{equation}

Having calculated the navigation cubes for all the available observations, we have available for each pixel the \(\frac {I}{F}(\lambda )\) and the local photometric angles \(i\), \(e\) and \(g\), including the impact of local topography. Due to the fact that a VIMS pixel can contain multiple facets of the triangular tessellated model comprising the DSM, one VIMS pixel is associated with the mean of the contained topography. In order to address the variability of the photometric angles within one pixel, we calculate the incidence, emission and phase angle associated with each pixel as the mean of the respective values at the pixel corners. We can then correct for Enceladus’ surface photometric behavior by calculating the equigonal albedo per pixel and wavelength, such as:

\begin{equation} A_{eq}(g,\lambda )=\frac {\frac {I}{F}(i, e, g,\lambda )}{D(i,e,g)}.\tag {7} \end{equation}

For every wavelength we calculate the coefficients \(a_0\) and \(a_1\) for a linear regression (simplest but robust approach as demonstrated by Robidel et al. (2020)) of all the points in the graph of the equigonal albedo (\(A_{eq}\)) as a function of the phase angle (\(g\)), also called the phase curve:

\begin{equation} A_{eq}(g,\lambda ) = a_1(\lambda ) g + a_0(\lambda ).\tag {8} \end{equation}

Then, the photometrically corrected values of the radiance factor at reference geometry \(i=e=g=0{^\circ }\) can be calculated by implementing the model-agnostic normalization by Xu et al. (2020) as follows (\(D(0,0,0)=1\)):

\begin{equation} \label {eq:rf_cor} r_{f_{cor}} = \left ( \frac {I}{F} \right )_{cor}(0, 0, 0,\lambda ) = r_f(i,e,g,\lambda )\frac {r_{f_{model}}(0, 0, 0,\lambda )}{r_{f_{model}}(i,e,g,\lambda )}=\frac {r_f(i,e,g,\lambda )}{D(i,e,g)\left [1+\frac {a_1(\lambda )}{a_0(\lambda )}g\right ]}.\tag {9} \end{equation}

In order to ensure that we are implementing the above-described model within an acceptable range of application (observation under extreme conditions are not eligible) and at the same time maintain global coverage and sufficient overlapping of the pixels to apply SR at the next step, we have further filtered our data set at the pixel level. With regards to the incidence and emission angles, we follow the approach of Robidel et al. (2020), where \(i, e \leq 80{^\circ }\). With regards to the values of the radiance factor, we filter using a lower and an upper limit \(0.001 \leq \frac {I}{F} \leq 1\), to achieve sufficient coverage in all wavelengths, including those with fainter signals. Similarly after testing gradually increasing upper limits for the spatial sampling of the pixels, we conclude that \(s_p \leq 78\) km/pixel achieves a good balance between sufficient coverage, number of overlapping observations and the core assumption in our SR approach that each pixel is small enough to represents a spectrally homogeneous part of the surface. At the end of this pixel-level filtering our data set is reduced to 156 294 pixels (at standard wavelength 1.80401 \(\mu \)m. In Figure 5, we create a scatter plot of the equigonal albedo (\(A_{eq}\)) as a function of the phase angle (\(g\)), for all the pixels in the finalized data set. The red line represents the best linear fit of the data (phase curve) while the color scale quantifies the Point Density Function (PDF), that indicates the relative data density at each point. The coefficients of the best linear fit for all wavelengths can be found in the supporting material (Appendix C).

Figure 5: Equigonal albedo (\(A_{eq}\)) as a function of the phase angle (\(g\)) for wavelength 1.80401 \(\mu \)m, with the best linear fit shown in red. Color scale represents the Point Density Function (PDF).

Using the coefficient of the slope of the best linear fit in Equation (11) for every pixel in our data set, we calculate the new photometrically corrected radiance factor value for each pixel and construct the photometrically corrected global mosaics of Enceladus, presented in Figures 6a, 6b and 6c. The calculation and propagation of the radiometric error is discussed in the supporting material (Appendix E). It is therefore possible to compare with the uncorrected mosaics (Figs. 4a, 4b and 4c) and appreciate the significant improvements of the photometrically corrected ones.

Figure 6: Photometrically corrected global mosaic of the VIMS observations of Enceladus at wavelength 1.80401 \(\mu \)m in (a) equirectangular projection and orthographic projection centered at the (b) north and (c) south pole.

Comparing the uncorrected (Fig. 4a) and photometrically corrected (Fig. 6a) global mosaics of Enceladus, we can enumerate three qualitative points that at a first glance validate the photometric model we have implemented: 1) visible diminishing of artifacts-seams in the transition between two spatially neighboring observation (especially when these are under very different photometric conditions), 2) a more homogeneous global view, as it is expected from a target such as Enceladus, and 3) gain in contrast at the local scale, that better visualizes geomorphological features, as is the case for the tiger stripes in the south polar terrain. A quantitative validation of the photometric correction step relies on comparing the spectra of overlapping pixels between the original uncorrected map and the photometrically corrected map. This can be seen in Figure 7 with two overlapping pixels from two different cubes, recorded under different observation conditions:

• pixel 1: cube=1660437591_1, index=103, i=55.91\(^\circ \), e=63.09\(^\circ \), g=24.32\(^\circ \), sampling=26.94 km/pixel
• pixel 2: cube=1858921724_1, index=101, i=74.66\(^\circ \), e=36.65\(^\circ \), g=78.24\(^\circ \), sampling=29.52 km/pixel

As we can observe from Figure 7, before applying the photometric correction the spectra from the overlapping pixels are vastly different even though they describe similar areas on the surface of Enceladus with similar sampling, albeit under quite different observation conditions. After the photometric correction the two spectra present a better match (Fig. 7d), validating the photometric model in the spectral dimension. The matching of the spectra after the photometric correction is reflected in the spatial dimensions by the diminishing of the visible seams between the neighboring cubes. Of course the photometric correction is not perfect, with the limitations of the method manifesting as small diversions between the two spectra. These small differences can be attributed to the inherent limitation of the photometric model and also to the fact that we are performing the correction at the global level, even though the photometric behavior, like that of any planetary body, exhibits local variability. Additionally, the two overlapping pixels presented do not depict the exact same area, a fact that can also contribute to the small difference observed in the photometrically corrected spectra.

Figure 7: Global uncorrected (a) and corrected (b) map of Enceladus at 1.80401 \(\mu \)m with overlapping pixels cube=1660437591_1 (pixel 1), index=103 and cube=1858921724_1, index=101 (pixel 2) marked and spectra comparison of the uncorrected (c) and corrected (d) pixels along with their two sigma envelopes (the calculation of which is explained in Appendix E of the supporting material) and the residuals (bottom panels).

In the spatial domain it is also very interesting for validation purposes to present a transect comparison at an area with very high photometric angle contrasts, hence a very dramatic transition of the radiance factor values of neighboring pixels due to the different observation and illumination conditions. Such is the case for the transect extending from longitude 100\(^\circ \)E to 180\(^\circ \)E at latitude 1\(^\circ \)N (using 2000 equally spaced points) as it can be seen in Figure 8. Here we present the uncorrected (top left) and corrected (top right) global maps with the transect marked in red and green respectively, the radiance factor values of the two transects with the vertical dotted lines marking transition points from one cube to another (middle) and the values of the photometric angles along the transect (bottom). The efficiency of the photometric correction becomes obvious in this case, if we observe that the sharp changes in the reflectance value at transition points 101\(^\circ \)E, 106\(^\circ \)E and 173\(^\circ \)E (corresponding to seams in the global map) almost disappear and the original signal with very sharp dips is smoothed out. At the same time, within the same cube, the photometric model manages to correct the downwards trend of the reflectance factor transect of the uncorrected map (red curve), associated with the gradual increase of the incidence angle, while conserving and even augmenting the contrast of real features within a given cube. All the above validate the efficiency of the photometric model used to photometrically correct the Enceladus global mosaic in order to be ingested by the next part of the pipeline.

Figure 8: Uncorrected map with transect (latitude 1\(^\circ \)N, longitude 100\(^\circ \)E to 180\(^\circ \)E) in red (top left), photometrically corrected map with same transect in green (top right) and radiance factor values of transects with two sigma envelopes (middle), with dotted gray lines marking the transition between cubes. In the bottom panel we present the values of the incidence (blue), emission (orange) and phase (purple) angles along the transect.

3.3 Step 3: Weighted Average Super-Resolution (WASR) and final SR map ofEnceladus

The third and final step of our SR mapping pipeline consists of applying the SR calculation to the filtered, navigated and photometrically corrected data set. As we have described in Section 1, the SR approaches found in the literature cannot address the distinct nature of the problem in the context of the images collected during a complex and varied orbital tour and therefore we have to design our own merging and averaging technique to create global maps of augmented resolution. We start with the mathematical formulation of our new approach to SR, namely the Weighted Average Super-Resolution (WASR). We assume a pixel of random size and shape, with radiance factor \(r_f=\frac {I}{F}\) and area \(A_i\). The core assumption in our approach is that the radiance factor \(r_f\) is constant within the pixel. In the general case, each pixel covers an area of varying geological features and surface texture and composition, and therefore of varying spectral units, but since we have no sub-pixel information on their distribution, our assumption is the best possible case. The smaller the pixels used, the more homogeneous the area they represent is and the closer our assumption reflects the true nature of the problem.

The remaining question is how to combine the contribution of the various intersecting sub-section of all the overlapping pixels in order to calculate the radiance factor of the new SR pixel. The simplest approach would be to define the radiance factor of the SR pixel (\(r_{f_{SR}}\)) as a linear combination of the radiance factor of its overlapping pixels. After testing different models, we deduce that a simple and at the same time robust way to combine the radiance factor values of the overlapping pixels, while preserving the signal in the spectral dimension (SR should only have an effect on the spatial dimension) is to calculate the weighted average, hence the WASR method. The new radiance factor value of each SR pixel is then calculated by the following equation:

\begin{equation} \label {eq:sr} r_{f_{SR}} = \left ( \frac {I}{F} \right )_{SR} = \frac {\sum _{}^{} w_i r_{f_i}}{\sum _{}^{} w_i} = \frac {\sum _{}^{} w_i \left (\frac {I}{F}\right )_i}{\sum _{}^{} w_i}\tag {10} \end{equation}

where \(w_i\) is the weight associated with each overlapping pixel, defined as the ratio of the area of the original pixel contained in the SR element (henceforward called a “super-cell”) (\(a_i\)) over the total area of the pixel (\(A_i\)):

\begin{equation} w_i = \left ( \frac {a_i}{A_i} \right )^n\tag {11} \end{equation}

where the exponent \(n=1, 2, 3 ...\) allows for an empirical shifting of the weight contribution to allow for higher influence of the more resolved overlapping pixels, the higher the value of \(n\) is. This technique is implemented to address the diffusion of information observed for \(n=1\) in areas where one HR pixel is mixed with much bigger pixels. Higher values of \(n\) do not recover more super-resolved features, but rather allow for final cartographic products with higher contrast. The higher the value of \(n\), the more the most resolved overlapping pixel dominates over the contribution of the rest of the overlapping pixels approaching the photometrically corrected signal. Naturally, the higher values of \(n\) also increase the contrast of the residual from the photometric correction seams and therefore a value of \(n=2\) has been selected as a reasonable compromise.

In order to overcome the complexity of calculating overlapping geometries in 3D space, namely on the surface of a triaxial ellipsoid representing the target body, we are implementing the SR calculation on the 2D projection of the photometrically corrected map. Without loss of generality we chose an equirectangular projection. Areas on the target body contained in specific overlapping pixels, will be contained in the same pixels no matter the projection. Moreover, the calculated weights (\(w_i\)) are not affected by the area distortion potentially introduced by the selected projection since they are defined as ratios of the areas. This conclusion has been validated by performing the SR calculation both in 2D and 3D (where we calculate the area of spherical polygons) with the resulting maps being identical. The small differences in the SR radiance factor between the 2D and 3D case (less than 1%) are associated with the size of the polygon (small difference in area when considering flat versus spherical polygons) and show no correlation with latitude. Both cases (2D and 3D) have been incorporated in the corresponding part of the pipeline. The ideal theoretical approach to SR, would be to create a map, where each element is in the general case an irregular polygon, representing the common overlapping area of a distinct list of pixel from the data set, hence a distinct combination of \(\frac {I}{F}\) values. With a number of pixels in the order of \(10^5\) for Enceladus, the repetitive calculation of the common overlapping polygons becomes a computationally arduous process. Moreover, the complicated shape of the pixels along with their random orientation, can cause the repetitive intersection operations to produce sliver polygons (highly elongated polygons) that are unable to be evaluated further by the intersection operation. In order to bypass this computational challenge we implement instead the SR using a regular rectangular grid, calculating the list of original pixels contained fully or partially within each super-cell of the new rectangular grid. As the SR grid becomes finer, the resulting SR map becomes an oversampled version of the theoretical case. The SR strategy described above, is presented schematically in Figure 9. In Figure 9a we present the SR regular rectangular grid in red (size of super-cells is exaggerated for clarity) on top of all the different layers of the photometrically corrected global map of Enceladus at 1.80401 \(\mu \)m (binned by spatial sampling). One super-cell of the grid is marked in bold red and traverses all the different layers of the map to signify the contribution of all the overlapping pixels of different spatial sampling. An example of one super-cell of size 0.25\(^\circ \)×0.25\(^\circ \) with all its overlapping pixels corrected at 1.80401 \(\mu \)m is presented in Figure 9b. In this 3D view the overlapping pixels of all the different layers are collapsed on the same plane, with the super-cell in red. A zoom-out version of the super-cell is extruded with its radiance value being calculated as the weighted average of the values of all the overlapping pixels.

Figure 9: Schematic representation of the SR strategy with (a) an example of a SR grid in red and the different layers of the photometrically corrected global map of Enceladus at 1.80401 \(\mu \)m (binned by spatial sampling) and (b) an example of one super-cell (in red) of size 0.25\({^\circ }\times 0.25{^\circ }\) with all its overlapping pixels corrected at 1.80401 \(\mu \)m.

The question that naturally arises at this point is which grid size guarantees the best possible SR map without undersampling the highly resolved areas of the original map and at the same time it is not computationally too time-consuming. Any further increase to the grid size would not yield any better results, since we will be oversampling within the smallest element of the original map. By testing gradually finer rectangular grids in the specific case of our VIMS data set for Enceladus, we deduce that a grid of 1440\(\times \)720 cells, where each super-cell has a size of 0.25\(^\circ \)\(\times \)0.25\(^\circ \), is the optimum compromise between resolving the smallest pixels in the data set (found in cube 1637463200_1 with sampling 1.72 km/pixel and angular size of approximately 0.6\(^\circ \)\(\times \)0.2\(^\circ \)) and computational cost. Regarding the algorithmic implementation of the WASR methodology, the initial approach consists of the sequential evaluation of the intersection of the whole super-cell grid against each one of the pixels in the data set. The computational time required (on a MacBook Pro with an Apple M1 Max chip and 10 cores) increases exponentially with the size of the SR grid and for 1440\(\times \)720 cells is more than 46 hours. As it becomes clear, any SR grid finer than that (to be implemented on icy-moons data sets with smaller pixels e.g. the ones expected to be acquired by JUpiter ICy moons Explorer (JUICE) on Ganymede ; Grasset et al., 2013) becomes impractical and requires the implementation of parallel computation. It is important to note, that the WASR method we have developed is very well suited for parallelization, since the regular global grid of 360\(^\circ \)\(\times \)180\(^\circ \) can be segmented into independent panels (e.g. 9 panels of 120\(^\circ \)\(\times \)60\(^\circ \)), which can then be passed in parallel through the SR pipeline and the results consolidated in a final global product. Another way to make the calculation more time-efficient would be to implement an adaptive grid, instead of the regular grid we have selected, which will automatically increase or decrease density of the super-cells according to the size of the pixels available in each area of the map. This approach has not been investigated in the current work and can be part of future expansions of the methodology. Instead we have tested and implemented the Sort-Tile-Recursive (STR) method (Leutenegger et al., 1997), a simple and efficient algorithm for a common indexing technique for spatial data (R-tree), which decreases the computational time to a remarkable 15 minutes.

Figure 10: SR global mosaic of the VIMS observations of Enceladus at wavelength 1.80401 \(\mu \)m for the 1440x720 grid (super-cell size 0.25\(^\circ \)\(\times \)0.25\(^\circ \)) and \(n=2\) in (a) equirectangular projection and orthographic projection centered at the (b) north and (c) south pole.

The final SR cartographic products can be seen in Figure 10a, where we present the global SR map of Enceladus in equirectangular projection at 1.80401 \(\mu \)m for a regular rectangular SR grid of 1440\(\times \)720 cells and \(n=2\). Similarly we present the orthographic projections of the two maps centered at the north and south pole in Figures 10b and 10c respectively. A first visual inspection of the new Enceladus SR map at the global scale (macroscale) already demonstrates the capacity of our methodology to further improve the photometric correction implemented in the previous step, as visible from the further diminishing of visible seams between different observations. A more continuous signal is achieved in the global scale without harsh transitions/artifacts, while the natural variation due to surface features is enhanced, with large scale structures becoming more prominent, better defined and continuous between different observations. New spectral variability is revealed, as is the case for the dark (high-absorption) features south to a similar prominent feature centered at latitude 30\(^\circ \)N and longitude 90\(^\circ \)W, that has been also identified in previous studies (Robidel et al., 2020), matching with the geological sub-unit LH smooth as defined by Crow-Willard and Pappalardo (2015). We can also observe that due to the limited number of overlapping pixels at the polar regions, the SR results are not as striking as in the equatorial regions. The few seams that remain in the SR map can be traced back to imperfections in the photometric correction. Where the photometric model can not achieve good integration of neighboring observation, the SR map will also present some visible seams between those observations. These imperfections can also be associated with areas where we have a very well resolved observation on top of significantly lower sampling observation, creating a disparity between the two that can not be compensated for by the SR model.

In Figure 11 we are visualizing for the same grid size (1440\(\times \)720 cells) the number of overlapping pixels for each super-cell, providing a view of the density of overlaps at a global scale. As we can see, these values vary substantially with surface location, from 224 to 2 pixels per super-cell for the highest and lowest density. It is worth noting that the leading hemisphere has consistently a higher coverage while the lowest number of overlapping pixels can be seen, as expected from the design of Cassini’s orbital tour, at the polar regions. We should also underline that overlaps do not occur exclusively between different cubes, but due to the “whisk broom” operation mode of the VIMS-IR instrument, it is common to observe overlaps between the pixels of different rows of the same cube (“self-overlapping”).

Figure 11: Global map of number of overlapping pixels for SR grid size 1440\(\times \)720 cells (super-cell size 0.25\(^\circ \)\(\times \)0.25\(^\circ \)).

4 Validation of the WASR method

Beyond the visual improvement of the final cartographic product of Enceladus achieved by the SR methodology, as described in the previous section, we have validated our approach by visual comparison with the ISS mosaic at the mesoscale and in a more quantitative way by comparing spectra (spectral dimension) and transects (spatial dimension) between the photometrically corrected, the new SR and the ISS maps.

In Figure 12 we inspect at a smaller scale the global cartographic products presented in Figures 6 and 19 at wavelength 1.80401 \(\mu \)m for two regions: for longitude 100\(^\circ \)W to 50\(^\circ \)E and latitude 50\(^\circ \)S to 10\(^\circ \)S (left column) and for longitude 50\(^\circ \)W to 90\(^\circ \)E and latitude 40\(^\circ \)S to 0\(^\circ \) (right column). The color scale has been artificially constrained to the [0.5, 1] range for improved contrast of the surface features. Same color arrows note corresponding features between the photometrically corrected map (Figs. 12a and b), the SR for \(n=2\) map (Figs. 12c and d) and the high-resolution ISS mosaic by Bland et al. (2018) (12e and f) that become more prominent or are completely retrieved by SR from the sub-pixel information available in the different layers of the photometrically corrected map. In the first region presented (Figs. 12a, c and e), the orange arrows point to the wide trough creating the border of the geomorphological unit defined as Central Leading Hemisphere (clh) terrain by Crow-Willard and Pappalardo (2015), which although already visible in the photometrically corrected maps it becomes more prominent in the SR map. Even more striking are the results of SR in retrieving parts that are not visible in the photometrically corrected map from the wide trough at the south border of the Trailing Hemisphere curvilinear (cl1) terrain (red arrows) and the narrow trough (green arrow) at the edge of the western border of the transitional (t) unit. Lastly, we mark with the white arrow the location of the only named crater in this area (Kamar crater), which can also be studied in more detail thanks to SR. In the second region under investigation (Figs. 12b, d and f), we focus on the area around Labtyat Sulci. Here we can observe how the SR method retrieves the northern part of the sulci (cyan arrows) from the combination of the overlapping pixels (Fig. 12d), information completely obscured in the photometrically corrected map (Fig. 12b), while better constraining spatially its southern extension. In this scale (global and mesoscale) we can already observe the impressive results of SR in improving the consolidation of different observations into mosaics and in retrieving at the same time surface features at finer detail. A more thorough inspection of the resulting cartographic product in the microscale will be presented in the following paragraphs.

Figure 12: Mesoscale comparison for two regions (left and right columns) of photometrically corrected (a and b), SR for \(n=2\) (c and d) and Enceladus ISS global mosaic by Bland et al. (2018) (e and f). Colored arrows mark surface features that become more prominent or retrieved altogether by SR, while in white we mark named surface features.

We will now check the validity of the SR approach in the spectral dimension by comparing the spectra of the selected super-cell resulting from the SR calculation with that of the most resolved pixel among all the overlapping pixels used for the SR calculation. In Figure 13 we present the spectra comparison for super-cell with index 374 040 (super-cells in the map are identified by an index from 0 to 1440\(\times \)720-1 in compliance with the python numbering scheme). In Figure 13a we mark, in the photometrically corrected map, the most resolved pixel from the list of overlapping pixels in blue (depicted in the zoom-in panel for clarity). In the SR map (Fig. 13b) we mark the corresponding super-cell in orange for grid size 1440\(\times \)720 cells (depicted in the zoom-in panel for clarity). In Figure 13c we present the 189 photometrically corrected spectra of all the overlapping pixels for the selected super-cell, the linear combination of which comprises the new spectrum of the super-cell. In the comparison of the spectra of the super-cell with the most resolved of its photometrically corrected overlapping pixels (Fig. 13d) we also include the two sigma envelopes (the calculation of which is explained in Appendix E of the supporting material), which for most of the wavelengths are negligible compared to the spectral values (as it can be seen in the zoom-in panel). This comparison shows that the two spectra almost perfectly coincide, evident from the near zero residuals provided in the adjacent panel, validating that the SR calculation conserves the information in the spectral dimension. The black arrow in the bottom left figure points at part of a saturated spectrum that evaded the filtering method presented in Section 3.1, but which nonetheless does not significantly effect the final SR spectrum. As already mentioned, the SR method operates only in the spatial dimension of the hyperspectral cartographic products. However, by inspecting the comparison between the SR spectrum and the photometrically corrected spectrum of the most resolved (top most) of the VIMS overlapping pixels (Fig. 13d), we observe that the SR approach introduces indirectly two major improvements in the spectral dimension: 1) SR achieves through the weighted averaging of all the overlapping spectra the retrieval of the part of the spectrum of the most resolved pixel that is missing due to saturation in the shorter wavelengths and 2) SR achieves, again through the weighted averaging of the overlapping pixels, lower noise levels in the longer wavelengths than what would have been possible by using only the most resolved pixel to construct our map.

Figure 13: Spectra comparison of photometrically corrected map (a) with the most resolved of the overlapping pixels show in blue (depicted in the zoom-in panel for clarity) and SR map at grid size 1440\(\times \)720 cells (b) with the super-cell with index 374 040 in orange (depicted in the zoom-in panel for clarity) at wavelength 1.80401 \(\mu \)m. Bottom panels present the 189 photometrically corrected spectra of all overlapping pixels with the black arrow pointing at spectrum saturation residuals (c) and the comparison of the spectra of the most resolved pixel against the corresponding super-cell accompanied by the respective residuals (d), with a zoom-in panel added for better visibility of the two-sigma error envelope (the calculation of which is detailed in Appendix E of the supporting material).

The validation in the spatial dimension is presented in Figure 14, where we plot a transect extending from longitude 100\(^\circ \)E to 180\(^\circ \)E at latitude 1\(^\circ \)N, over the photometrically corrected (Fig. 14a) and SR for \(n=2\) (Fig. 14b) maps, in red and green respectively. In Figure 14c we present the radiance factor values over the transect for the photometrically corrected and SR (for \(n=\)1, 2, 3 and 4) maps, with the vertical dotted lines marking the transition points between different cubes. Figure 14d presents the difference between the SR radiance factor for the values of \(n\) and the corrected radiance factor, to demonstrate that the SR values approach the photometrically corrected values for increasing \(n\). In a similar fashion we plot, in Figure 15, a transect extending from longitude 100\(^\circ \)W to 50\(^\circ \)E at latitude 80\(^\circ \)S, in the highly interesting and diverse south polar area with its tiger stripe (in this case we present only the SR transect for \(n=2\) for clarity). In both cases, the transect on the SR map traces a much smoother profile between transitions compared to the transect on the photometrically corrected map manifesting as further diminishing of the visible seams between neighboring observations. Most importantly for validating the SR calculation we can see that the SR transect not only maintains surface features, but introduces details at finer spatial scale than the original VIMS pixel sampling. Hence, SR provides a global cartographic product of improved quality and greater detail at smaller scale, especially in locations where the photometrically corrected map is not well resolved.

Figure 14: Transect comparison (latitude 1\(^\circ \)N ranging from longitude 100\(^\circ \)E to 180\(^\circ \)E) of photometrically corrected map (a) and SR map for \(n=2\) at grid size 1440\(\times \)720 cells (b) at wavelength 1.80401 \(\mu \)m. Panel (c) presents the comparison of the radiance factor along the transect for the photometrically corrected map and the SR map for \(n\) = 1, 2 and 3, with the dashed vertical lines marking the transition points between cube. Panel (d) presents the difference between the SR radiance factor for the different values of \(n\) with the photometrically corrected radiance factor.

Figure 15: Transect comparison (latitude 80\(^\circ \)S ranging from longitude 100\(^\circ \)W to 50\(^\circ \)E) of photometrically corrected map (a) and SR map for \(n=2\) at grid size 1440\(\times \)720 cells (b) at wavelength 1.80401 \(\mu \)m. Panel (c) presents the comparison of the two transects with the dashed vertical lines marking the transition points between cube. Panel (d) presents the difference between the SR radiance factor with the photometrically corrected radiance factor.

As we have demonstrated above from the transect comparison, the SR method introduces new features (manifesting as variation in the radiance factor value) at a finer scale than the original VIMS pixel size. The question that remains to be addressed to complete the validation of the methodology is whether these variations can be associated with real features and are not artifacts resulting from the combination of the signal from overlapping pixels. The validation strategy employed involves the recalculation of the photometrically corrected and SR resolution maps of Enceladus (as described in Sections 3.2 and 3.3 respectively), after removing a selected HR cube from the data set, namely cube 1487302144_1, to be used separately as the HR reference image. If the SR map without the HR reference manages to reproduce small scale surface features available in the HR reference and not in the photometrically corrected map without the HR reference, then the validation test is deemed successful. In Figure 16 we present regional maps of the area covered by the HR reference cube 1487302144_1 at standard wavelength 1.80401 \(\mu \)m: the photometrically corrected without the HR reference cube (top left), the 1440×720 super-cells SR map resulting from the photometrically corrected map without the HR reference for \(n=2\) (top middle) and the photometrically corrected HR reference cube (top right). Each map includes a transect at latitude 1\(^\circ \)N ranging from longitude 79\(^\circ \)E to 85\(^\circ \)E, the radiance factor profile of which is presented in the middle panel (red for corrected without the HR reference, green for SR without the HR reference and purple for HR reference), with the vertical dotted lines marking the borders between cubes. Because the small scale features retrieved by the SR approach can be subtle, we present in the bottom panel the absolute gradient values of the transect profiles (red and purple lines shifted vertically by +0.4 and –0.4 respectively for clarity), in order to reveal the location of feature borders by the sharp peaks (peaks aligning with the dotted vertical lines representing cube borders are naturally interpreted as seams-artifacts, artificially enhanced here by the gradient operation). In the bottom panel we offer the validation of the SR approach as we can observe that peaks in purple (VIMS HR reference) are associated with peaks in green (SR without the reference) in locations where the red line is flat (no feature in the corrected map without the HR reference) as indicated by the black arrows. Naturally, the green peaks are not as prominent as their corresponding purple peaks, since the SR method can retrieve the HR features but of course not with the same contrast-clarity as a HR observation. The slight shift of the location of the corresponding peaks can be attributed to the fact that HR features are not retrieved in the same resolution as the HR acquisition and therefore their borders identified by the transect peaks present at slightly different longitude. We can also observe that SR further improves the results of the photometric correction, in the sense that the sharp peaks associated with seams (e.g. between longitude 79\(^\circ \)E and 80\(^\circ \)E) further diminish.

Figure 16: Regional maps of Enceladus at 1.80401 \(\mu \)m encompassing HR reference cube 1487302144_1: photometrically corrected without the HR reference cube (a), 1440\(\times \)720 super-cells SR map resulting from the photometrically corrected map without the HR reference for \(n=2\) (b) and the photometrically corrected HR reference cube (c). Transects at latitude 1\(^\circ \) ranging from longitude 79\(^\circ \)E to 85\(^\circ \)E in red, green and purple respectively, with the vertical lines marking the segment of the transect presented in panel (d) and the vertical dotted lines marking the borders between cubes. Panel (e) presents the absolute gradient values of the transects (red and purple lines shifted vertically by +0.4 and –0.4 respectively for clarity), with arrows pointing to HR features retrieved in the SR map.

Similarly, in Figure 17 we validate the results of the SR methodology by comparing the transects at latitude 1\(^\circ \)N extending from longitudes 72\(^\circ \)E to 90\(^\circ \)E between the photometrically corrected (Fig. 17a), SR (Fig. 17b) and ISS maps (Fig. 17c) by Bland et al. (2018). Since we are solely interested in feature localization between maps of different signals and different treatment pipelines, we only present the gradient of the three transects (Fig. 17d) for the photometrically corrected map (red line, shifted vertically -1 for clarity), the SR map (green line) and ISS mosaic (purple line, divided by 1000 to scale it to the level of the VIMS signal transects and shifted vertically +0.5 for clarity). Similarly to Figure 16, where we compare within the same data set by removing a HR cube, the comparison with the ISS mosaic also reveals areas where the SR method retrieves features (peaks in the green line) associated with peaks in the ISS transect (purple line) that are not present in the photometrically corrected map (red line). Two such cases are highlighted by the dashed black boxes.

Figure 17: Regional maps of Enceladus at 1.80401 \(\mu \)m encompassing HR reference cube 1487302144_1: photometrically corrected (a), 1440\(\times \)720 super-cells SR map for \(n=2\) (b) and the ISS mosaic by Bland et al. (2018) (c) with transects at latitude 1\(^\circ \)N ranging from longitude 72\(^\circ \)E to 90\(^\circ \)E in red, green and purple respectively. Panel (d) presents the absolute gradient values of the transects (red and purple lines shifted by –1 and +0.5 respectively for clarity), with dashed boxes pointing to areas where features in the ISS map are retrieved in the SR map. The ISS signal has been divided by 1000 for clarity.

Finally, we define a new metric for the spatial sampling of the SR map, henceforth called “effective sampling” (\(s_{eff}\)), and through that the gain in spatial sampling, with regards to the original sampling of the initial map (\(s_{VIMS}\)), defined as:

\begin{equation} gain = \frac {s_{VIMS}}{s_{eff}}.\tag {12} \end{equation}

The method for calculating the effective sampling for each super-cell is detailed in the supporting material (Appendix D). In Figure 18 we present the global map of Enceladus at grid size 1440×720 cells with the color-map corresponding to the values of the gain in spatial sampling achieved by SR compared to the original VIMS sampling. From these results we showcase that 13% of the super-cells have gain up to twice the original VIMS sampling, while 32% of the super-cells have gain greater or equal to 5. Thus, we demonstrate that the majority of the super-cells in the SR map achieve a good improvement in the original VIMS sampling even in the conservative approach we have selected for the calculation of the effective sampling of the SR map. A question that naturally arises is that since there are super-cells with gains as high as 20 (20 times better sampling than the native VIMS sampling in that area), we should be able to resolve features at very small scale using SR. However, under closer examination of the global map (Fig. 18), we understand that these super-cells are very limited in number and generally isolated, surrounded by super-cell with significantly smaller gain. Beyond the extreme cases, the higher values in gain are associated with areas where the original VIMS sampling values are relatively low. Contrary, areas over well resolved cubes exhibit gain values closer to one. In this context, the gain values calculated are consistent with the visible spatial improvement between the photometrically corrected and the SR global maps.

Figure 18: Global map of Enceladus of gain in spatial sampling achieved by SR at grid size 1440×720 cells.

5 Spectral diversity of the new SR map of Enceladus

The main scope of this paper has been to present in detail the methodology we have developed for SR mapping of icy moons in order to provide improved cartographic and spectroscopic products with finer detail at smaller spatial scales. This methodology has been developed as a tool to offer the opportunity to reanalyze data sets from previous and current missions, optimize the planning of future observations, and ultimately help gain new insight in the connection between geomorphology, composition and physical state of icy moons surface. This information is crucial in constraining important scientific questions regarding their internal structure, evolution and potential habitability.

In order to investigate the spectral diversity and the variations in the physical state of Enceladus surface regolith with our new SR mapping technique, we create a new false color composite SR maps of Enceladus with different characteristic wavelengths at each one of the different color channels Red (R), Green (G), Blue (B), similarly to Robidel et al. (2020): R: 3.09689 \(\mu \)m/ 1.65567 \(\mu \)m, G: 2.00167 \(\mu \)m, B: 1.80401 \(\mu \)m. In Figure 19a we show the global RGB composite map of Enceladus in equirectangular projection for a SR grid of 1440\(\times \)720 cells and in orthographic projection centered at the north and south poles in Figures 19b and 19c respectively. The blue channel is normalized for the maximum and minimum values in the corresponding wavelength, while the green and red channels are normalized between their minimum values and 0.4, in order to better portray their dynamic range.

The selection of the wavelengths for each color channel in the RGB composite is based on the characteristics of the water-ice spectrum, which largely dominates the surface of Enceladus. The water-ice spectrum is characterized by absorption bands at 1.25, 1.5, 1.65 and 2 \(\mu \)m, with their existence, position, shape, intensity and width being indicative of the structure of the water-ice regolith as well as its temperature and grain size (Brown et al., 2006; Jaumann et al., 2008; Newman et al., 2008; Robidel et al., 2020; Scipioni et al., 2017). As a general rule the depth of the absorption bands increases with water-ice abundance and particle size (Filacchione et al., 2012; Jaumann et al., 2008). Finally, a strong indicator of the crystallinity of the water-ice compared to amorphous ice are the 1.65 \(\mu \)m absorption band and the so called “Fresnel” peak at 3.1 \(\mu \)m, that are more prominent for crystalline ice (Newman et al., 2008). The degree of crystallinity is an important characteristic, since it can be associated with younger areas or plume depositions if we consider that crystalline ice converts to amorphous ice under the influence of high-energy particle bombardment (Brown et al., 2006).

Figure 19: Composite RGB SR global mosaic of the VIMS observations of Enceladus for the 1440×720 grid (super-cell size 0.25\(^\circ \)\(\times \)0.25\(^\circ \)) and \(n=2\) in (a) equirectangular projection and orthographic projection centered at the (b) north and (c) south pole. Color channels are controlled by R: 3.09689 \(\mu \)m / 1.65567 \(\mu \)m, G: 2.00167 \(\mu \)m and B: 1.80401 \(\mu \)m.

As we can observe the composite map is dominated by the green-blue channels (controlled by wavelengths 2.00167 \(\mu \)m and 1.80401 \(\mu \)m), an expected result considering that Enceladus is primarily composed of water-ice. The red channel (controlled by the 3.09689 \(\mu \)m/ 1.65567 \(\mu \)m ratio) is very prominent at the location of the fractures close to the south pole, indicating crystalline water-ice. Lighter hues of green-blue are overlaid on the cratered terrain, whereas stronger reddish tones cover in a spatially contained way the tiger stripes. Similarly other fractures on the surface of Enceladus and large regions among the youngest terrains bear the characteristic reddish hue also visible in the area surrounding the south polar terrains and somewhat centered at leading and trailing hemispheres. In particular, the reddish prominent feature in the maps is a high absorption area, both in the G (2.00167 \(\mu \)m) and B (1.80401 \(\mu \)m) channels, centered at latitude 30\(^\circ \)N and longitude 90\(^\circ \)W, that has been also identified in previous studies (Robidel et al., 2020), matching with the geological sub-unit Leading Hemisphere (LH) smooth as defined by Crow-Willard and Pappalardo (2015). For the same channels we can also observe higher \(\frac {I}{F}\) values (lower absorption) in the Cratered Terrains dividing at the sub- and anti-Saturn hemispheres, as well as north of the Leading Hemisphere Terrain (LHT) and Trailing Hemisphere Terrain (THT) units and in pockets at the border between the THT and the South Polar Terrain (SPT) (units as defined by Crow-Willard & Pappalardo, 2015).

In Figure 20 we present the same color composite map overlaid on the ISS global mosaic (Bland et al., 2018), allowing for a direct association of the different colors of spectral information with the highly resolved surface geomorphology of the ISS map. In this manner we can observe how the higher radiance factor values in the R (3.09689 \(\mu \)m/ 1.65567 \(\mu \)m) channel, associated with crystalline water-ice, cover the tiger stripes, where in white we mark the location of plumes identified by Kempf et al. (2010) (I, II, III, etc.), the area surrounding the SPT and then spread to higher latitudes coating the leading and trailing hemispheres all the way to the north polar area. In the same figures we can observe the prominent comparatively low absorption of the cratered terrain between the leading and the trailing hemispheres, consistent with an older in age terrain as justified by the existence of the numerous impact craters.

Figure 20: Composite RGB SR map of Enceladus at grid size 1440\(\times \)720 cells, overlaid over the ISS global mosaic by Bland et al. (2018). Color channels are controlled by R: 3.09689 \(\mu \)m/ 1.65567 \(\mu \)m, G: 2.00167 \(\mu \)m and B: 1.80401 \(\mu \)m. In white we mark five geomorphological sub-units (ls-LH smooth, clh-central LH, cl3-southern curvilinear, nl-northern lineated and part of the cl1-TH curvilear) from Crow-Willard and Pappalardo (2015). Numbered points in the south polar terrain indicate the location of plumes identified by Kempf et al. (2010) (I, II, III, etc.).

In Figure 20 we moreover mark in dashed white lines five units from Crow-Willard and Pappalardo (2015): LH smooth (ls), central LH (clh) southern curvilinear (cl3), northern lineated (nl) and part of the TH curvilinear (cl1). The unit ls is visibly associated with a prominent spectral feature centered at latitude 30\(^\circ \)N and longitude 90\(^\circ \)W. This feature is defined by a combination of higher absorption in the green (2.00167 \(\mu \)m) and blue (1.80401 \(\mu \)m) channels and a stronger signal in the red channel (3.09689 \(\mu \)m/ 1.65567 \(\mu \)m). These indications characterize a potentially younger surface with larger ice grains and in crystalline form. Similar in appearance on the composite map, we identify new sub-units within the geomorphological unit clh that can be observed in more detail thanks to the new SR map. The location of these sub-units (redder hues in the composite map) is consistent with the location of wide trough presented by (Crow-Willard & Pappalardo, 2015).

Moreover, we examine the relationship between the geomorphological unit cl3 of Crow-Willard and Pappalardo (2015), that surrounds the south polar terrain, and the information of the SR composite map. As we can observe in Figure 20, the western border (between longitudes 180\(^\circ \)W and 50\(^\circ \)W) of the geomorphological unit cl3 coincides almost perfectly with the fractured area encircling the south pole exhibiting a strong signal in the red channel (3.09689 \(\mu \)m/ 1.65567 \(\mu \)m) indicative of crystalline water-ice. In contrast and as we follow the border of geomorphological unit cl3 to the east, we can identify intrusions that appear spectrally similar to the older inter-hemisphere cratered terrain, which overlay areas shown in the ISS background map to be significantly less fractured than the western parts of the cl3 unit that encloses the south pole and are highlighted with dark red hues. Finally, the units cl1 and nl marked in Figure 20 enclose the spectrally distinct corridors separating the cratered plains to the west from the striated planes to the east and two crated plains in the north (distinctly brighter) respectively.

6 Discussion

In the preceding sections we have described our novel SR methodology for icy moons by sequentially describing the method, validating it and presenting the results of each one of its three main steps: 1) Navigation, 2) Photometric correction and 3) SR. In this section we will discuss the improvements and choices made compared to existing methodology and explain its limitations.

Following the well-established paradigm for the format of the VIMS navigation cubes (Le Mouélic et al., 2019b; Penteado & Trilling, 2015), we have created our own version of the navigation cubes introducing four important improvements in the navigation code: 1) we develop the code in Python aligning with the state-of-the-art in scientific coding and allowing for the exploitation of the abundance of online resources and open-source packages, 2) each pixel is treated as a polygon with a predefined number of vertices and not as a rectangle (as is the common approach) for increased accuracy in the projection, a characteristic necessary for the accurate merging of the overlapping pixels, 3) we include a DSM for Enceladus for increased accuracy in the calculation of the navigation and 4) we introduce a new more flexible metric for calculating the spatial sampling of the pixels. Despite these improvements, the accuracy of the navigation remains bounded by the accuracy of the SPICE kernels employed for the calculation of the location of the pixels on the surface of the target. Nonetheless, the inspection of the global mosaics of all steps (uncorrected, photometrically corrected and SR) presents a good continuity of the surface features among different observations, suggesting that any potential navigation error is smaller than the VIMS resolving capacity. As noted before, an accurate navigation is an integral part of the SR methodology as it ensures that we are always mixing pixels corresponding to the same areas. As the spatial sampling of the observation improves (as is expected to be the case for the observations of JUICE), so does the need for increased accuracy in the navigation of the observations (e.g. more accurate kernels).

Even though we have implemented an existing methodology for the photometric correction of the Enceladus global mosaic, we have further expanded on the pivotal work of Robidel et al. (2020) in the following ways: 1) evaluation of the whole Enceladus data set available at the Nantes VIMS Data Portal and not just 355 selected cubes, allowing for increased statistics, 2) filtering of the data to be used at the pixel level, 3) holistic treatment of the saturation issue, 4) treatment of pixels as polygons with a predetermined number of vertices instead of rectangles, 5) addition of a DSM for Enceladus, instead of the triaxial ellipsoid, allowing for more accurate values for the photometric angles used in the photometric model, 6) calculation of the linear coefficients for all the bands of the instrument, and 7) calculation of radiometric error for each photometrically corrected \(\frac {I}{F}\) value. Similarly to the navigation step, an effective photometric correction of the overlapping observations is critical for the implementation of SR, ensuring that we are combining pixels under the same observation and illumination conditions. The limitations of this step are dictated by the photometric model which offers a statistical treatment of the photometric behavior at the global scale. This approach works sufficiently well for a relatively homogeneous surface as is the case for Enceladus, but might require segmentation or complete substitution by a radiative-transfer based model for different targets. The limitations of the photometric correction model are manifesting in the photometrically corrected map as a few remaining visible seams between neighboring observations.

In the third and final step we have demonstrated how the combination of all the different layers of overlapping pixels can produce the SR global mosaic with further diminished seams and finer detail at a smaller scale that is not available in the photometrically corrected map and corresponds to real features. Isolating the SR step from the limitations propagated by the previous steps, a first order constraint appears to be the number of overlapping pixels associated with each super-cell and therefore the number of observations over the same area. A first attempt to quantify the effectiveness of SR has been made by defining the effective sampling, which as explained in Appendix D can have certain limiting assumptions. At the same time, the distribution of the spatial sampling and the radiance factor of the overlapping pixels can potentially have an effect on the identification of the areas of the map where SR works better compared to others. This multi-variable dependency of the effectiveness of SR requires further investigation, that can not only contribute to further understand the limitations of implementing our SR methodology on existing data sets, but can also help define observation strategies for maximizing the scientific returns of future missions such as JUICE.

7 Conclusion

In this work, we have developed and validated a complete and modular pipeline for producing SR cartographic products of icy moons, with a specific application to Cassini-VIMS observations of Enceladus. Our methodology integrates three key steps: 1) accurate pixel-level navigation, 2) global photometric correction, and 3) a novel Weighted Average Super-Resolution (WASR) approach, designed to address the unique challenges posed by heterogeneous, hyperspectral data sets acquired during multiple flybys.

By treating VIMS pixels as polygons rather than idealized rectangles and by incorporating a digital shape model of Enceladus, we significantly improve the geometric accuracy of the projected data. The global photometric correction, based on the Kaasalainen–Shkuratov formalism, effectively mitigates illumination-related artifacts and produces homogeneous mosaics suitable for subsequent super-resolution processing. The WASR method then exploits the sub-pixel spatial offsets among overlapping observations to construct global maps with enhanced spatial detail without introducing artificial spectral information. Finally, we validate the WASR approach through multiple complementary tests, including spectral conservation, transect analyses, a high-resolution reference cube and high-resolution ISS mosaic comparison. The effective sampling gain analysis further quantifies the spatial improvement achieved, showing that a substantial fraction of the surface benefits from significant spatial sampling, hence resolution, enhancement.

The resulting SR maps of Enceladus reveal improved continuity across observations, reduced photometric seams, and enhanced visibility of geological features, including tectonic structures and terrain boundaries previously obscured by resolution limitations. These hyperspectral cartographic products are now ready to be used in combination with spectral analysis models to retrieve compositional and physical properties of the surface of Enceladus at finer spatial scale than anticipated with the current Cassini VIMS dataset.

The pipeline presented here is fully transferable to other icy moons and to future hyperspectral data sets, including those expected from the JUICE mission. More broadly, this work demonstrates that meaningful SR mapping can be achieved for planetary bodies with limited and heterogeneous observational coverage, providing a powerful tool for maximizing the scientific return of existing and future missions.

Acknowledgements. This work is the output of the corresponding PhD project “Development of super-resolution mapping methods for icy moons: Application to Enceladus, in preparation for future observations by the ESA JUICE mission” (Ntinos, 2024), carried out in the Planétologie et Sciences Spatiales (PSS) research group of the Institut de Physique du Globe de Paris (IPGP), part of the Université Paris-Cité (UPC) and has been co-funded by the European Space Agency (ESA) in the context of the Open Space Innovation Platform (OSIP), project ID: I-2020-00336. Authors also acknowledge funding support from the French National Research Agency (Project RaD3-net: ANR-21-CE49-0020-567 04/RAD3-NET and EOLE: ANR-23-CE56-0008) and from the CNES (Centre national d’études spatiales).

Open science statements

Author contributions. Conception: SR, SLM, TC, NA, GT, CV. Data Preparation: CN, SR, SLM, RR, BS. Methodology: CN, SR. Software: CN. Validation: CN, SR. Formal analysis: CN, SR. Investigation: CN, SR. Writing-original draft: CN. Writing-review and editing: CN, SR, SLM, RR, BS, TC. Funding acquisition: SR, NA. All authors contributed to the article and approved the submitted version. Code availability. All code is available on Gitlab (https://pss-gitlab.math.univ-paris-diderot.fr/ntinos/super-resolution-for-icy-moons). In the same repository we provide the list of the VIMS cubes used in this study (_Enceladus_cubes_ds4.txt) and the list of metakernels containing the SPICE kernels used (most updated at the time of the development). All versions of the code can be cited by using the DOI 10.5281/zenodo.21037587. Data availability. The two super-resolution maps presented in this work (monochromatic at 1.80401\(\mu {}\)m and RGB composite) are available on https://zenodo.org/records/23669861. The calibrated VIMS cubes (.cub) used in this study can be retrieved from the Nantes Cassini VIMS Data Portal by the Laboratoire de Planétologie et Géosciences (LPG) of Nantes, through its publicly available interface (https://vims.univ-nantes.fr/). The five additional cubes used can be publicly accessed from the Planetary Data System (PDS) cartography and imaging sciences node (https://pds-imaging.jpl.nasa.gov/volumes/iss.html) and calibrated using the third version of the Integrated Software for Imagers and Spectrometers (USGS-ISIS3) found in https://isis.astrogeology.usgs.gov/8.3.0/. Finally the DSM that results from the photogrammetric inversion can be found in the Annex of the PDS Cartography and Imaging Sciences Node, provided by the USGS in support of the public delivery of NASA-funded geospatial data (provided in jpg and csv format in https://astrogeology.usgs.gov/search/map/Enceladus/enceladus_cassini_iss_shapemodel_bland_2019/enceladus_2019pm_point_cloud). Competing interests. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. BS is the treasurer and a board member of Planetary Research Cooperative but was not involved in the peer-review nor in the editorial decisions for this publication.

References

Acton, C. H. (1996). Ancillary data services of NASA’s navigation and ancillary information facility. Planetary and Space Science, 44(1), 65–70. https://doi.org/10.1016/0032-0633(95)00107-7

Alonso-Fernandez, F., Farrugia, R. A., Bigun, J., Fierrez, J., & Gonzalez-Sosa, E. (2018). A survey of super-resolution in iris biometrics with evaluation of dictionary-learning. IEEE Access, 7, 6519–6544. https://ieeexplore.ieee.org/abstract/document/8586871

Archinal, B. A., Acton, C., A’hearn, M., Conrad, A., Consolmagno, G., Duxbury, T., Hestroffer, D., Hilton, J., Kirk, R. L., Klioner, S., et al. (2018). Report of the IAU working group on cartographic coordinates and rotational elements: 2015. Celestial Mechanics and Dynamical Astronomy, 130(3), 22. https://doi.org/10.1007/s10569-017-9805-5

Bland, M. T., Becker, T. L., Edmundson, K. L., Roatsch, T., Archinal, B. A., Takir, D., Patterson, G. W., Collins, G. C., Schenk, P. M., Pappalardo, R. T., et al. (2018). A new Enceladus global control network, image mosaic, and updated pointing kernels from Cassini’s 13-year mission. Earth and Space Science, 5(10), 604–621.

Bland, M. T., Weller, L. A., Mayer, D., & Archinal, B. (2020). A global shape model for Saturn’s moon Enceladus from a dense photogrammetric control network. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 3, 579–586. https://doi.org/10.5194/isprs-annals-V-3-2020-579-2020

Brossier, J. F., Rodriguez, S., Cornet, T., Lucas, A., Radebaugh, J., Maltagliati, L., Le Mouélic, S., Solomonidou, A., Coustenis, A., Hirtzig, M., et al. (2018). Geological evolution of Titan’s equatorial regions: Possible nature and origin of the dune material. Journal of Geophysical Research: Planets, 123(5), 1089–1112. https://doi.org/10.1029/2017JE005399

Brown, R., Baines, K., Bellucci, G., Bibring, J.-P., Buratti, B., Capaccioni, F., Cerroni, P., Clark, R., Coradini, A., Cruikshank, D., et al. (2004). The cassini visual and infrared mapping spectrometer (vims) investigation. Space Science Reviews, 115(1), 111–168.

Brown, R. H., Clark, R. N., Buratti, B. J., Cruikshank, D. P., Barnes, J. W., Mastrapa, R. M., Bauer, J., Newman, S., Momary, T., Baines, K. H., et al. (2006). Composition and physical properties of Enceladus’ surface. Science, 311(5766), 1425–1428. https://www.science.org/doi/10.1126/science.1121031

Brown, R. H., & Soderblom, L. (n.d.). VIMS IR pixel timing (Report). PDS Cartography and Imaging Sciences Node. https://pds-atmospheres.nmsu.edu/data_and_services/atmospheres_data/Cassini/logs/VIMS%20IR%20Pixel%20Timing_final-a.pdf

Capel, D., & Zisserman, A. (2000). Super-resolution enhancement of text image sequences. Proceedings 15th International Conference on Pattern Recognition. ICPR-2000, 1, 600–605. https://ieeexplore.ieee.org/document/905409

Clark, R. N., Brown, R. H., Lytle, D. M., & Hedman, M. (2018). The VIMS wavelength and radiometric calibration 19, final report. NASA Planetary Data System.

Clark, R. N., Curchin, J. M., Barnes, J. W., Jaumann, R., Soderblom, L., Cruikshank, D. P., Brown, R. H., Rodriguez, S., Lunine, J., Stephan, K., et al. (2010). Detection and mapping of hydrocarbon deposits on Titan. Journal of Geophysical Research: Planets, 115(E10). https://doi.org/10.1029/2009JE003369

Coutelier, M., Cordier, D., Seignovert, B., Rannou, P., Le Gall, A., Cours, T., Maltagliati, L., & Rodriguez, S. (2021). Distribution and intensity of water ice signature in south Xanadu and Tui Regio. Icarus, 364, 114464. https://doi.org/10.1016/j.icarus.2021.114464

Crow-Willard, E. N., & Pappalardo, R. T. (2015). Structural mapping of Enceladus and implications for formation of tectonized regions. Journal of Geophysical Research: Planets, 120(5), 928–950. https://doi.org/10.1002/2015JE004818

Dalle Ore, C. M., Long, C. J., Nichols-Fleming, F., Scipioni, F., Valentín, E. G. R., Oquendo, A. J. L., & Cruikshank, D. P. (2021). Dione’s wispy terrain: A cryovolcanic story? The Planetary Science Journal, 2(2), 83. https://doi.org/10.3847/PSJ/abe7ec

Es-Sayeh, M., Rodriguez, S., Coutelier, M., Rannou, P., Bézard, B., Maltagliati, L., Cornet, T., Grieger, B., Karkoschka, E., Le Mouélic, S., et al. (2023). Updated radiative transfer model for Titan in the near-infrared wavelength range: Validation against Huygens atmospheric and surface measurements and application to the Cassini/VIMS observations of the Dragonfly landing area. The Planetary Science Journal, 4(3), 44. https://iopscience.iop.org/article/10.3847/PSJ/acbd37/meta

Filacchione, G., Capaccioni, F., Clark, R., Cuzzi, J., Cruikshank, D., Coradini, A., Cerroni, P., Nicholson, P., McCord, T., Brown, R., et al. (2010). Saturn’s icy satellites investigated by Cassini–VIMS: II. results at the end of nominal mission. Icarus, 206(2), 507–523. https://doi.org/10.1016/j.icarus.2009.11.006

Filacchione, G., Ciarniello, M., d’Aversa, E., Capaccioni, F., Cerroni, P., Buratti, B., Clark, R., Stephan, K., & Plainaki, C. (2018a). Photometric modeling and VIS-IR albedo maps of Tethys from Cassini-VIMS. Geophysical Research Letters, 45(13), 6400–6407. https://doi.org/10.1029/2018GL078602

Filacchione, G., Ciarniello, M., d’Aversa, E., Capaccioni, F., Cerroni, P., Buratti, B., Clark, R., Stephan, K., & Plainaki, C. (2018b). Photometric modeling and VIS-IR albedo maps of Dione from Cassini-VIMS. Geophysical Research Letters, 45(5), 2184–2192. https://doi.org/10.1002/2017GL076869

Filacchione, G., Capaccioni, F., Ciarniello, M., Clark, R. N., Cuzzi, J. N., Nicholson, P. D., Cruikshank, D. P., Hedman, M. M., Buratti, B. J., Lunine, J., et al. (2012). Saturn’s icy satellites and rings investigated by Cassini–VIMS: III–radial compositional variability. Icarus, 220(2), 1064–1096. https://doi.org/10.1016/j.icarus.2012.06.040

Filacchione, G., Ciarniello, M., D’Aversa, E., Capaccioni, F., Clark, R. N., Buratti, B. J., Helfenstein, P., Stephan, K., & Plainaki, C. (2022). Saturn’s icy satellites investigated by Cassini-VIMS. V. spectrophotometry. Icarus, 375, 114803. https://doi.org/10.1016/j.icarus.2021.114803

Filacchione, G., D’Aversa, E., Capaccioni, F., Clark, R. N., Cruikshank, D. P., Ciarniello, M., Cerroni, P., Bellucci, G., Brown, R. H., Buratti, B. J., et al. (2016). Saturn’s icy satellites investigated by Cassini-VIMS. IV. daytime temperature maps. Icarus, 271, 292–313. https://doi.org/10.1016/j.icarus.2016.02.019

Grasset, O., Dougherty, M., Coustenis, A., Bunce, E., Erd, C., Titov, D., Blanc, M., Coates, A., Drossart, P., Fletcher, L., et al. (2013). JUpiter ICy moons Explorer (JUICE): An ESA mission to orbit Ganymede and to characterise the Jupiter system. Planetary and Space Science, 78, 1–21. https://doi.org/10.1016/j.pss.2012.12.002

Griffith, C. A., Penteado, P. F., Turner, J. D., Neish, C. D., Mitri, G., Montiel, N. J., Schoenfeld, A., & Lopes, R. M. (2019). A corridor of exposed ice-rich bedrock across Titan’s tropical region. Nature Astronomy, 3(7), 642–648. https://doi.org/10.1038/s41550-019-0756-5

Hansen, C., Shemansky, D. E., Esposito, L. W., Stewart, A., Lewis, B., Colwell, J., Hendrix, A., West, R. A., Waite Jr, J., Teolis, B., et al. (2011). The composition and structure of the Enceladus plume. Geophysical Research Letters, 38(11). https://doi.org/10.1029/2011GL047415

Hapke, B. (1981). Bidirectional reflectance spectroscopy: 1. Theory. Journal of Geophysical Research: Solid Earth, 86(B4), 3039–3054. https://doi.org/10.1029/JB086iB04p03039

Hapke, B. (1984). Bidirectional reflectance spectroscopy: 3. Correction for macroscopic roughness. Icarus, 59(1), 41–59. https://doi.org/10.1016/0019-1035(84)90054-X

Hapke, B. (1986). Bidirectional reflectance spectroscopy: 4. The extinction coefficient and the opposition effect. Icarus, 67(2), 264–280. https://doi.org/10.1016/0019-1035(86)90108-9

Hapke, B. (2002). Bidirectional reflectance spectroscopy: 5. The coherent backscatter opposition effect and anisotropic scattering. Icarus, 157(2), 523–534. https://doi.org/10.1006/icar.2002.6853

Hapke, B. (2008). Bidirectional reflectance spectroscopy: 6. Effects of porosity. Icarus, 195(2), 918–926. https://doi.org/10.1016/j.icarus.2008.01.003

Hapke, B. (2012). Bidirectional reflectance spectroscopy 7: The single particle phase function hockey stick relation. Icarus, 221(2), 1079–1083. https://doi.org/10.1016/j.icarus.2012.10.022

Hapke, B., & Wells, E. (1981). Bidirectional reflectance spectroscopy: 2. Experiments and observations. Journal of Geophysical Research: Solid Earth, 86(B4), 3055–3060. https://doi.org/10.1029/JB086iB04p03055

Jaumann, R., Stephan, K., Hansen, G., Clark, R., Buratti, B., Brown, R., Baines, K., Newman, S., Bellucci, G., Filacchione, G., et al. (2008). Distribution of icy particles across Enceladus’ surface as derived from Cassini-VIMS measurements. Icarus, 193(2), 407–419. https://doi.org/10.1016/j.icarus.2007.09.013

Kaidash, V., Gerasimenko, S. Y., Shkuratov, Y. G., Opanasenko, N., Velikodskii, Y. I., & Korokhin, V. (2009). Photometric function variations observed on the near side of the Moon: Mapping. Solar System Research, 43, 89–99. https://doi.org/10.1134/S0038094609020014

Karkoschka, E., McEwen, A., & Perry, J. (2017). Producing the best global mosaic of Titan’s surface albedo using Cassini images. 48th Annual Lunar and Planetary Science Conference, (1964), 2518. https://ui.adsabs.harvard.edu/abs/2017DPS....4930106K/abstract

Karkoschka, E. (2015). Uranus’ southern circulation revealed by Voyager 2: Unique characteristics. Icarus, 250, 294–307. https://doi.org/10.1016/j.icarus.2014.12.003

Kempf, S., Beckmann, U., & Schmidt, J. (2010). How the Enceladus dust plume feeds Saturn’s E ring. Icarus, 206(2), 446–457. https://doi.org/10.1016/j.icarus.2009.09.016

Kennedy, J. A., Israel, O., Frenkel, A., Bar-Shalom, R., & Azhari, H. (2006). Super-resolution in PET imaging. IEEE transactions on medical imaging, 25(2), 137–147. https://ieeexplore.ieee.org/document/1583760

Ku, H. H. (1969). Notes on the use of propagation of error formulas (use of propagation of error formulas). 1969., 331–341.

Le Mouélic, S., Cornet, T., Rodriguez, S., Sotin, C., Seignovert, B., Barnes, J. W., Brown, R. H., Baines, K. H., Buratti, B. J., Clark, R. N., et al. (2019a). The Cassini VIMS archive of Titan: From browse products to global infrared color maps. Icarus, 319, 121–132. https://doi.org/10.1016/j.icarus.2018.09.017

Le Mouélic, S., Seignovert, B., Rodriguez, S., Cornet, T., Sotin, C., Barnes, J., Brown, R., & Clark, R. (2019b). The Nantes’ Cassini VIMS data portal for Titan and Saturn’s icy satellites. EPSC-DPS Joint Meeting 2019, 2019, EPSC–DPS2019.

Leutenegger, S. T., Lopez, M. A., & Edgington, J. (1997). STR: A simple and efficient algorithm for R-tree packing. Proceedings 13th international conference on data engineering, 497–506. https://ieeexplore.ieee.org/abstract/document/582015

Matson, D. L., Spilker, L. J., & Lebreton, J.-P. (2002). The Cassini/Huygens mission to the Saturnian system. Space Science Reviews, 104(1), 1–58. https://doi.org/10.1023/A:1023609211620

Mccord, T. B., Coradini, A., Hibbitts, C., Capaccioni, F., Hansen, G. B., Filacchione, G., Clark, R. N., Cerroni, P., Brown, R. H., Baines, K. H., et al. (2004). Cassini VIMS observations of the Galilean satellites including the VIMS calibration procedure. Icarus, 172(1), 104–126. https://doi.org/10.1016/j.icarus.2004.07.001

McCord, T. B., Hayne, P., Combe, J.-P., Hansen, G. B., Barnes, J. W., Rodriguez, S., Le Mouélic, S., Baines, E. K. H., Buratti, B. J., Sotin, C., et al. (2008). Titan’s surface: Search for spectral diversity and composition using the Cassini VIMS investigation. Icarus, 194(1), 212–242. https://doi.org/10.1016/j.icarus.2007.08.039

Mermy, G. C., Schmidt, F., Andrieu, F., Cornet, T., Belgacem, I., & Altobelli, N. (2023). Selection of chemical species for Europa’s surface using Galileo/NIMS. Icarus, 394, 115379. https://doi.org/10.1016/j.icarus.2022.115379

Nasrollahi, K., & Moeslund, T. B. (2014). Super-resolution: A comprehensive survey. Machine vision and applications, 25, 1423–1468. https://doi.org/10.1007/s00138-014-0623-4

National Academies of Sciences, E., & Medicine. (2022). Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032 (2022) (Report). Committee on the Planetary Science, Astrobiology Decadal Survey; Space Studies Board; Division on Engineering, Physical Sciences; National Academies of Sciences, Engineering, and Medicine. https://doi.org/10.17226/26522

Newman, S. F., Buratti, B., Brown, R., Jaumann, R., Bauer, J., & Momary, T. (2008). Photometric and spectral analysis of the distribution of crystalline and amorphous ices on Enceladus as seen by Cassini. Icarus, 193(2), 397–406. https://doi.org/10.1016/j.icarus.2007.04.019

Nguyen, K., Fookes, C., Sridharan, S., Tistarelli, M., & Nixon, M. (2018). Super-resolution for biometrics: A comprehensive survey. Pattern Recognition, 78, 23–42. https://doi.org/10.1016/j.patcog.2018.01.002

Nicholson, P. D., Ansty, T., Hedman, M. M., Creel, D., Ahlers, J., Harbison, R. A., Brown, R. H., Clark, R. N., Baines, K. H., Buratti, B. J., et al. (2020). Occultation observations of Saturn’s rings with Cassini VIMS. Icarus, 344, 113356. https://doi.org/10.1016/j.icarus.2019.06.017

Ntinos, C. (2024). Development of super-resolution mapping methods for icy moons: Application on enceladus, in preparation for future observations by the ESA JUICE mission [Doctoral dissertation, Université Paris Cité]. https://theses.hal.science/tel-05044288v1

Penteado, P., & Trilling, D. (2015). Titanbrowse: Using a new paradigm for access to hyperspectral data. Second Planetary Data Workshop, 1846, 7044.

Peter, J. S., Nordheim, T. A., & Hand, K. P. (2024). Detection of HCN and diverse redox chemistry in the plume of Enceladus. Nature Astronomy, 8(2), 164–173. https://doi.org/10.1038/s41550-023-02160-0

Porco, C. C., Helfenstein, P., Thomas, P., Ingersoll, A., Wisdom, J., West, R., Neukum, G., Denk, T., Wagner, R., Roatsch, T., et al. (2006). Cassini observes the active south pole of Enceladus. Science, 311(5766), 1393–1401. https://www.science.org/doi/abs/10.1126/science.1123013

Postberg, F., Schmidt, J., Hillier, J., Kempf, S., & Srama, R. (2011). A salt-water reservoir as the source of a compositionally stratified plume on Enceladus. Nature, 474(7353), 620–622. https://doi.org/10.1038/nature10175

Robidel, R., Le Mouélic, S., Tobie, G., Massé, M., Seignovert, B., Sotin, C., & Rodriguez, S. (2020). Photometrically-corrected global infrared mosaics of Enceladus: New implications for its spectral diversity and geological activity. Icarus, 349, 113848. https://doi.org/10.1016/j.icarus.2020.113848

Scipioni, F., Schenk, P., Tosi, F., D’Aversa, E., Clark, R., Combe, J.-P., & Dalle Ore, C. (2017). Deciphering sub-micron ice particles on Enceladus surface. Icarus, 290, 183–200. https://doi.org/10.1016/j.icarus.2017.02.012

Seignovert, B., Le Mouélic1, S., Heslar, M., Cordier, D., Coutelier, M., & MacKenzie, S. (2023, January). PyVIMS (Version v1.0.4). Zenodo. https://doi.org/10.5281/zenodo.7535017

Shi, Z., Li, P., Jin, H., Tian, Y., Chen, Y., & Zhang, X. (2017). Improving super-resolution mapping by combining multiple realizations obtained using the indicator-geostatistics based method. Remote Sensing, 9(8), 773. https://doi.org/10.3390/rs9080773

Shkuratov, Y., Kaydash, V., Korokhin, V., Velikodsky, Y., Opanasenko, N., & Videen, G. (2011). Optical measurements of the Moon as a tool to study its surface. Planetary and Space Science, 59(13), 1326–1371. https://doi.org/10.1016/j.pss.2011.06.011

Spencer, J., Pearl, J., Segura, M., Flasar, F., Mamoutkine, A., Romani, P., Buratti, B., Hendrix, A., Spilker, L., & Lopes, R. (2006). Cassini encounters Enceladus: Background and the discovery of a south polar hot spot. Science, 311(5766), 1401–1405. https://www.science.org/doi/10.1126/science.1123013

Sun, M., & Chen, S. (2022). Deep learning-based super-resolution reconstruction and algorithm acceleration of Mars hyperspectral CRISM data. Remote Sensing, 14(13), 3062. https://doi.org/10.3390/rs14133062

Tacconi, C. S. A., Alessandra Buonanno, M. C., Olivier Grasset, A. H., Luciano Iess, E. K., Jérémy Leconte, J. L., Jesús Martín-Pintado, R. N., & Watson, D. (2021). Voyage 2050 final recommendations from the Voyage 2050 Senior Committee (Report). Senior Committee. https://www.cosmos.esa.int/web/voyage-2050

Taffin, C., Grasset, O., Le Menn, E., Bollengier, O., Giraud, M., & Le Mouélic, S. (2012). Temperature and grain size dependence of near-IR spectral signature of crystalline water ice: From lab experiments to Enceladus’ south pole. Planetary and Space Science, 61(1), 124–134. https://doi.org/10.1016/j.pss.2011.08.015

Tao, Y., Conway, S. J., Muller, J.-P., Putri, A. R., Thomas, N., & Cremonese, G. (2021). Single image super-resolution restoration of TGO CaSSIS colour images: Demonstration with perseverance rover landing site and Mars science targets. Remote sensing, 13(9), 1777. https://doi.org/10.3390/rs13091777

Tao, Y., & Muller, J.-P. (2016). A novel method for surface exploration: Super-resolution restoration of Mars repeat-pass orbital imagery. Planetary and Space Science, 121, 103–114. https://doi.org/10.1016/j.pss.2015.11.010

Tao, Y., Xiong, S., Muller, J.-P., Michael, G., Conway, S. J., Paar, G., Cremonese, G., & Thomas, N. (2022). Subpixel-scale topography retrieval of Mars using single-image DTM estimation and super-resolution restoration. Remote Sensing, 14(2), 257. https://doi.org/10.3390/rs14020257

Van Reeth, E., Tham, I. W., Tan, C. H., & Poh, C. L. (2012). Super-resolution in magnetic resonance imaging: A review. Concepts in Magnetic Resonance Part A, 40(6), 306–325. https://doi.org/10.1002/cmr.a.21249

Waite, J., Lewis, W., Magee, B., Lunine, J., McKinnon, W., Glein, C., Mousis, O., Young, D., Brockwell, T., Westlake, J., et al. (2009). Liquid water on Enceladus from observations of ammonia and 40Ar in the plume. Nature, 460(7254), 487–490. https://doi.org/10.1038/nature08153

Waite, J. H., Combi, M. R., Ip, W.-H., Cravens, T. E., McNutt Jr, R. L., Kasprzak, W., Yelle, R., Luhmann, J., Niemann, H., Gell, D., et al. (2006). Cassini Ion and Neutral Mass Spectrometer: Enceladus plume composition and structure. Science, 311(5766), 1419–1422. https://www.science.org/doi/10.1126/science.1121290

Waite, J., Glein, C. R., Perryman, R. S., Teolis, B. D., Magee, B. A., Miller, G., Grimes, J., Perry, M. E., Miller, K. E., Bouquet, A., et al. (2017). Cassini finds molecular hydrogen in the Enceladus plume: Evidence for hydrothermal processes. Science, 356(6334), 155–159. https://www.science.org/doi/10.1126/science.aai8703

Wang, C., Zhang, Y., Zhang, Y., Tian, R., & Ding, M. (2021). Mars image super-resolution based on generative adversarial network. IEEE Access, 9, 108889–108898. https://ieeexplore.ieee.org/document/9503382

Xie, H., West, R. A., Seignovert, B., Jewell, J., Kurth, W., & Averkamp, T. (2021). Compression algorithms for high-data-volume instruments on planetary missions: A case study for the Cassini mission. Journal of astronomical telescopes, instruments, and systems, 7(2), 028002–028002. https://doi.org/10.1117/1.JATIS.7.2.028002

Xu, X., Liu, J., Liu, D., Liu, B., & Shu, R. (2020). Photometric correction of Chang’E-1 Interference Imaging Spectrometer’s (IIM) limited observing geometries data with Hapke model. Remote Sensing, 12(22), 3676. https://doi.org/10.3390/rs12223676

Zhang, H., Yang, Z., Zhang, L., & Shen, H. (2014). Super-resolution reconstruction for multi-angle remote sensing images considering resolution differences. Remote Sensing, 6(1), 637–657. https://doi.org/10.3390/rs6010637

A Digital Shape Model (DSM)

During the incorporation of the DSM in the calculation of the navigation cubes for VIMS, we have investigated the difference associated with modeling the target body with a DSM instead of a triaxial ellipsoid. In Figure A1 we present the topography (difference between the radius of the DSM and the triaxial ellipsoid), while in Figures A2 and A3 we present the difference in the calculated incidence and emission angles respectively. The phase angle does not change significantly between the DSM and the ellipsoid as expected (only five pixels in the whole data set present an absolute difference greater than 0.001\(^\circ \)).

Figure A1: Topography of the DSM model of Enceladus.

Figure A2: Difference in incidence angle between the DSM and ellipsoid shape models in Enceladus global color map.

Figure A3: Difference in emission angle between the DSM and ellipsoid shape models in Enceladus global color map.

B VIMS data wavelength shift and saturation

As introduced in Section 3.1, during the course of the Cassini mission it has been discovered that the VIMS instrument suffers from a wavelength shift, the origin of which has still not been entirely understood. The VIMS wavelength shift is fully documented in Clark et al. (2018). The original VIMS calibration pipeline used a wavelength grid and the Full Width at Half Maximum (FWHM) as measured during its thermal vacuum testing before launch (Brown et al., 2004). This means that the VIMS response function in previous Radiometric Calibrations (notably RC17), that was measured preflight on the ground, does not take into account the gradual shift of wavelengths, resulting in artifacts in radiometrically calibrated spectra that increase with time (Clark et al., 2018). In contrast, the last iteration of the Radiometric Calibration (RC19) takes into account the observed wavelength shift of VIMS, making the RC19 a time-dependent radiometric calibration. Each calibrated cube contains in its header a list (annotated with the keyword “Center”) that contains the new shifted wavelengths associated with each spectral band saved in the cube. These wavelengths deviate more and more as the mission progresses with the first band drifting from a central wavelength of 0.88461 \(\mu \)m, at the beginning of the mission, to 0.89461 \(\mu \)m, at the end of the mission, and the last band from 5.12382 \(\mu \)m to 5.13382 \(\mu \)m, hence more than half of the typical spectral band FWHM in the VIMS infrared channel. In the context of implementing our SR approach, which requires the combination of the spectral information of overlapping pixels acquired at different times, we need to make sure that we are not mixing different wavelengths contained in the same band as the result of the shift. The solution is for each spectrum, namely for each pixel used in our maps, to calculate a linear interpolation and resample it on a standard, reference wavelength grid defined at the beginning of the mission (Clark et al., 2018), as already implemented on VIMS cubes of Titan by Le Mouélic et al. (2019a). Since the wavelength associated with the first band shifts to longer wavelengths as the mission progresses and in order to interpolate the spectra from the whole mission over the same wavelength range, we remove the first standard wavelength and interpolate the spectra over a range of 0.900753 \(\mu \)m to 5.122500 \(\mu \)m with 255 bands. During the wavelength shift correction, from a total of 201 044 pixels we identify and remove 8709 pixels that have recorded no spectrum (all 256 infrared bands missing from the calibrated cube). From the remaining 192,335 pixels, 133 977 pixels have bands missing, 59\(\%\) of which are only missing up to 10 bands. In order to maximize the number of pixels, we decided to include them in the data set by maintaining the parts of the original spectrum that are missing, ensuring that we do not introduce artifacts through the interpolation of the missing bands. One of the advantages of the SR pipeline we have developed is that it can consolidate spectra with internally missing values by averaging at each wavelength only the available information. This in turn guarantees maximum utilization of the available data.

Another issue that was discovered upon careful examination of the VIMS Cassini data set is that of pixel saturation. The VIMS data set of Enceladus, similarly to the other icy moons observed by Cassini VIMS (Dalle Ore et al., 2021), contains RC19 calibrated cubes (.cub files), which can exhibit saturation in the shorter wavelengths of the infrared channel. This is particularly relevant for a highly reflective surface, as in the case of Enceladus. Contrary to previous calibration pipelines, where the saturated part of the spectrum has been removed (RC17), the current available data set exhibits a reconstruction of the spectrum which introduces spectral artifacts, which left untreated will inexorably lead to the introduction of artifacts into the SR cartographic products and any further spectral analysis we will implement using them. Hence it is of outmost importance, before moving to the next steps of our methodology, to guarantee the exclusion from the data set of the saturated parts of the spectrum. The initial strategy tested relied on the horizontal exclusion of a specific range of the spectrum at shorter wavelengths, where saturation effects seemed to be most prominent, for all pixels included in the data set. However, this approach greatly increases the potential useful parts of the spectra that are being discarded, since it does not consider the dependency of the saturation on the exposure time. Moreover, a relatively high exposure time can also not be used as an indicator of saturation for a given pixel, because in the same cube, where all pixels have the same exposure time, there can be found a variety of illumination and observation conditions, resulting in both saturated and non-saturated pixels. We should also consider that shorter wavelengths can be critical in implementing radiative transfer models to retrieve surface composition and physical properties (Mermy et al., 2023) and therefore removing completely these parts of the spectrum from our cartographic products can greatly limit their applicability for spectral analysis. As described by Mccord et al. (2004) VIMS uses 12-bit digitization, and therefore can record up to \(2^{12}-1=4095\) digital numbers (DN). However, taking into account the background noise (bns) that needs to be subtracted from the total capacity of the detector, saturation occurs at \(\approx \)3550 DN for the VIMS-IR channel. This approach properly quantifies the saturation limit instead of relying on a qualitative evaluation of each spectrum separately, but also fails to take into account the dependency of the saturation limit on wavelength and exposure time. Instead of using a mean background noise and hence the same saturation limit across different pixels and wavelengths, we propose to subtract from the maximum capacity of the detector the background noise values that can be found in the backplanes of each VIMS cube (limit=4095-bns). In this way we ensure that the implemented limit is specifically tailored to each pixel as a function of wavelength and therefore ensures maximum identification/utilization of the non-saturated part of each spectrum. In Figure B1 we present the raw (bottom left) and corresponding calibrated data (top left) for an example of a saturated spectrum from cube=1487294146_1, sample=45, line=15. The top left figure presents the results of removing the saturated part of the spectrum, with the continuous line representing the remaining spectrum after removal of the saturated part and the dotted line the original calibrated spectrum, for three different saturation limits: constant 3550 DN, limit=4095-bns and 0.97\(\cdot \)(4095-bns). The spectra are shifted vertically for clarity. For reference, we also present on the right panels of the same figure the raw and corresponding calibrated data for a non-saturated spectrum from cube=1489049741_1, sample=1, line=24. Taking the more conservative approach (considering the artifacts in the spectrum presenting as absorption bands can lead to erroneous detection of species) we add a 3% margin on the saturation limit as defined above to ensure the exclusion of the saturated part of the spectrum.

Figure B1: Examples of saturated (left) and non-saturated pixels (right) in the data set: cube=1487294146_1, sample=26, line=15 and cube=1489049741_1, sample=1, line=24 respectively. The top row presents the calibrated spectra and the bottom row their corresponding raw data. The top left figure presents the results of removing the saturated part of the spectrum (continuous line represent the remaining spectrum and the dotted line the original calibrated spectrum) using three different saturation limits (constant 3550 DN, 4095-bns and 0.97\(\cdot \)(4095-bns)) shifted vertically for clarity. The bottom left figure presents the corresponding raw data of the saturated spectrum with the red dotted line marking the constant saturation limit of 3550 DN, the dotted orange line the 4095-bns and the orange envelope the 3% margin. The embedded panel in the bottom left figure presents a zoomed-in view of the saturated raw signal and its relationship with the different saturation limits examined.

By implementing the above mentioned strategy to the entirety of the Enceladus data set, we can produce global cartographic products free of saturation in every wavelength, while maximizing the number of pixels and wavelength range preserved. In Figure B2 we demonstrate the global statistics, hence the significant influence, of the saturation issue by presenting the percentage (left vertical axis) and absolute number (right vertical axis) of the saturated pixels for different exposure groups (from top to bottom: all pixels, exposure > 0.6 sec, 0.4 sec < exposure \(\leq \) 0.6 sec, 0.2 sec < exposure \(\leq \) 0.4 sec and exposure \(\leq \) 0.2 sec). The identification of the saturated part of the spectrum is presented for two different saturation limits: for the constant across wavelength limit = 3550 DN (continuous blue line) and for the proposed dynamic limit = 0.97\(\cdot \)(4095-bns). From Figure B2 we can deduce that for short wavelengths up to more than 50% of the pixels in the data set can exhibit saturation, while the effect of saturation can extend up to longer wavelengths up to 2.5 \(\mu \)m. At the same time, the binning of the pixels for different exposure time ranges reveals the relationship between saturation and exposure time, with both the main lobe at short wavelengths and secondary lobes (centered at 1.8 \(\mu \)m and 2.25 \(\mu \)m) gradually diminishing the lower the exposure time. Finally, we note that implementing a constant DN limit fails to capture the saturation occurring around 1, 1.25 and 1.34 \(\mu \)m, which can only be addressed by implementing the dynamic filter we are proposing.

Figure B2: Percentage (left vertical axis) and absolute number (right vertical axis) of the saturated pixels for different exposure groups (from top to bottom: all pixels, exposure > 0.6 sec, 0.4 sec < exposure \(\leq \) 0.6 sec, 0.2 sec < exposure \(\leq \) 0.4 sec and exposure \(\leq \) 0.2 sec) and different saturation limits (for the continuous blue line limit = 3550 DN and for the dashed orange line limit = 0.97\(\cdot \)(4095-bns)).

C Photometric correction coefficients

As described in Section 3.2, we perform the photometric correction for each wavelength separately and therefor calculate the coefficient of the linear fit (\(a_0, a_1\)) as a function of wavelength (\(\lambda \)). In Figures C1 and C2 we present these coefficients as a function of wavelength. From Figure C1 the physical meaning of the constant in the linear fit becomes evident as the consolidated equigonal albedo at zero phase angle for the whole surface of Enceladus (\(a_0\)), exhibit the characteristic absorption bands of a surface dominated by water-ice. In Figure C3 we present the correlation of the coefficients of the linear fit (\(a_0, a_1\)) with the color map representing the corresponding wavelength.

Figure C1: Photometrical correction model intercept \(a_0\) of the best linear fit as a function od standard wavelength \(\lambda \) (also identified as the equigonal albedo at zero phase angle).

Figure C2: Photometrical correction model slope \(a_1\) of the best linear fit as a function of standard wavelength \(\lambda \).

Figure C3: Photometrical correction model slope \(a_1\) with respect to the photometrical correction model intercept \(a_0\) of the best linear fit. The color map represents the corresponding standard wavelength \(\lambda \)

D Effective sampling

Quantifying the spatial sampling for non-rectangular pixels is not a straight-forward process and therefore we have defined a relevant metric for the effective sampling of the SR maps. The effective sampling can be thought of as the smallest element in the SR map, that introduces a difference in the brightness value. In the theoretical case, this corresponds to a SR map made up of polygons, each one of which corresponds to a common overlapping area of pixels from the original data set. The super-cells consist of an interpolation grid and cannot be used to define the spatial sampling of the resulting SR calculation. Hence we need to calculate the effective sampling of each super-cell in our grid. In the general case, every super-cell is associated with two types of overlapping pixels: 1) fully overlapping pixels, that completely contain within their borders the super-cell and 2) partially overlapping pixels, that have a common overlapping area with the super-cell without fully containing it. In Figure D1 we present an example of the fully (top panel in red) and partially (bottom panel in orange) overlapping pixels for super-cell with index 691,639. In green we note the super-cell while the shaded blue area corresponds to the common overlapping area of the fully overlapping pixels.

Figure D1: Example of fully (in red) and partially (in orange) overlapping pixels for super-cell with index 691,639 (top and bottom panel respectively). In green we note the super-cell while the shaded blue area corresponds to the common overlapping area of the fully overlapping pixels.

Opting for the most conservative approach, the calculation of the effective sampling for each super-cell of the SR map, is carried out by evaluating, from the list of overlapping pixels related to each super-cell, the pixels that fully overlap with it. By definition, the pixel-polygons that fully contain the super-cell will always define a common overlapping polygon (shaded blue in the example of Figure D1). We define this common overlapping polygon as the “effective” polygon and the maximum distance between any two of each of its edges the effective sampling (as calculated for the VIMS pixels in Section 3.1) associated with the specific super cell. From the 1 036 800 super-cells in the 1440\(\times \)720 grid 37 have only partially overlapping pixels and therefore do not define a common overlapping polygon. In these cases we set \(gain=1\), meaning that SR preserves the original VIMS sampling of the specific area. Moreover, due to the conservative approach selected in calculating the effective sampling, 20 423 super-cells (less than 2% of the total super-cells) result in \(gain<1\), due to the exclusion of the partially overlapping pixels from our calculation, which would further increase the gain from the SR method. In this case we also define \(gain=1\). This does not mean that these super-cells are less resolved than the most resolved VIMS observation over that same area, rather than the calculation of their effective sampling is limited by the definition of the metric used.

E VIMS radiometric error

Throughout the implementation of our pipeline, we are interested not only in the spectral values themselves, but also in their accompanying error. Various methodologies have been employed in the literature to quantify the radiometric error of the VIMS data, such as the adoption of the error as a fixed percentage of the signal (Mermy et al., 2023) or statistical evaluation of the error as the standard deviation of the ensemble of spectra contributing to specific bins (Coutelier et al., 2021; Filacchione et al., 2010). In this work we have followed the approach by Es-Sayeh et al. (2023) to calculate the radiometric noise (\(\sigma _{r_f}\)) associated with the spectral values of each pixel in the Enceladus VIMS data set by adding the various sources of noise. The radiometric error of the radiance factor of the calibrated cubes can then be propagated to the photometrically corrected and SR radiance factor values for each pixel. From the theory of statistics we know that the error of a function \(f\)(\(\sigma _f\)) with \(n\) variables \(x_i\), is a function of the error of each parameter \(\sigma _{x_i}\) and the partial derivatives \(\frac {\partial f}{\partial x_i}\) as follows (Ku, 1969):

\begin{equation} \label {eq:propagation} \sigma _{f}=\left [ \sum _{i=1}^{n} \left ( \frac {\partial f}{\partial x_i} \right )^2 \sigma _{x_i}^2 \right ]^{1/2}\tag {E1} \end{equation}

By implementing Equation E1 to Equation 11 and Equation 12, we get respectively:

\begin{equation} \sigma _{r_{f_{cor}}}= r_{f_{cor}}\left [ \left ( \frac {\sigma _{r_f}}{r_f} \right )^2 + \left ( \frac {\sigma _{a_0}}{a_0} \right )^2 + \frac {\sigma ^2_{a_0}+\sigma ^2_{a_1}g^2}{\left (a_0+a_1g\right )^2} \right ]^{1/2}\tag {E2} \end{equation}

and

\begin{equation} \sigma _{SR}=\left [ \frac {\sum _{i}^{} w_i^2 {\sigma ^2_{r_{f_{cor_{i}}}}}}{\left (\sum _{i}^{} w_i \right )^2} \right ]^{1/2}\tag {E3} \end{equation}

where \(\sigma _{r_{f_{cor}}}\) the one sigma error of the photometrically corrected radiance factor (\(r_{f_{cor}}\)), \(\sigma _{a_0}\) and \(\sigma _{a_1}\) the one sigma errors of the \(a_0\) and \(a_1\) coefficient of the linear fit respectively as performed in Section 3.2 and \(\sigma _{SR}\) the one sigma error of the radiance factor of the SR maps. When pixels have missing values in the calculated one sigma error, they are set to zero in the one sigma error of the super-cells overlapping with them. In Figure E1 we present a global color map of the one sigma error of the SR radiance factor (\(\sigma _{SR}\)) at (1.80401 \(\mu \)m) associated with each super-cell in our final SR cartographic product. We observe that SR operation does not significantly increase the one sigma error, which remains at the same level of the original radiometric noise (\(\sigma _{r_{f}}\)), namely less than 1% of the signal. We recognize that this approach can potentially underestimate the radiometric noise and we are in the process of updating our method.

Figure E1: Global SR map of Enceladus at grid size 1440\(\times \)720 cells with the color scale representing the one sigma error of the SR radiance factor (\(\sigma _{SR}\)) at 1.80401 \(\mu \)m.