A GLOBAL SHAPE MODEL FOR SATURN’S MOON ENCELADUS FROM A DENSE PHOTOGRAMMETRIC CONTROL NETWORK

A planetary body’s global shape provides both insight into its geologic evolution, and a key element of any Planetary Spatial Data Infrastructure (PSDI). NASA’s Cassini mission to Saturn acquired more than 600 moderateto high-resolution images (<500 m/pixel) of the small, geologically active moon Enceladus. The moon’s internal global ocean and intriguing geology mark it as a candidate for future exploration and motivates the development of a PSDI. Recently, two PSDI foundational data sets were created: geodetic control and orthoimages. To provide the third foundational data set, we generate a new shape model for Enceladus from Cassini images and a dense photogrammetric control network (nearly 1 million tie points) using the U.S. Geological Survey’s Integrated Software for Imagers and Spectrometers (ISIS) and the Ames Stereo Pipeline (ASP). The new shape model is near-global in extent and gridded to 2.2 km/pixel, ~50 times better resolution than previous global models. Our calculated triaxial shape, rotation rate, and pole orientation for Enceladus is consistent with current International Astronomical Union (IAU) values to within the error; however, we determined a new prime meridian offset (Wo) of 7.063. We calculate Enceladus’ long-wavelength topography by subtracting the best-fit triaxial ellipsoid from our shape model. The result is comparable to previous global models but can resolve topographic features as small as 5-7 km across in certain areas. To evaluate the spatially varying quality of the model, we calculate the point density (variable from 5 to more than 50 points per pixel), normalized median absolute deviation of the points within each pixel (typically less than 100 m), and the minimum expected vertical precision of each point (ranging from 29 m to 2 km).


The importance of planetary shape and topography
The global shape and topography of a planetary body provides insight into its interior structure, orbital dynamics, and thermal history. Characterizing a body's shape is therefore necessary for understanding its evolution over time. On a more practical level, topography is also critical to the development of numerous higher order data products (e.g., orthoimages), and along with geodetic control and orthoimages, constitutes one of the three foundational data sets necessary for a Planetary Spatial Data Infrastructure (PSDI) (Laura et al. 2017). PSDI is the collection of spatial data, access mechanisms, standards, policies, and data users that enables effective application of spatial data to scientific problems (Laura et al. 2018). The three foundational products can be used in combination with other data sets to generate numerous higher order spatial framework data. An effective PSDI must also include a rigorous characterization of product quality and uncertainty for all data products.
Here we describe a new global shape and topography (i.e., the deviation in shape from the best-fit ellipsoid) model for Saturn's moon Enceladus, which provides the final element of a PSDI for this moon (see section 1.2). In section 2, we describe our photogrammetric approach to shape model generation, which conceptually follows that of Archinal et al. (2005) and Becker et al. (2016), who applied similar techniques to the Moon and Mercury, respectively. Our results are described in section 3, and a characterization of the spatially variable quality of the shape model follows in section 4. The shape model is publicly available as both a point cloud and gridded (interpolated) products at the Annex of NASA's Planetary Data System (PDS) Cartography and Imaging Science Node, which is supported by the U.S. Geological Survey (USGS), where we also provide updated pointing kernels (ck), and additional metadata (measures of quality).

