Determining the Optimal Grid Resolution for Topographic Analysis on an Airborne Lidar Dataset

. Digital Elevation Models (DEMs) are a gridded representation of the surface of the earth and typically contain uncertainties due to data collection and processing. The topographic metrics slope and aspect contain errors and uncertainties inherited both from the representation of a continuous surface as a grid (referred to as truncation error, TE), and from any DEM uncertainty. We analyze in detail the impacts of TE and propagated elevation uncertainty (PEU) on slope and aspect. Using synthetic data as a control, we deﬁne functions to quantify both TE and PEU for arbitrary grids. We then develop a 5 quality metric which captures the combined impact of both TE and PEU on the calculation of topographic metrics. Our quality metric allows us to examine the spatial patterns of error and uncertainty in topographic metrics, and to compare calculations on DEMs of different sizes and accuracies. Using lidar data with point density of ∼ 10 pts/m 2 covering Santa Cruz Island in southern California, we are able to generate DEMs and uncertainty estimates at several grid resolutions. Slope (aspect) errors on the one-meter dataset are on average 10 0.3 ◦ (0.9 ◦ ) from TE, and 5.5 ◦ (14.5 ◦ ) from PEU. We calculate an optimal DEM resolution for our SCI lidar dataset of four meters that minimizes the error bounds on topographic metric calculations due to the combined inﬂuence of TE and PEU for both slope and aspect calculations over the entire SCI. Average slope (aspect) errors from the four meter DEM are 0.25 ◦ (0.75 ◦ ) from TE and 5 ◦ (12.5 ◦ ) from PEU. While the smallest grid resolution possible from the high-density SCI lidar is not necessarily optimal for calculating topographic metrics, high point-density data are essential for measuring DEM uncertainty 15 across a range of resolutions.