Why Enceladus?
Despite its small size (251 km in radius), Saturn's moon Enceladus is geologically active, with plumes of ice and gas emanating from four warm, parallel fractures at its south pole (Hansen et al. 2006, Porco et al. 2006, Spencer et al. 2006 (Fig.  1). Salts detected in the plume by NASA's Cassini spacecraft suggest the fractures are conduits to a liquid reservoir (Postberg et al. 2009) at a depth of just 4-14 km (Hemingway and Mittal 2019). Detection of a large physical libration strongly suggests that the liquid reservoir is global in extent (Thomas et al. 2016) with an average ice shell thickness of 20-30 km (McKinnon 2015;Hemingway and Mittal 2019). Although the south pole is the locus of current activity, much of Enceladus' surface is geologically young, as indicated by extensive tectonic deformation and low impact crater density (Crow-Willard and Pappalardo 2015). Many of its craters have been extensively modified, indicating unusually high heat flow in Enceladus' past (Bland et al. 2012). The current elevated heat flow at the south pole, global ocean, and past periods of high heat flow must be maintained by tidal heating, but the exact mechanism by which sufficient heat is generated within Enceladus remains unclear (e.g., Meyer and Wisdom 2007, Roberts and Nimmo 2008, Roberts 2015, Souček et al. 2019. Given Enceladus' remarkable geology, its status as a confirmed ocean world, and its tendency to spew its interior into space (where it can "easily" be sampled by spacecraft), the moon is a focus of scientific interest and a target for future exploration. NASA's Cassini mission returned a wealth of image data of Enceladus (section 2.2). Unfortunately, inaccurate knowledge of the spacecraft position and pointing when the images were acquired make the data challenging to use: image locations on the surface can be inaccurate by several to tens of kilometers. This issue can be addressed by the development of a Planetary Spatial Data Infrastructure (PSDI) to enable the scientific community to effectively utilize geospatially accurate data returned by Cassini.
Previously, Bland et al. (2018) described a photogrammetric control network for Enceladus that included a total of 586 Cassini Imaging Science Subsystem (ISS) (Porco et al. 2004) images of Enceladus (additional images were added later to bring the total number of images to 621). That work provided two of the three foundational PSDI products: photogrammetric control and orthorectified images (section 1.1). We distinguish between photogrammetric control, in which there is no "ground," and geodetic control, which uses laser altimetry data or some other means to determine absolute control. To date, no such data have been acquired for Enceladus, so photogrammetric control is the highest level of control achievable. We note that the control network of Bland et al. (2018), following IAU recommendations (Archinal et al. 2018, Table 2), did fix the location of the crater Salih at 5 o W, which defines Enceladus' longitude system and provides one geodetic reference point. We also note that the orthoimages produced by Bland et al. (2018) are projected to Enceladus' triaxial shape rather than local topography. Figure 1: The surface of Saturn's small moon Enceladus (251 km radius) is highly deformed, indicating a complex geologic history. Plumes of dust (mostly ice) and gas (mostly water vapor) emanate from its south pole (inset). Images were acquired by the ISS on NASA's Cassini spacecraft and are publicly available through the Planetary Data System (PDS).

Previous shape models
Several global, or semi-global topographic data sets have previously been created for Enceladus. These data sets were derived from two different techniques: limb fitting, and stereo imaging. Porco et al. (2006) provided the first measurements of Enceladus triaxial shape from Cassini images of Enceladus' limb (the approach was described more fully by Thomas et al. (2007)), and these measurements were subsequently updated by Thomas (2010) at the end of Cassini's nominal mission, with the only change being a decrease in the uncertainty in the tidal (a) semi-axis. The Thomas (2010) analysis derived the triaxial shape of Enceladus (a, b, c semi-axes) to within 200-300 m, but did not report on Enceladus' regional or local topography. Schenk and McKinnon (2009) provided the first semi-global topography model for Enceladus by using stereo imaging to derive topography for ~50% of Enceladus' surface. The model has a relatively small horizontal pixel scale of 200-950 m/pixel and a stated vertical precision of 50 to 140 m. Significantly, their analysis identified six large-scale depressions, 800-1500 m deep, that are not correlated with obvious surface geology, and confirmed the presence of a previously identified south polar depression (cf. Porco et al. 2006;Collins and Goodman 2007). The depressions might correspond to thin regions of the ice shell (Schenk and McKinnon 2009) or local regions of convectiondriven compaction (Besserer et al. 2013).
The first truly global shape model of Enceladus was provided by Nimmo et al. (2011), who performed a spherical harmonic expansion of the Thomas (2010) limb data up to degree 8. Although the degree-8 expansion corresponds to a spatial scale of just 45 o or ~200 km at the equator (λ~2πR/L, where R is Enceladus' mean radius of 251 km, and L is the harmonic degree), the model enabled a more detailed evaluation of Enceladus' long-wavelength global topography. Where they overlap, the shape derived by Nimmo et al. (2011) is in general agreement with the stereo model of Schenk and McKinnon (2009) and generally confirms, with some variation, the existence of the large-scale topographic depressions.
Using a similar approach, but with the benefit of a more extensive data set, and the additional use of control points, Tajeddine et al. (2017) performed a spherical harmonic expansion of Enceladus limb data to degree 16, improving the spatial "resolution" by a factor of 2 relative to the Nimmo et al. (2011) model (to 22.5 o or ~100 km at the equator). The results are largely consistent with the lower degree shape model of Nimmo et al. (2011) and again reveal a chain of topographic basins that Tajeddine et al. (2017) argue are evidence of a relic equator, and thus true polar wander. However, the two models are different enough in their details that analysis of their implications for Enceladus' interior structure yields modestly different results (Hemingway and Mittal 2019).
The existing shape models for Enceladus thus have their strengths and weaknesses. The spherical harmonic expansions of Nimmo et al. (2011) and Tajeddine et al. (2017) are global in extent, but relatively low resolution. In contrast, the stereo model of Schenk and McKinnon (2009) is 2 orders of magnitude higher resolution (at least locally) but covers just a portion of Enceladus. Here we describe a new shape model for Enceladus created with a photogrammetric approach that balances resolution and coverage, providing a shape model that is nearly global (>92% coverage), relatively high resolution (2 pixel/degree or 2.2 km/pixel at the equator), and vertically precise (root mean square radius uncertainty of 57 m).

Overview
We have created a global shape model of Enceladus using 625 Cassini ISS images of Enceladus following the conceptual approach applied to Mercury MESSENGER data by Becker et al. (2016). We first established a dense network of image tie points and photogrammetrically solved for point latitude, longitude, and radius using a least squares bundle adjustment. The resulting point cloud was then interpolated to a regularly gridded 2.5D global shape model. We subtracted a best-fit triaxial ellipsoid from the model to reveal the moon's long-wavelength topography (local deviation from the ellipsoid). The data set and each of these steps are described in more detail below. We describe the shape model quality in section 4.

The Cassini image data set
NASA's Cassini mission spent 13 years (2004 to 2017) in the Saturn system. The long mission duration and high-performing spacecraft enabled 23 targeted flybys of Enceladus. These flybys, and an additional ~30 more-distant flybys, returned more than 22,000 images, including more than 600 images with a pixel scale better than 500 m/pixel, from Cassini's ISS cameras. The ISS included a 2-m focal length, 0.35 o field of view (FOV) narrow angle camera and a 0.2-m focal length 3.5 o FOV wide angle camera that shared a 1024 by 1024 pixel charged couple device (CCD) detector (Porco et al. 2004).
Although Cassini achieved global image coverage, the orbital dynamics of the mission, the scientific focus on Enceladus' south pole, and the seasonal illumination constraints resulted in highly variable image quality across the satellite. For example, only a few images of the northern polar region were acquired due to it being dark during northern winter until late in the extended mission. Image coverage is also sparse in the trailing hemisphere due to the orbital dynamics of the mission: Cassini only flew over these areas early in the mission. Even where multiple flybys returned repeat coverage of the same areas, differences in flyby altitude and the location of the closest approach of each flyby resulted in large differences in image pixel scale at the same location. Illumination and viewing geometry also varied substantially between flybys. Many images include the moon's limb, terminator, or both (incidence and emission angles up to 90 o ) and phase angles reach a maximum of 178 o .
The multiple-flyby nature of the data set poses a substantial challenge for the image-to-image matching that is required to create a photogrammetric control network for the purposes of improving image locations and/or deriving shape information. High incidence angles create problematic shadows and high emission angles distort features and compromise measurement accuracy. Changes in illumination also substantially affect the appearance, and even visibility, of surface features (see  for examples). Matching images with substantially different pixel scales is especially challenging.
Despite these challenges, Bland et al. (2018) created a relatively sparse (compared to the present work) tie-point network and photogrammetrically updated the location of 586 images with pixel scale between 50 and 500 m/pixel, phase angle less than 120 o , and filter settings of CLR, GRN, IR3, and UV3 (Porco et al. 2004). Least squares bundle adjustment resulted in root-meansquare (RMS) residuals of 0.45 pixels, corresponding to ground point uncertainties of 66, 51, and 46 m in latitude, longitude, and radius, respectively. Subsequently, 35 images of the north pole acquired late in the Cassini mission were tied to the network, creating a global set of 621 images. Additional high-resolution images were then tied to the global basemap on an image-byimage basis. The resulting mosaics are available through the USGS Astropedia data portal at https://astrogeology.usgs.gov/search/map/Enceladus/Cassini/En celadus_Cassini_mosaic_global_110m. These updated images form the basis for the work described here. Although points from the sparse network itself were not utilized, the updated image locations facilitated more efficient matching during construction of our new dense tie-point network.

The dense photogrammetric control network
To create a dense, global network of tie points we used an imageby-image approach such that every image in the data set has its own individual tie-point network. These 621 individual networks were then merged to create the global network (Fig. 2). Four additional images were subsequently added to provide additional coverage in the trailing hemisphere, bringing the total number of images used to 625.
For every image, we first established a regularly spaced grid of tie points with an initial density of one point every 20 line/sample (Fig. 2b,c). Increasing the density lead to computational difficulty due to the network size. Automated tie-point matching was performed between the reference image in question and every other overlapping image. To do so, we used an area-based approach with a maximum correlation algorithm via ISIS's pointreg application (Garcia et al. 2015). A weighted centroiding approach was used to improve the match to subpixel accuracy. Matches below a Goodness-of-Fit threshold were ignored, and false matches were removed based on bundle adjustment residuals (section 2.4). Every successfully matched image provided an image "measure" for that tie point. Most tie points are tied to multiple other images, providing numerous measures and substantial network "depth," which yields a photogrammetrically robust solution ) and increases stereo strength. Because every image is associated with its own small network, matching between each overlapping image pair is attempted twice. That is, image-B is included in image-A's network and image-A is include in image-B's network. Tie points are identified in both directions but because images often only partially overlap and are of different resolution, points from each network occur at different line/sample locations and are therefore generally unique. The 625 individual image networks were combined into a single global network using ISIS's cnetmerge application. This approach has the benefit of correlating tie-point density with depth of image coverage at a location. That is, where many images overlap, point density is high, and where image coverage is limited, point density is lower. Spatial variation in tie-point density is therefore naturally indicative of local data density (although not necessarily quality, see section 4). The resulting tie-point network contains 892,457 points and more than 30 million image measures. The network is nearly ~90x more spatially dense than the network used by Bland et al. (2018).
Our approach to control network generation differs from that of Becker et al. (2016) who used the OpenCV library (opencv.org) of feature-based matching algorithms (as implemented in ISIS's findfeatures application) to create a dense network for Mercury. We found that our initial attempts at feature-based matching did not provide adequate spatial coverage (i.e., some areas had insufficient tie-point density), and that matching images with highly variable illumination conditions were more challenging then when area-based matching was used. However, given the large number of detector and extractor algorithms available, utilizing feature-based matching on challenging flyby data sets warrants additional investigation.

Bundle adjustment
Once the global network of tie points was established and point matching was complete, we performed a least squares bundle adjustment (Brown 1958) using the ISIS application jigsaw (Edmundson et al. 2012) to update camera pointing and image locations. All observations were weighted equally. In the adjustment, we solved for the camera pointing (three rotation angles) and the 3D coordinates (latitude, longitude, and radius) of each tie point in the network. We did not solve for spacecraft position, as doing so added additional parameters without improving the solution. A priori camera pointing (angles) was constrained to ± 0.5 o , and a priori point radii were constrained to 1000 m. In our final solution, we also constrained point latitude and longitude to 1000 m. This additional constraint improved convergence of our bundle solution and decreased point uncertainties. We used jigsaw iteratively to identify and remove false matches from the network. Measures with very high residuals after bundle adjustment often indicate a false match. These measures are removed from the network, and a new bundle adjustment is performed. This process is repeated until either all point residuals are within an acceptable tolerance or visual inspection indicates quality matches despite some high individual residuals. Our final bundle adjustment resulted in RMS ground point uncertainties of 37 m, 36 m, and 57 m, in latitude, longitude, and radius, respectively. The total RMS uncertainty in image location was 0.3 pixels. A more detailed description of the shape model quality is provided in section 4. Figure 2: Illustration of our approach to control network generation for Enceladus. a) The blue rectangle is the footprint of a single ISS image (our example image is 167 km across). Also shown are the outlines of every image footprint that overlaps our example image (equirectangular projection, west longitude). b) Each image has its own line/sample tie-point network (green points on image footprint, which is the highlighted footprint in 'a'). c) The same tie points as seen on the actual unprojected image. d) Merging two networks locally increases network density. e) The final distribution of tie points on our example image from all the individual merged networks (see panel a).
Once a final bundle adjustment has been performed, the collection of 3D coordinates of each tie point (latitude, longitude, radius) constitute a point cloud that defines Enceladus' shape. As a final step, we calculate the minimum expected vertical precision (EP) of each point (see section 4) and filter out any points with EP greater than ~2 km, which is the maximum topographic relief on Enceladus. Snapshots of the point cloud are shown in Fig. 3.

The triaxial shape solution
The ISIS jigsaw application can also be used to calculate the bestfit triaxial shape and orientation of Enceladus. Doing so requires specifying an initial "guess" and uncertainty values for each value. In total, we solved for the triaxial shape (a, b, c semi-axes), average radius, spin rate, pole orientation (right ascension and declination), and prime meridian offset (Wo). Because images were acquired over a range of true anomaly, we assumed the effect of libration on our solution is small. The results are described in section 3.

Interpolation to a 2.5D global shape model
The point cloud resulting from our photogrammetric solution is not uniform in spatial density. We therefore interpolated to a uniformly gridded 2.5D shape model using the Ames Stereo Pipeline (ASP) point2dem application (Beyer et al. 2018). We chose a grid spacing of 2 pixel/degree (ppd) (2.2 km/pixel at the equator) based on an analysis of areal coverage and point density (points per pixel) (Figure 4). At 2 ppd we achieve 92.5% coverage of Enceladus. Decreasing the resolution to 1 ppd reduces the resolution by a factor of 2 (4.4 km/pixel) but provides only a small increase in coverage (from 92.5% coverage to 96.3%). Alternatively, increasing the resolution to 3 ppd decreases the coverage by 7 percentage points (from 92.5% coverage to 85.6%), which we believe is unacceptably low for a global product. We also find that a grid spacing of 2 ppd yields an average point density of 10.6 points per pixel (Fig. 4b), which enables more robust evaluation of point statistics in the gridded product (section 4). Using a resolution of 3 ppd decreases the number of points per pixel by a factor of 2. A mean density of 1 point per pixel is achieved at 6 ppd (733 m/pixel) but such a high resolution only provides coverage over 57% of Enceladus. Thus, a grid spacing of 2 ppd provides a good balance of resolution, coverage, and point density.
Radius values at each pixel were determined by a Gaussian weighted average of the points within each pixel. We used a 1pixel search window in point2dem, which eliminates any smoothing that results from including points outside the pixel in the average. No additional gap-filling or smoothing was applied. The gridded shape model is shown in Fig. 5.