Introduction
Continuous surfaces are often projected onto an evenly sampled grid -digital elevation models (DEMs) are a common example.
The accuracy of the gridded representation of the underlying surface is controlled by the spacing of the grid, the variability of the surface (i.e., the terrain itself), and the amount of uncertainty added during data collection and processing. 20 In recent years, an ever-growing array of DEM datasets have become available across a range of resolutions and spatial scales. As new data acquisition and processing strategies -such as lidar, stereo photogrammetry, and structure from motion -mature, the number and variability of features represented in a DEM will grow dramatically. In this manuscript, we refer to DEMs as high resolution when their grid spacing is low (e.g., a one meter DEM), and low resolution when their grid spacing is high (e.g., a 30 meter DEM). Across all DEM resolutions, it is important to quantify the uncertainties in DEMs, and their 25 (e.g., Roering et al., 1999;Dietrich et al., 2003;Pelletier, 2008;Grieve et al., 2016), hydrologic applications (e.g., Band, 1986;Zhang and Montgomery, 1994;Tarboton, 1997;Pelletier, 2010Pelletier, , 2013 including flood risk modeling (e.g., Ouma and Tateishi, 2014), tectono-geomorphic modeling (e.g., Whipple and Tucker, 1999;Snyder et al., 2000;Tucker and Hancock, 2010;Kirby and Whipple, 2012;Bookhagen and Strecker, 2012), ecologic modeling (e.g., Franklin, 1995;Guisan and Zimmermann, 2000;Thompson et al., 2001), and vegetation analysis (e.g., Pierce et al., 2005;Kent, 2011).

10
A wide range of methods have been developed to accurately derive slope and aspect from elevation data using a range of approaches, each optimized for different use cases (e.g., Evans, 1980;Zevenbergen and Thorne, 1987;Skidmore, 1989;Tarboton, 1997;Dunn and Hickey, 1998;Schmidt et al., 2003). For example, the methods D-8 and D-Infinity are optimized towards hydrologic flow-routing problems, particularly on DEMs of low spatial resolution (Tarboton, 1997), and the methods of Evans (1980) and Zevenbergen and Thorne (1987) include some smoothing of the underlying DEM to minimize the impacts 15 of DEM noise. However, all of these methods contain some intrinsic error due to truncation error (TE) -the deviation of the gridded sample from the continuous original surface (Figure 1). The TE describes the error made by truncating an infinite sum and approximating it by a finite sum. It often includes a discretization error, which arises from taking a finite number of steps to approximate an continuous surface. Several authors have explored the theoretical limitations of estimations of slope and aspect calculations on gridded surfaces (e.g., Heuvelink et al., 1989;Hunter and Goodchild, 1997;Florinsky, 1998;Abarbanel 20 et al., 2000;Zhou and Liu, 2004); in general, TE can be related directly to the grid-sampling -more tightly sampled grids will deviate less from the original surface. In addition to TE, real-world DEMs will have some degree of measurement uncertainty. The magnitude of that uncertainty varies across data collection methods and data post-processing, but also depends on the terrain itself and is often difficult to estimate (Heuvelink et al., 1989;Lee et al., 1992;Bolstad and Stowe, 1994;Florinsky, 1998;Fisher, 1999;Smith and Sandwell, 2003;Farr et al., 2007;Wechsler and Kroll, 2006;Wechsler, 2007;Mukul et al., 2017;Bookhagen, 2017, 2018;Wessel et al., 2018). Any uncertainty in elevation estimates will propagate into the calculated topographic metrics.
The impacts of TE and DEM uncertainty on slope and aspect estimations can be quantified for any gridded data. We first 5 use synthetic data with known properties as a control to define generalized functions applicable to any DEM, and to develop a quality metric for slope and aspect calculations. Following this analysis, we turn to a high-resolution lidar dataset covering complex terrain to study the spatial structure of uncertainty in topographic metrics propagating from both TE and DEM uncertainty. This novel approach allows us to identify the optimal grid resolution that minimizes the error bounds from the combination of TE and PEU, and to analyze the implications of using sub-optimal DEM resolutions for calculating slope and 10 aspect.

Data and Methods
In this analysis, we demonstrate the limitations of calculating slope and aspect using synthetic data surfaces. These surfaces serve as a valuable control dataset, as the precise analytical values for slope and aspect at each grid cell are known. We further add normally distributed random noise with a known mean and standard deviation to our synthetic data to investigate the 5 impacts of DEM uncertainty on topographic metrics. Following the discussion of synthetic data, we apply the same methods to high-resolution lidar data interpolated to DEMs with varying spatial resolutions.

Synthetic Data
We generate regular n × n grids of size n=11, 101, 1001 describing a Gaussian hill with a maximum height of one, with elevations defined as z = e (−x 2 −y 2 ) ( Figure 2). A second surface describing a sphere is included in the Supplement. It is 10 important to note that the horizontal and vertical units for the synthetic data are arbitrary -the grid cell spacings generated by using differently sized grids are an expression of the width-to-height ratio of the surface. The python functions used to generate our synthetic surfaces can be found on GitHub: https://github.com/UP-RS-ESP/TopoMetricUncertainty.

Topographic Metrics
In this study, we focus on the topographic metrics slope S(x, y) and aspect A(x, y) which we define as, 15 S( with x and y being the spatial coordinates. Key in calculating slope and aspect is the calculation of the directional derivatives ∂z ∂x and ∂z ∂y . There exist a wide range of methods for calculating the directional derivatives on gridded data that broadly fall into three classes: (1) four-neighborhood methods, (2) eight-neighborhood methods, and (3) steepest descent methods (e.g. Evans, 1980;5 Zevenbergen and Thorne, 1987;Skidmore, 1989;Tarboton, 1997;Dunn and Hickey, 1998;Schmidt et al., 2003;Zhou and Liu, 2004;Oksanen and Sarjakoski, 2005). In our analysis, we use the standard second-order finite difference approximation as implemented in Python Numpy (Fornberg, 1988;Durran, 1999), which is a four-neighborhood method. We tested additional methods using eight-neighborhoods (e.g., singular-value decomposition plane fitting and the method of Evans (1980)) -which provide some underlying DEM-smoothing and are well-suited to noisy DEM data -and using steepest descent (e.g., D-8 10 and D-Infinity (Tarboton, 1997)), which identify the steepest slope that water would take and are optimized for hydrologic use cases. However, these methods do not necessarily provide the most representative slope of the terrain, and show distinct differences from the analytical solutions of slope and aspect (cf. Supplemental Figures 2 and 3). We choose to use the simple four-neighbor method in our analysis as it results in the lowest relative error when compared with the analytical solution for slope and aspect (Supplemental Figure 2) and does not include any smoothing of the underlying DEM. 15

Sources of Uncertainty in Topographic Metrics
Biases in topographic metrics can be determined both numerically, by comparing the results of calculations to known values, and analytically, by deriving the impacts of both TE and DEM uncertainty on topographic metrics.

Truncation Error
We derive the TE from the formulas of topographic metrics by propagating the uncertainty of the second-order finite difference 20 approximation for the directional derivatives ∂z ∂x and ∂z ∂y . Neglecting rounding errors, we assume the uncertainty ε of the second-order finite difference to be bounded by where a represents either x or y, and h is the grid spacing. The considered topographic metrics slope S(x, y) and aspect A(x, y) depend on both directional derivatives ∂z ∂x and ∂z ∂y . As 25 the corresponding uncertainties ε x and ε y are not independent, propagated TE T for slope (T S ) and aspect (T A ) are given by These can be more explicitly expressed using Equation (3) as The sign of the TE is derived from symmetry considerations in the following way:   As can be seen in Figure 3, the magnitude of the slope and aspect error decreases with increasing n. This implies that as more tightly sampled models of a surface are used, the accuracy of slope and aspect estimations increases. This follows the work of Abarbanel et al. (2000), which shows that finite-difference derivative estimations have error magnitudes that scale with the inverse of the grid spacing for evenly spaced grids. However, even with increasingly accurate surface models, there remains some TE simply due to the fact that a smooth surface cannot be perfectly represented by gridded data. At very fine scales, computer rounding will also lead to errors in slope and aspect calculations. In this study, however, we do not use fine 5 enough grids for this to impact our error estimates.
As the precise mathematical definitions of our synthetic surfaces are known, we can analytically derive their slope and aspect values at each point and compare those results to numerical calculations (see Supplemental Figure 4). In the absence of DEM uncertainty, the difference between analytically-and numerically-derived topographic metrics will be dominated by the TE from the second-order finite difference approximation. Additionally, modifying the absolute magnitude of the shape heights (e.g., a maximum height of 100 instead of one), does not modify the absolute magnitudes of slope and aspect TE (see Supplemental Figure 5).
While the magnitudes of the errors are small, they have important impacts on distributions drawn from larger grids (e.g., entire surfaces, see Figure 4). In the case of aspect, the alternating spatial pattern of its TE pushes aspect values towards the cardinal directions -leading to distinct positive and negative spikes in the aspect distribution. This error is not confined only to the cardinal directions -there exists symmetry across the aspect distribution that is easiest to see with very fine binning (see  Slope and aspect calculations have very differently spatial error patterns (see Figure 3). However, as with aspect errors, slope errors depend on the orientation of the grid despite the radial symmetry of the surface. This is easiest to see in the slope difference map of the Gaussian hill, where error magnitudes are not simply scaled by slope (see Figure 3E).

10
Elevation models are never perfect -there are always errors and uncertainties due to data collection or processing. The uncertainty in calculated topographic metrics can be constrained by propagating DEM uncertainties into their calculations (e.g. Heuvelink et al., 1989;Hunter and Goodchild, 1997;Zhou and Liu, 2004;Oksanen and Sarjakoski, 2005); we refer to the uncertainty introduced into slope and aspect calculations from elevation uncertainty as propagated elevation uncertainty (PEU).
For our synthetic data, we assume the simplistic but straight forward error model of spatially uncorrelated white noise. We de-15 fine noise to be normally distributed and consider two cases: (1) homogeneous white noise, i.e., noise drawn from one unique normal distribution N (µ = 0, σ 2 = const), and (2) slope dependent, spatially varying white noise, i.e., noise drawn from a family of normal distributions N (µ = 0, σ 2 S ). Since for each grid cell noise is drawn independently, the PEU for slope (E S ) and aspect (E A ) is given by As before, this can be expressed more explicitly using Equation (3) as As with TE, PEU introduces distinctive spatial patterns into the slope and aspect estimates ( Figure 5). Both slope and aspect 5 uncertainties are higher on shallow slopes, as has been shown in previous work (Florinsky, 1998;Zhou and Liu, 2004). On flat or nearly flat terrain, the elevation uncertainty will more strongly impact the estimated partial derivatives, and thus lead to erroneous slope and aspect estimates. It is also important to note that aspect varies between [0, 360) and slope only between [0, 90), which is why we expect standard deviations to be different by a factor of four (see Figure 5). Slope STD (degree) Figure 5. Gaussian hill (n=101) slope and aspect standard deviations for spatially invariant noise (std=1e −4 ), derived from elevation uncertainty propagation (Equations 12-13). Aspect has higher standard deviations than slope, as aspect varies between [0, 360) and slope only between [0, 90). This is evident especially on flat terrain (e.g., the edges of the Gaussian hill).
When we compare our analytically-derived topographic metric uncertainty patterns ( Figure 5) to those generated from actual 10 random noise -using an ensemble of 10,000 synthetic surfaces, each generated with normally distributed noise -the spatial patterns and magnitudes are nearly identical between the two methods (cf. Supplemental Figure 7). Both TE and PEU yield distinct spatial patterns, implying that the error and uncertainty in topographic metrics will vary throughout a DEM. Using both the TE and PEU magnitudes, we can examine whether TE or PEU is dominant at an arbitrary grid spacing, given a known DEM uncertainty ( Figure 6). As can be seen in Figure 6, TE will approach zero for sufficiently small grid spacings, while PEU will rise as the amount of noise relative to the grid spacing increases. Using both TE and PEU magnitudes across grid spacings, an optimal grid spacing 5 for a given noise level which minimizes their combined influence can be chosen. Note that we do not assume that the quality of the data is maximized when the combination of TE and PEU is minimized, but that the error bounds on topographic metric calculations are minimized. Thus, the optimal grid spacing is preferred for calculating slopes and aspects for use in applications where uncertainty across the entire grid should be minimized.

10
The trade off between TE and PEU can also be thought of as a quality ratio. Simply put, the impact of DEM uncertainty on topographic metrics is modulated by the grid spacing -i.e., 1 cm of vertical noise on a 10 m grid will have a very different impact than 1 cm of vertical noise on a 1 cm grid. As both TE and PEU have distinct spatial patterns, a combined Quality Ratio where m is a normalization factor which accounts for the four-fold difference in the range of possible values between slope and aspect. On a surface without noise, the QR will approach one as grid spacing and TE approach zero. The influence of TE and PEU on calculated QRs can be seen by comparing the same synthetic grid with two noise levels (see Supplemental Figure   5 8). High noise levels reveal similar spatial patterns to those seen in Figure 5, while low noise levels show the influence of TE (cf. Figure 3). The QR can also used as a normalized metric to compare DEMs across grid spacings and noise levels ( Figure 7, and Supplemental Figure 9). Each noise level/grid spacing pair has a point at which the QR of the DEM is maximized. This point represents the optimal trade off between TE and PEU over the entire grid.
As can be seen in Figure 7, as noise levels increase, the optimal grid spacing increases, and the QR of that optimal spacing decreases. For very noisy datasets (purple and yellow lines, Figure 7), very large grid cell spacings -i.e., aggregation -can help reduce the combined influence of TE and PEU on slope and aspect estimates.

5
There are distinct differences between the optimal slope and aspect grid spacings for the same dataset. For any given noise level, the optimal grid spacing to calculate aspect is higher than that for slope. The difference in optimal grid spacing is driven by the relative magnitudes of the PEU between slope and aspect calculations -the PEU for aspect calculations is always higher than that for slope due to differences in the formulas for slope and aspect (see Figure 6 and Equations 12-13).

Heterogeneous Noise
In practice, datasets will not have homogeneous noise across grid spacings -coarser parameterizations of a surface will include more uncertainty as fine-scale features are aggregated into single pixels. Additionally, noise in real-world datasets is influenced by non-uniform landscape features such as slope, aspect, terrain relief, and vegetation cover (e.g., Kraus and Pfeifer, 1998;Kyriakidis et al., 1999;Holmes et al., 2000;Smith and Sandwell, 2003;Carlisle, 2005;Oksanen and Sarjakoski, 2006;Rodriguez et al., 2006;Shortridge and Messina, 2011;Bookhagen, 2017, 2018). While several authors (e.g., Hunter 15 and Goodchild, 1997; Zhang and Goodchild, 2002;Fisher and Tate, 2006) have examined complex error models in real data, we choose to use a simplistic error model based on terrain slope -e.g., higher slopes have higher noise levels -to look at the first-order impacts of heterogeneous noise on the calculation of topographic metrics.
When we compare the heterogeneous noise result to the analysis shown in Figure 7 (cf. Supplemental Figure 10), the results are very similar, albeit with higher optimal grid spacings. We can also compare the optimal grid spacings for homogeneous 20 and heterogeneous noise directly for a set of noise levels from 1e −8 to 1e −2 ( Figure 8).
As can be seen from Figure 8, slope can be calculated on a smaller grid spacing with slope-biased noise, while aspect should be calculated with a larger grid spacing. The divergence in the optimal grid spacings between homogeneous and heterogeneous noise are driven by the differences in the influence of terrain slope on slope and aspect QRs. While both slope and aspect calculations have similar grid-spacing dependent TE magnitudes, the differences in the relative magnitudes of the PEU lead 25 to different responses to the distribution of noise on the DEM (see Figure 6). The simplistic slope-dependent error model used here is unlikely to translate to a real-world setting; however, it is clear that the spatial distribution of noise on the DEM influences the choice of optimal grid spacing for both slope and aspect calculations.
Heterogeneous noise leads to smaller optimal grid spacing Homogeneous noise leads to smaller optimal grid spacing On the line, the noise distribution does not change the optimal grid spacing A B Figure 8. Optimal grid spacing for slope (A) and aspect (B) across a range of noise levels (1e −2 to 1e −8 ). Points in the upper triangle have smaller optimal grid spacings with homogeneous noise. Points in the lower triangle have smaller optimal grid spacings with heterogeneous noise. Points on the dashed line have the same optimal grid spacings regardless of the noise distribution. Heterogeneous noise leads to smaller slope grid spacings but larger aspect spacings for the same average noise level.

Impacts of Noise on Topographic Distributions
As can be seen in Figures 5-7, the introduction of noise has non-trivial impacts on slope and aspect calculations. As was evident in Figures 1 and 4, TE can significantly modify the ensemble distributions of topographic metrics. To examine the impact of both TE and PEU on the frequency distributions of slope and aspect, we generate 1,000 member ensembles of slope and aspect for each combination of two different grid sizes (n=101, 1001) and three different noise levels (homogeneous noise, std=1e −2 , 1e −3 , 1e −4 ). We then bin the ensemble data into one-degree slope and aspect bins. When homogeneous noise is added to synthetic surfaces, the large spikes in aspect distributions related to points moving in and out of cardinal directions (e.g., at 45/90/135 degrees, see Figure 4) are smoothed out into neighboring bins. However, this does not imply that the calculated aspects are 'better', just that they have a more even distribution -the PEU overprints the TE as seen in Figure 3. This can clearly be seen in Supplemental Figure 11, which shows the impact of noise on slope and aspect grids.
Slope distributions are systematically shifted towards higher values as more noise is added to the synthetic surfaces. This 5 effect is due to the presence of more large 'steps' in the synthetic data, where smooth transitions across elevation gradients are replaced with more stepped hills and dips. Across noise levels and grid sizes, distributions with similar QRs maintain the same general shape. This indicates that the QR captures how 'wrong' the aggregate distribution is when compared to the original synthetic surface.
If different window sizes are used to calculate slope and aspect (e.g., 5x5 and 7x7), slope estimates are higher and the 10 aspect 'spikes' at cardinal directions are slightly diminished. Larger window sizes somewhat -but not completely -remove the impacts of noise on aspect distributions. They do not, however, help minimize differences in the slope distributions (see Supplement Figure 12).
6 Case Study: Multi-Resolution Lidar on Santa Cruz Island

Dataset Description
The Santa Cruz Island (SCI, see Figure 10) is part of the Channel Islands off the coast of south-central California. It is a tectonically active island (Neely et al., 2017) with steep topography and pronounced erosion (Perroy et al., 2010(Perroy et al., , 2012, diverse lithologic cover (Dibblee, 2001), and moderate and species-rich vegetation cover (Baguskas et al., 2014). The lidar point cloud data used to derive DEMs at a range of spatial resolutions were obtained from the 2010 US Geological Survey

Deriving Elevations and Uncertainties from Point Cloud Data
Using LAStools (LAStools, 2017), we use a triangulated irregular network approach to derive average point cloud elevations at each grid cell center, for each of our chosen grid resolutions (1-30 meters, in one meter steps, and 90 meters). We have tested simpler interpolation schemes, including the lowest elevation and the central location for each grid cell, but the resulting uncertainties in topographic metrics remain the same. Thus, we decided to use the more commonly applied triangulated irregular network approach. In addition to producing DEMs, we generate pixel-wise standard deviations estimated from the point 5 cloud. That is, we determined the standard deviation of all ground-classified lidar points falling into each grid cell for each grid resolution. While vegetation cover may impact ground-classification results, our calculations are based on ground-classified points only. Because of the high point-cloud density, this standard deviation metric should reflect terrain variability and not vegetation cover. This measure serves as a DEM uncertainty and is highly spatially heterogeneous. Elevation and standard deviation maps for additional grid resolutions can be found in the Supplement (Supplemental Figures 14-16). An island-wide 10 canopy-height model can be found in Supplemental Figure 13.

Mapping Uncertainty on SCI
Using our standard deviation grids, we can then calculate both TE and PEU on SCI (see Supplemental Figure 17). Using these two grids as inputs, we can also derive the QR across the entire SCI ( Figure 11). Note that maximal QRs are low when compared to the theoretical models shown in Figures 7 and 9; this is due to both the much higher inherent noise in the real-world data, 15 and to the more complex terrain represented by the lidar data.

Identification of Optimal Grid Resolution
Using our pixel-wise standard deviations -and derived QRs -at each grid resolution, we can determine the island-wide optimal grid resolution for the calculation of topographic metrics which minimizes the combined influence of TE and PEU on our lidar dataset.
Unlike the synthetic data examined in Figure 7, the optimal grid resolution for aspect calculations is the same as that for slope. This is due to the distribution of both elevations and elevation uncertainties in the SCI DEM, and the different 5 responses of slope and aspect calculations to the diverse terrain on SCI. Whole-island median TE, PEU, and QR can be found in Supplemental Table 1 for grid spacings from one to ten meters. While four meters is the calculated optimal grid resolution, the QR peak is quite flat (see Figure 12), implying that both three-meter and five-meter DEMs will yield slope and aspect results of nearly equal quality. Indeed, the TE and PEU magnitudes are similar for both metrics between three and five meters (see Supplemental Table 1). It should be noted that a four-meter optimal grid resolution is likely too large for many applications, particularly those interested in micro-topography or other small-scale features. In those cases, a higher-resolution DEM should be used with caution -PEU will continue to grow as grid spacing shrinks, and the rate of that increase will accelerate (see Figure 6). This effect will occur even if the DEM uncertainty decreases at smaller scales simply due to the larger impacts of PEU at very small grid spacings. This means that the uncertainty bounds on slope and aspect estimates on very high resolution DEMs will increase quickly, which may limit some analyses. For high-resolution analyses, very high quality DEMs are required.

5
In this analysis of the airborne lidar dataset from SCI, we find that a four-meter DEM resolution minimizes the combined influence of TE and PEU. Slope and aspect analyses performed at this grid resolution have the smallest error bounds, and thus provide the most reliable data for analysis of topographic metrics over the entire island.
While a single, whole-island, grid resolution is useful in some applications, it is clear from Figure 11 that there are large spatial variations in the DEM uncertainty and the resulting QR, and hence in the optimal grid resolution. We can extend our 10 method of identifying the optimal grid resolution on a catchment-by-catchment basis ( Figure 13, map view in Supplemental Figure 21).
There exist large spatial variations in the optimal grid resolution across SCI, driven by differences in terrain slope and DEM quality. When these differences are compared to median catchment slope and median catchment standard deviation, a clear pattern emerges (Figure 13). While there remains significant scatter in the data, catchments with high median slopes tend to 15 have both higher standard deviations and lower optimal grid resolutions. The data is bimodal, with the majority (83% for slope, 95% for aspect) of catchments having optimal resolutions below ten meters. resolution, colored by the standard deviation of that point. Stars show the whole-island optimal resolution and median slope. There exists a clear trend where lower median slopes lead to lower optimal grid resolutions. Standard deviations also tend to be smaller in these catchments, due to higher DEM fidelity over flat terrain. Catchments with high (>10 meter) optimal grid spacings are small and have complex topography.
For a map view of these results, we refer to Supplemental Figure 21.
It is worth noting that the vast majority of the catchments with optimal slopes and aspects above ten meters are small catchments concentrated on the northwestern edge of SCI (see Supplemental Figure 21). These catchments have two unique aggregating several small catchments into one larger analysis unit (to improve the statistical reliability of the median QR) result in smaller optimal grid spacings. We posit that the high optimal grid spacings are a product of the unique two-stage topography in this region. However, an in-depth analysis of the specific factors which drive the optimal terrain resolution in each individual 5 zone is outside of the scope of this manuscript. In applying this method of deriving optimal grid spacings, it is important to take into account the possibility of large differences in the optimal grid spacing over relatively small spatial scales.

Island-Wide Slope and Aspect Distributions
As has been shown in Figures 7 and 12, there is clearly an optimal spatial resolution for the calculation of slope and aspect on our SCI lidar dataset. Furthermore, in Figure 9, it was shown that changes in QR can have a substantial impact on the slope 10 and aspect distributions of a synthetic surface. This analysis can be extended to the SCI lidar data at a range of spatial scales ( Figure 14).  Figure 11 and Supplemental Figures 13-15). Optimal slope and aspect distributions plotted in black. (C) Quantile-quantile plot showing difference between the entire SCI slope distribution vs the 'optimal' (four meter) resolution. We observe a clear pattern in slope distributions where resolutions below the optimal resolution tend to overestimate slope, those below tend to underestimate slope.
There are clear differences in the slope and aspect distributions with increasingly poor QRs. This effect can also be seen in quantile-quantile plots of the entire SCI slope distribution when compared to the optimal grid resolution (see Figure 14C). Grid resolutions higher than the optimal result in slope distributions that are dominated by higher slopes, while lower resolutions 15 tend towards lower slopes. This effect can be seen even when the grid resolution is only slightly removed from the optimal resolution (Supplemental Figure 22). Notably, grid resolutions close to the optimal (e.g., three meters and five meters), have QRs and slope distributions that are very similar to the optimal four meter distribution.
As a final test of grid-resolution dependence, we resampled several grids to nominally match the optimal grid resolution -in this case we chose three meters for the Pozo catchment (see Figure 10, Supplemental Figure 23) as the reference distribution.

5
As can be seen in Figure 15A, there are differences between test (one to five meter) resolutions and the reference three meter histogram. When the test datasets are resampled to the three meter reference resolution (using bilinear resampling), the histograms do not converge to the reference three meter distribution ( Figure 15B). Slope distributions from resampled data should be used with caution; simply resampling elevation data to the optimal grid resolution before calculating slope will not yield better results -the DEM has to be created from the ungridded point data to fully take advantage of optimal gridding for slope  Ideal Slope Distribution 3m --QR: 0.141 Figure 15. (A) Divergence of higher and lower grid-resolution slope distributions from that at the optimal resolution for the Pozo catchment (see Figure 10). Excepting the one meter histogram, the differences between resolutions are small. (B) Slope distributions taken from one, two, four, and five meter data resampled to three meter resolution using bilinear resampling. The resampled data distributions do not converge on the reference three meter distribution. This study presents a detailed accounting of uncertainties and errors in the topographic metrics slope and aspect derived from both truncation errors and uncertainty in the underlying source digital elevation model (DEM). We first develop our analysis on synthetic data, which acts as a control dataset. We then extend our methodology onto a high point-density lidar dataset, which allows us to compare the relative impacts of gridding truncation errors and propagated elevation uncertainty on the calculation of topographic metrics.

5
From the analysis of both synthetic and real-world data, we identify the following key points: (1) the relative impact of truncation error and propagated elevation uncertainty can be captured in a single quality metric, which we coin the Quality Ratio. This metric can be used to compare the accuracy of topographic metrics across DEM spatial resolutions and uncertainty distributions.
(2) There exists an optimal grid resolution at which to calculate the topographic metrics slope and aspect for a given dataset that minimizes the total impact of both truncation error and propagated elevation uncertainty; the distribution of 10 DEM uncertainties leads to spatial variation in the optimal grid resolutions at which to calculate slope and aspect. For the Santa Cruz Island in southern California, we find an optimal grid resolution of four meters, with island-average slope (aspect) errors of 0.25 • (0.75 • ) for truncation error, and 5 • (12.5 • ) for propagated elevation uncertainty.
(3) Topographic metrics calculated at sub-optimal grid resolutions will have potentially large errors in their slope and aspect distributions. Resampling grids to the optimal resolution does not completely ameliorate these issues; it is important to generate the DEM at the optimal resolution 15 from the underlying lidar dataset. (4) Slope and aspect calculations at high DEM resolutions (>3 meters) are significantly impacted by the nonlinear relationship between grid spacing and propagated elevation uncertainty. A DEM with very low uncertainties is required to support robust analysis of fine-scale topography.
Given that grid-resolution-driven effects on regional slope and aspect distributions could have significant impacts on the interpretation of landscape morphology, we recommend that region-specific optimal DEM resolutions be determined before 20 the calculation of topographic metrics.
Code and data availability. All codes and data for reproducing the results can be found in on GitHub