Calculating Enceladus' long-wavelength topography
In addition to Enceladus' global shape, which is strongly triaxial, we also calculated the long-wavelength topography, which we define as local deviation in radius from the best-fit triaxial shape. To do this, we calculate the expected radius at each point (latitude/longitude) in our point cloud assuming the smooth triaxial shape derived in section 2.5 and then subtract it from the local ("true") radius value for that point derived from our bundle solution. The result is a set of latitude/longitude/elevation points, where the elevation is relative to the best-fit ellipsoid. As with the shape, we interpolated the topography data to a fixed grid using point2dem with a 2 ppd grid size and Gaussian weighting. The gridded topography is shown in Fig. 6.

Enceladus' triaxial shape
Direct interpolation of the point cloud to a gridded 2.5D product yields a nearly continuous shape model for Enceladus (Fig. 5). Data gaps (mostly in the trailing hemisphere) appear as black spots or regions where the underlying grayscale basemap shows through. As expected, the model is dominated by Enceladus' triaxial shape which stems from the tidal distortion of the moon by Saturn and results in a long tidal axis (a) that is oriented toward/way from Saturn, a short polar axis (c), and an intermediate axis (b). The jigsaw solution (section 2.7) yields a, b, and c semi-axes of 256.3 km, 251.2 km, and 248.3 km, respectively. The values are consistent with the previous determination by Bland et al. (2018), who used a similar method but much sparser network, and are within the uncertainty of the current IAU values , their table 5).
Figure 5: A) The gridded 2.5D global shape model of Enceladus. The model is dominated by the moon's triaxial shape. Values are relative to a mean sphere with radius 251.5 km (total range of +5300 to -3900 m): yellows are high, and blues are low. Equirectangular projection in west longitude. B) As in A but overlain on an image mosaic of Enceladus.
Again using jigsaw, we calculate a spin rate that is effectively identical to the IAU value ( Bland et al. (2018). Jigsaw also yields Wo of 6.088 o , which is somewhat smaller than the IAU value of 6.32 o   Table 2). Following Bland et al. (2018), we also calculate Wo using a more physical approach, in which we set a tie point within the center of the crater Salih (which is located at 5 o W longitude and defines the Enceladus coordinate system) and determine the longitude adjustment necessary to maintain Salih's correct position. This approach yielded a Wo of 7.063 o , similar to that of the sparse network (7.089 o ) reported in Bland et al. (2018), but significantly larger than the IAU value. The reason for the inconsistency between the two approaches remains unclear and will require further validation of the jigsaw body solver. Figure 6 shows our gridded 2.5D topographic model. Here the topography is relative to our best-fit ellipsoid. Several features are notable in the topography. The expected asymmetry (based on the earlier shape determination, e.g., Porco et al. (2006)) between the northern and southern hemispheres is apparent. The south pole is 200-400 m below the triaxial shape and ~500-1000 m lower than the north pole. The model also clearly shows a chain of depressions stretching across the anti-Saturnian hemisphere from southeast to northwest. The deepest depression (near 150 o W, -12 o S) sits ~800 m beneath the reference ellipsoid, and more than 1 km below the topography directly to the west. The depressions become shallower moving to the northwest, with the northern-most depression (near 210 o W, 45 o N) having a maximum depth of ~500 m. In our model, obvious basins do not extend farther north or west, although the region is modestly lower (~100 m) than its surrounding. This contrasts with the lower-resolution model of Tajeddine et al. (2017), which identified a small, shallow basin extending farther west.

Long-wavelength topography
Our model also resolves a deep and extensive basin in the northern portion of the trailing hemisphere (300 o W, 28 o N) that extends west to the sub-Saturnian hemisphere (0 o W and 50 o N). The quality of the shape model is relatively poor in this region, with sparse data and limited stereo (section 4); however, the broad outline of the basin is easily distinguished. The basin in our model is ~700 m deep. The region was identified as two separate basins by Schenk and McKinnon (2009), but this is likely due to a data gap in their model (where our data are also poor). Interestingly, the model of Tajeddine et al. (2017) also resolves a single deep basin centered at 330 o W (farther west relative to our current model and that of Schenk and McKinnon (2009)) with much smaller longitudinal extent. Figure 6: The gridded topography of Enceladus (values relative to the best-fit ellipsoid). Panels are as in Fig. 5.
The long-wavelength topography of Enceladus' leading hemisphere is also notable. The northwestern half consists of a broad depression 250-300 km in diameter and ~600 m deep, whereas the southeastern half consists of a broad high that sits ~800-1000 m above the reference. This large-scale topographic structure is not resolved in the Tajeddine et al. (2017) model, which shows a relatively flat region with a few small topographic depressions and highs. The region was not covered in the highresolution stereo model of Schenk and McKinnon (2009).
In some regions, our shape model clearly resolves smaller-scale topography associated with actual surface features. Relatively small craters (~7 km in diameter) are easily distinguished near 180 o W, 0 o N, even without the underlaying base image. We also resolve high-standing crater rims in these regions. Farther north, the rim of the large crater Dunyazad stands out above the depression to its west and south. Topography around the south polar terrain is also clearly associated with features, although the relationship is more complex. The band of ridges at 105 o W is high standing, as is Cashmere Sulci: the pronounced set of ridges south of the deep trough Labtayt Sulci. Labtayt Sulci itself is also resolved, although the model is noisy in the region. Between those two sets of ridges, the shape model has a clear artefact, as indicated by an abrupt change in elevation across what appears to be an image seam. This type of artefact is rare in the model, and the root cause is not yet clear. Polar projections of the south pole reveal that the model resolves (barely) the large fractures associated with plume eruptions. The observations of such smallscale features reveal the power of this new data set.

EVALUATING SHAPE MODEL QUALITY
Generating foundational data sets is necessary for the construction of a PSDI; however, a measure of the quality of that data set is also necessary to help ensure that data are used appropriately. Below we describe measures of the quality of the shape model. Given the highly variable image data set from which it was derived (section 2.2), shape model quality is also highly spatially variable.

Data density and variability
One of the simplest measures of data quality is data density: literally how many points were included in each pixel. As discussed in section 2.2, image coverage of Enceladus is highly variable, resulting in highly variable point density. Figure 7a shows a map of the point density of the model, which ranges from 162 points per 2.2-km pixel to 0 points per pixel. Point density is highest near the sub-and anti-Saturnian points, especially in the equatorial region near 210 o W where typical point counts are 40-60. This region includes both a few relatively high-resolution images and numerous, overlapping lower-resolution images (Fig.  2). Point density outside of these regions (dark grey) are typically much lower: often just 5-10 points. The mean point density for the model is 10.6 points per pixel (mean is 7 points per pixel). Point density is a useful measure of where the data are but is of limited use in understanding the actual quality of the data. Pixels with high point density might still be low quality if, for example, all the images have similar viewing geometry (low stereo strength). Alternatively, pixels with low point density may be high-quality if all the points were derived from excellent stereo views.
To better assess the variability of the points within each pixel we calculated the normalized median absolute deviation (NMAD) (Fig. 7b), which is a robust measure of the statistical dispersion of the data (i.e., more resilient to outliers than, e.g., standard deviation). Yellow pixels indicate that the NMAD is less than the precision of our bundle solution (57 m). Blue pixels indicate NMAD > 100 m, and thus highlight regions where point measurements are more variable. We note that some of this variability is "real": NMAD values are high near resolved craters (e.g., 200 o W, 3 o N), which results from actual topographic variation. NMAD values are typically less than 250 m, except for small isolated regions with limited data coverage.

Evaluating stereo strength
In general, stereo imaging was not a primary driver during the acquisition of Cassini ISS images. The availability of "good" stereo is therefore as much a function of happenstance as planning. In order to evaluate the spatial variability of the stereo strength of our model, we calculated the expected vertical precision (EP) of every pairwise combination of images that contribute to a given point in the model and report the minimum EP as a characteristic value for each point. We follow Kirk et al. (2003) and define EP = ρS/(p/h), where S is the RMS pixel scale of the two images, (p/h) is the parallax to height ratio of the images, and ρ is the quality of the matching in fractions of a pixel. We conservatively take ρ to be 1 (pixel, rather than subpixel matching) due to the often large disparity in pixel scale between the stereo images for each point. The result is that the EP reported here is likely an upper limit. The resulting minimum EP is shown in Fig. 8. The best 10% of the data have EP between 29 and 111 m, and 70% of the data have EP < 200 m. EP is highest at the poles where similar viewing geometry (south pole) and lack of data (north pole) severely limits stereo quality.

CONCLUSION
We have generated a publicly available, global, high-resolution (relative to existing products) shape model for Saturn's moon Enceladus, which is the target of ongoing scientific research and future exploration. The spatially variable quality of the model has been well-characterized. The model provides the final foundational data product necessary for an Enceladus PSDI and demonstrates the applicability of the technique to future outer Solar System missions, such as NASA's Europa Clipper and the European Space Agency's Jupiter Icy moons Explorer (JUICE) mission.