Amante, C.J., 2018. Estimating coastal digital elevation model uncertainty.
Integrated bathymetric–topographic digital elevation models (DEMs) are representations of Earth's solid surface that extend across the coastal land–water interface. DEMs are fundamental to the modeling of coastal processes, including tsunami, storm-surge, and sea-level-rise inundation. Vertical errors in coastal DEMs are deviations in elevation values from the actual seabed or land surface, which originate from the (1) elevation measurements, (2) datum transformation that converts bathymetric and topographic measurements to a common vertical reference system, (3) spatial resolution of the DEM, and (4) interpolative gridding technique that estimates elevations in areas unconstrained by measurements. The magnitude and spatial distribution of the vertical errors are typically unknown, and a DEM uncertainty surface is a statistical assessment of the likely magnitude of these errors. The National Oceanic and Atmospheric Administration National Centers for Environmental Information develops DEMs for United States' coastal communities. This study describes a methodology to derive uncertainty surfaces that estimate coastal DEM vertical errors at the DEM cell–level. A coastal DEM south of Sarasota, Florida is the case study for deriving uncertainty surfaces. Results indicate that large vertical uncertainty exists in deeper waters offshore with sparse echo-sounder measurements, and small vertical uncertainty exists on flat terrains with dense light detection and ranging measurements. The estimated uncertainty can be propagated into the modeling of coastal processes that utilize DEMs by deriving numerous plausible DEM realizations within the uncertainty bounds. The numerous DEMs realizations can then produce an ensemble of coastal modeling results, and in turn, better-informed coastal management decisions.
INTRODUCTION
Integrated bathymetric–topographic digital elevation models (DEMs) are representations of Earth's solid surface that extend across the coastal land–water interface by seamlessly merging subaerial topography with adjacent bathymetry (Danielson et al., 2016; Eakins and Grothe, 2014; Gesch and Wilson, 2001; Thatcher et al., 2016). The National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (NCEI) develops DEMs for United States' coastal communities to support numerous coastal modeling efforts, including the modeling of tsunami propagation and coastal inundation (Eakins and Taylor, 2010).
Vertical errors in DEMs are defined in this manuscript as deviations in elevation values from the actual seabed or land surface (Hunter and Goodchild, 1997; Li et al., 2018). Such vertical errors originate from numerous sources, including the (1) elevation measurements (e.g., sonar, light detection and ranging [LIDAR]), (2) datum transformation that converts bathymetric and topographic measurements to a common vertical reference system, (3) spatial resolution of the DEM, and (4) interpolative gridding technique (e.g., spline, kriging) that estimates elevations in areas unconstrained by measurements. The magnitude and spatial distribution of DEM vertical errors are typically unknown. DEM uncertainty represents the lack of knowledge of the vertical errors, and a DEM uncertainty surface is a statistical assessment of the likely magnitude and spatial distribution of these errors (Hunter and Goodchild, 1997; International Hydrographic Organization, 2008; Li et al., 2018; Wechsler, 2007). Accuracy is defined in this manuscript as a general term for the agreement of values to known or accepted values (Amante and Eakins, 2016), and is typically assessed by statistical measures, such as root mean square error (RMSE).
DEM uncertainty affects the fidelity of coastal process modeling, such as tsunami propagation and coastal inundation (e.g., Gesch, 2013; Hare, Eakins, and Amante, 2011; Leon, Heuvelink, and Phinn, 2014). Consequently, it is important to estimate and incorporate DEM uncertainty in the modeling of coastal processes. The estimated uncertainty can be propagated into the modeling of coastal processes that utilize DEMs by deriving numerous plausible DEM realizations within the uncertainty bounds (e.g., Leon, Heuvelink, and Phinn, 2014). The numerous DEM realizations can then produce an ensemble of coastal modeling results, and in turn, better-informed coastal management decisions. Estimating the spatially varying DEM uncertainty also aids in prioritizing future elevation data collection, which will subsequently also improve the fidelity of coastal process modeling.
Numerous studies address uncertainty in topographic DEMs (e.g., Bater and Coops, 2009; Goulden et al., 2016; Leon, Heuvelink, and Phinn, 2014; Spaete et al., 2011; Su and Bork, 2006) and uncertainty in bathymetric DEMs (e.g., Amante and Eakins, 2016; Calder, 2006; Elmore et al., 2012; Jakobsson, Calder, and Mayer, 2002), but there is currently a lack of research on deriving uncertainty surfaces for coastal DEMs developed from multiple bathymetric and topographic data sets of disparate quality and age. Jakobsson, Calder, and Mayer (2002) and Calder (2006) create uncertainty surfaces that reflect potential measurement uncertainty for bathymetric data sets from different time periods, and the interpolation uncertainty between sparse measurements. However, these studies are limited in the context of coastal DEMs as they do not integrate bathymetric and topographic data sets into a seamless DEM, and do not consider the effect of the DEM spatial resolution on cell-level uncertainty. Zhang et al. (2015) improve the accuracy of coastal DEMs that integrate bathymetric and topographic data sets by incorporating the relative accuracy of three data sources to optimize the weighting of each data set in the interpolation process, but their research also does not consider the effect of the DEM spatial resolution on cell-level uncertainty. Furthermore, Zhang et al. (2015) do not derive an accompanying uncertainty surface. To the best of the author's knowledge, there is no published research on developing coastal DEMs from multiple topographic and bathymetric data sets of disparate quality and age with accompanying uncertainty surfaces that estimate potential DEM vertical errors at the cell-level that originate from the (1) elevation measurements, (2) vertical datum transformation, (3) DEM spatial resolution, and (4) interpolation technique.
Measurement Uncertainty
A primary contribution to DEM uncertainty is the uncertainty of the elevation measurements constraining the model. In the bathymetric realm, the International Hydrographic Organization's (IHO) standards for hydrographic surveys provide guidance on the allowable magnitude of depth measurement uncertainty that results from data collection and processing (Hare, Eakins, and Amante, 2011; International Hydrographic Organization, 2008). The IHO determines various orders of standards on the basis of the importance of under-keel clearance, with stricter standards, i.e. less allowable vertical uncertainty, in cases where under-keel clearance is critical (e.g., shipping lanes in shallow waters). All standards are provided as a function of depth, resulting in larger allowable uncertainty in deeper waters. Sources of probabilistic measurement uncertainty for modern single-beam and multibeam echo sounders originate from the platform, sensor, environment, integration, and calibration (Hare, Eakins, and Amante, 2011). Other technology, such as LIDAR utilizing blue-green wavelengths (∼532 nm), can penetrate shallow, clear water to measure depths near the land–water interface (Gao, 2009; Irish and Lillycrop, 1999). The uncertainty of bathymetric LIDAR is also typically estimated on the basis of the depth, with larger uncertainty in deeper waters, in accordance to IHO standards (Costa, Battista, and Pittman, 2009; International Hydrographic Organization, 2008).
Legacy bathymetric data sets have additional uncertainty including digitization, shoal biasing, and morphologic change, and estimating the uncertainty of legacy data sets creates additional challenges (Calder, 2006; Elmore et al., 2012; Hare, Eakins, and Amante, 2011; Jakobsson, Calder, and Mayer, 2002; Marks and Smith, 2008). Jakobsson, Calder, and Mayer (2002) estimate the vertical uncertainty of bathymetric data on the basis of the navigation system and depth measurement instrumentation (i.e. echo sounder) described in the metadata, and provide a worst-case scenario if the metadata are unavailable. Depth measurement uncertainty can also be estimated by comparison with presumed higher-accuracy data sets (Calder, 2006; Marks and Smith, 2008). Marks and Smith (2008) determine that the worst-case scenario (i.e. 95th percentile error) is approximately five times larger in pre-1969 sonar data than in post-1968 sonar data, and derive separate models to estimate depth measurement uncertainty on the basis of the depth and terrain slope for the two discrete time periods.
NOAA disseminates categorical zones of confidence (ZOC) with their nautical charts that are determined primarily by the age of the data that inform the chart depths (S. Legeer, personal communication). Some areas of nautical charts date back to the 19th century, in which the technology of the day (i.e. lead line surveys) results in large uncertainty in the chart depths. The ZOCs are also based on water depth, resulting in larger measurement uncertainty in deeper waters. More information on depth measurement uncertainty, including the IHO standards and NOAA ZOCs, is provided in Calder (2006), International Hydrographic Organization (2008), and Hare, Eakins, and Amante (2011).
In the topographic realm, LIDAR technology typically utilizes near-infrared wavelengths (∼1064 nm) to measure elevations of Earth's surface (Heritage and Large, 2009). Postcollection filtering is performed to remove LIDAR returns from vegetation and buildings, and the remaining ground returns are utilized in NOAA DEMs to represent the bare-earth conditions. The vertical accuracy of topographic LIDAR is often provided in the data set's metadata by a global statistic, such as RMSE. The RMSE represents the accuracy of a LIDAR data set containing millions of elevation measurements, but is commonly derived using a relatively small number (∼tens to hundreds) of colocated ground control points. The number of ground control points represents an extremely small percentage of the LIDAR data set, which brings into question the robustness of the accuracy assessment (Wechsler, 2007). Furthermore, a single, global metric of accuracy is limited as LIDAR accuracy is correlated with land cover and terrain (Bater and Coops, 2009; Goulden et al., 2016; Leon, Heuvelink, and Phinn, 2014; Spaete et al., 2011; Su and Bork, 2006). For example, LIDAR elevation measurements are typically biased toward higher elevations than the actual bare-earth surface in densely vegetated, coastal marshes because of poor laser pulse penetration (Hladik and Alber, 2012; Schmid, Hadley, and Wijekoon, 2011). LIDAR errors are also typically larger in areas of steep terrain, as any horizontal positional errors can result in large vertical errors (Goulden et al., 2016; Spaete et al., 2011; Su and Bork, 2006). Therefore, the estimated vertical uncertainty in DEMs constrained by LIDAR measurements should also vary spatially.
Vertical Datum Transformation Uncertainty
The development of integrated bathymetric–topographic DEMs typically requires the transformation of bathymetric and topographic measurements to a common vertical reference. Bathymetric data are usually referenced vertically to a tidal datum, such as mean lower low water (MLLW; Gill and Schultz, 2001), whereas topographic data are usually referenced to an orthometric datum, such as the North American Vertical Datum of 1988 (NAVD 88). Tidal datums are based on local observations of tidal variations over a period of time (Parker et al., 2003), whereas orthometric datums are based on Earth's gravity field (i.e. the geoid). For example, MLLW is defined as the arithmetic mean of the lower low water heights of the tide observed at a specific location over a 19-year period known as the National Tidal Datum Epoch (Gill and Schultz, 2001). A false vertical offset can result at the coastline where bathymetric and topographic data sets converge; if the data sets are not transformed to a common vertical reference, however, the transformation of bathymetric and topographic measurements to a common vertical datum adds additional vertical uncertainty into the DEM (Cooper and Chen, 2013; Gesch, 2013). The uncertainty originates from a combination of inaccuracies in the gridded fields used in the transformation, including the geoid, and in the source observation data used in the vertical datum transformation software, such as the elevation of the tidal datums or the height of the orthometric datum (NOAA, 2016). The incorporation of DEM uncertainty that originates from the vertical datum transformation is considered in multiple studies that evaluate the uncertainty of future flood risk due to sea-level rise (Gesch, 2013; Mitsova, Esnard, and Li, 2012; Schmid, Hadley, and Waters, 2013).
DEM Spatial Resolution and Cell-Level Measurement Uncertainty
Previous studies (e.g., Gao, 1997; Li, 1994; Shi, Wang, and Tian, 2014; Wechsler and Kroll, 2006) investigated the effect of DEM spatial resolution (i.e. cell size) on DEM vertical errors, and found that vertical errors generally increase at coarser resolutions (i.e. larger cell sizes). These previous studies quantified the errors by comparing the DEM elevation values to higher-accuracy, discrete ground control points, or a higher-resolution DEM. The ability of a DEM to accurately represent a terrain depends on the match between the DEM resolution and the spatial characteristics of the terrain (Fisher and Tate, 2006; Theobold, 1989). For example, a coarse-resolution DEM can more accurately represent a gently sloping beach than a steep beach cliff. Most current studies focus on the effect of DEM spatial resolution on the magnitude of vertical errors, and don't investigate the effect of DEM spatial resolution on the number of measurements per DEM grid cell and related information on subcell terrain variance. Hell and Jakobsson (2011) implement a multiple spatial resolution gridding approach to reduce interpolation artifacts in areas of sparse measurements, and subsequently improve the quality of the DEM, but do not produce uncertainty estimates using the number of measurements per DEM grid cell. The number of measurements per DEM grid cell, in conjunction with the subcell terrain variance, can estimate the DEM cell-level measurement uncertainty (Gleason, 2012; Wechsler, 2007).
Wechsler (2007) indicates that LIDAR technology typically provides multiple measurements per DEM grid cell, and that information on subcell terrain variance provided by these multiple measurements can be a useful component of DEM uncertainty estimations. Current state-of-the-art linear-mode LIDAR sensors have a data density of approximately 2–4 elevation returns/m2, and emerging single-photon LIDAR and Geiger-mode LIDAR have a data density of approximately 23 and 25 elevation returns/m2 for open terrain, respectively (Stoker et al., 2016). NOAA NCEI coastal DEMs have spatial resolutions that usually range from approximately 3 to 10 m (Amante and Eakins, 2016; Eakins and Grothe, 2014), resulting in multiple measurements per DEM grid cell where there is LIDAR coverage. DEM values typically represent a distance-weighted mean of all measurements located within an individual DEM grid cell when using an exact interpolation technique (Amante and Eakins, 2016; Caress and Chayes, 1996). The cell-level measurement uncertainty can therefore be expressed by the standard deviation of the mean, which is also commonly known as the standard error of the mean, or simply the standard error. The cell-level standard error depends on the measurement uncertainty and any vertical datum transformation uncertainty described in previous sections, the subcell terrain variance, and the number of measurements in a DEM grid cell at the defined spatial resolution.
Interpolation Uncertainty
A coastal DEM requires interpolation to estimate elevations in DEM grid cells not constrained by measurements to create a continuous surface and prevent instabilities while modeling coastal processes (Amante and Eakins, 2016). Interpolation techniques can be classified into general groups on the basis of the mathematical assumptions and features that estimate elevations for unmeasured locations using surrounding known measurements (Amante and Eakins, 2016). An important distinction between interpolation techniques in the context of DEM uncertainty is geostatistical vs. deterministic interpolation techniques.
Geostatistical Interpolation Techniques
Geostatistical techniques, such as kriging, are often utilized to generate surfaces from discrete measurements as they provide minimum-variance, linear unbiased estimations (Armstrong, 1998; Cressie, 1990; Matheron, 1963; Meyer, 2004). Kriging utilizes the semivariogram to estimate unknown elevations and to also predict their uncertainty (i.e. variance). A semivariogram captures the spatial correlation of the terrain by plotting the elevation variance of each pair of measurements as a function of the distance between the measurements, and then a mathematical model (e.g., linear, spherical, exponential) is fit to the semivariogram. There are numerous types of kriging, such as ordinary kriging, cokriging, and simple kriging, and each type has different statistical assumptions and constraints (Meul and Van Meirvenne, 2003). A main limitation of implementing geostatistical methods such as kriging is the large computational expense needed to create the semivariogram, especially with voluminous elevation data sets (Hell and Jakobsson, 2011). Geostatistical methods typically have computational costs that scale with the cube of the number of measurements (Cressie and Johannesson, 2008; Kleiber and Nychka, 2015). Attempts to optimize kriging methods, such as fixed-rank kriging (Cressie and Johannesson, 2008; Katzfuss and Cressie, 2011), may still not be feasible for developing coastal DEMs with tens to hundreds of millions of elevation measurements, which is common with dense LIDAR and multibeam sonar data sets (Katzfuss and Cressie, 2011).
Kriging is ideal when the terrain can be modelled as a stationary process with a constant variance (Detweiler and Ferris, 2010). Areas of coastal DEMs typically have different terrain morphologies, from relatively flat coastal plains to dynamic, coastal inlets with large terrain slope. Thus, the terrain is not a stationary process, and one model will not accurately capture the spatial structure of the entire DEM (Maune et al., 2007). Computational limitations and varying morphologies can necessitate dividing the region of interest into smaller sections; however, this approach can cause abrupt vertical offsets along the borders of the sections in the final, composite DEM (Memarsadeghi and Mount, 2007; Meyer, 2004). Another limitation of using kriging to develop coastal DEMs with accompanying uncertainty surfaces is the treatment of measurement uncertainty when integrating multiple data sets of disparate quality and age. The nugget of the semivariogram used in kriging represents microscale elevation variability and measurement uncertainty, and is calculated from the variance between elevation measurements at infinitesimally small distances apart (Clark, 2010; Cressie, 1993; Lythe and Vaughan, 2001). The measurement uncertainty typically varies throughout a coastal DEM when integrating multiple data sets of disparate quality and age. Thus, a global indicator of the measurement uncertainty, as provided by the semivariogram nugget, is of limited value for deriving coastal DEM uncertainty surfaces.
Zhang et al. (2015) recognize the limitation of ordinary kriging in that it does not consider measurement accuracy differences in areas where multiple data sources overlap. It is common practice to weight the contribution of data sets to the DEM value differently depending on their quality and age, with more recent, higher-quality data sets receiving a larger weight (Amante et al., 2011). Such data set weighting schemes are not easily implemented with geostatistical techniques. Cokriging could incorporate multiple data sets of disparate quality, and weight data sets on the basis of their semivariograms, but this version of kriging is even more computationally intensive as a multiple semivariograms require derivation (one for each data set). Furthermore, a weighting scheme based on the respective semivariograms would not incorporate the age of the data set. For example, a data set collected before a coastal storm could be more accurate (i.e. have a smaller variance), and, therefore, receive a larger weight than a hypothetically less-accurate (i.e. larger variance) poststorm data set. A DEM that aims to represent the present-day conditions should have the poststorm data set receive a larger weight, but this is not easily implemented with geostatistical techniques.
Deterministic Interpolation Techniques
Deterministic interpolation techniques, such as inverse distance weighting (IDW), triangulation, and spline, predict DEM values unconstrained by measurements, but notably, provide no estimates of their vertical uncertainty. Amante and Eakins (2016) and other previous research (e.g., Aguilar et al., 2005; Carlisle, 2005; Erdogan, 2009, 2010; Guo et al., 2010) indicate that terrain slope and curvature affect the accuracy of interpolation techniques. In general, interpolation techniques are less accurate in areas of large terrain slope and curvature. Amante and Eakins (2016) use Spearman's rank correlation to determine that IDW and triangulation deviations from measured depths are most positively correlated with terrain slope, whereas spline deviations are most positively correlated with terrain curvature. IDW and triangulation deviations are most positively correlated with slope because they are linear-weighted algorithms, and, therefore, local minima and maxima are not represented unless they are directly sampled. Spline deviations are most positively correlated with terrain curvature because its minimum curvature algorithm produces “overshoots” near areas of large curvature (Amante and Eakins, 2016).
Given the relationship between terrain slope and curvature and interpolation accuracy, the terrain can hypothetically predict interpolation uncertainty if there are dense data and interpolation is at short distances from measurements (i.e. a few DEM cells). Conversely, sparse depth measurements in coastal waters require interpolation at large distances (i.e. hundreds of cells) for the coastal DEM to retain the high spatial resolution of the topographic elevation measurements. In areas of large data gaps, the terrain is unknown and, therefore, it cannot be used directly to predict interpolation uncertainty. Instead of terrain, Amante and Eakins (2016) use the distance to the nearest measurement to derive predictive bathymetric interpolation uncertainty equations for Kachemak Bay, Alaska. These equations are limited because they do not incorporate measurement uncertainty from multiple, diverse bathymetric and topographic data sets, and because they are derived specifically for the terrain of Kachemak Bay (Amante and Eakins, 2016).
The relative accuracy of various interpolation techniques that generate DEMs, including deterministic techniques such as spline and geostatistical techniques such as kriging, varies depending on the terrain, data quality, and data density (Chaplot et al., 2006). NOAA NCEI develops coastal DEMs using spline interpolation for several reasons. Amante and Eakins (2016) determined that the accuracy of three deterministic methods (spline, IDW, triangulation) are approximately equivalent at short interpolation distances (one to two cells), but that spline is more accurate at large distances, and, therefore, is more appropriate for creating coastal DEMs with sparse bathymetric measurements. Spline interpolation also produces a smooth, gradually changing surface, which is representative of many coastal areas in the United States with gently varying terrain (Maune et al., 2007). LIDAR coverage in portions of this manuscript's study area in Florida (details forthcoming) indicates gently varying terrain, and, therefore, spline interpolation is an appropriate interpolation technique for the study area. The gradually changing elevation surface created by spline interpolation is also an important feature for coastal inundation models, as abrupt discontinuities can cause artificial barriers to water flow (Maune et al., 2007). The limitations of implementing kriging with multiple, diverse data sets previously described further justify utilizing spline interpolation instead of kriging to develop NOAA NCEI DEMs. A fundamental limitation with deterministic interpolation techniques, such as spline, is the lack of accompanying uncertainty estimates of the interpolated elevations. Methods to address this limitation of spline interpolation, in conjunction with methods to estimate the other sources of DEM uncertainty previously described, are the primary focus of this study.
Study Area and Coastal DEM Specifications
NOAA NCEI is developing coastal DEMs along the SW coast of Florida in a suite of 0.25° × 0.25° tiles. A DEM created for one 0.25° tile (bounding box: 26.75–27.00 N, 82.50–82.25 W) south of Sarasota, Florida is the case study to highlight the methods for deriving uncertainty surfaces (Figure 1). The DEM has a spatial resolution of one-ninth arc-second (∼3-m) and is referenced horizontally to the World Geodetic System (WGS 84) and vertically to NAVD 88. The study area has lowland elevations and shallow offshore depths that include the census-designated places of Englewood and Rotonda West, and the water bodies of the Gulf of Mexico, Lemon Bay, and the northern portion of Gasparilla Sound. The study area consists of mixed land use and land cover including residential development, marine, transportation, freshwater forested wetlands, mixed hardwood–coniferous, and improved pasture (Florida Fish and Wildlife Conservation Commission and Florida Natural Areas Inventory, 2016). The methods and source code are developed to create accompanying uncertainty surfaces for future NOAA NCEI DEMs in other locations in an automated manner.
METHODS
Methods to estimate potential DEM vertical errors at the individual cell-level that originate from the (1) elevation measurements, (2) vertical datum transformation, (3) DEM spatial resolution, and (4) interpolation technique are the primary focus of this study. Figure 2 provides an overview of the general methodology.
Figure 2.
General overview of the methodology for estimating the total vertical uncertainty at the individual DEM cell-level.
Software
MB-System (version 5.4.2220; Caress and Chayes, 1996) is the main software that generates the one-ninth arc-second coastal DEM. MB-System is a National Science Foundation-funded open-source software application specifically designed to manipulate multibeam sonar data, though it can utilize a wide variety of data types, including generic xyz data. The MB-System tool “mbgrid” applies spline interpolation to the xyz data to generate the coastal DEM (Amante et al., 2011). Several other open-source software programs including Generic Mapping Tools (GMT; version 4.5.13; Wessel et al., 2013), Geospatial Data Abstraction Library (version 2.1.0), Python computer language (version 2.7), as well as Unix utilities, including Grep, Awk, and Sed, in a Bash environment, aid in the derivation of the coastal DEM uncertainty surface.
Data Sources and Measurement Uncertainty
Several bathymetric and topographic point data sources are integrated to generate the coastal DEM (Table 1, Figure 3). Quality assessment and quality control are performed on each data set to identify and correct or remove obvious anomalies. Furthermore, newer data sets supersede older data sets where there is spatial overlap. Older data sets are spatially “masked” to newer data sets so that the DEM represents the most-recent elevations, and, therefore, the best approximation of the present-day terrain. The topographic and bathymetric data sets are collected with a data buffer 10% larger than the 0.25° × 0.25° extents to eliminate any potential interpolation edge effects (Amante and Eakins, 2016). Table 1 indicates the year and the vertical uncertainty of each data set at 1 standard deviation. The uncertainty is obtained from the published metadata for modern topographic and bathymetric–topographic LIDAR data sets, which are typically derived from ground control points or assumed technology standards. The data hierarchy used in the mbgrid gridding algorithm, as relative gridding weights, is also listed in Table 1. The weights are assigned on the basis of the overall quality assessment and the age of the data sets. Higher-quality and more recent data sets receive larger weights, and have greater influence on the predicted DEM value (Schmidt, Chayes, and Caress, 2006).
Table 1.
Data sets used to generate the DEM, their vertical uncertainty at 1 standard deviation, and their relative gridding weight in MB-System's mbgrid tool; larger weights have more influence on the predicted DEM value. The vertical uncertainty of the data sets vary, particularly between modern, LIDAR-derived data sets and legacy bathymetric data sets that date back to 1951.
Figure 3.
Spatial extent of the DEM source data. The data density varies throughout the study area, particularly with dense FDEM 2007 topographic LIDAR and sparse NOS Hydrographic Surveys. Note: The spatial footprints of the sparse bathymetry point data are enlarged by a factor of ∼9× to be visible at the scale of the map. (Color for this figure is available in the online version of this paper.)
Table 2 characterizes the uncertainty for bathymetric data sets derived from the NOAA ZOCs. The NOAA ZOC classification is generalized for this research solely on the basis of the age of the survey, and consequently, the presumed technology and quality standards, to estimate the uncertainty of bathymetric data sets. The National Oceanic Service (NOS) Hydrographic Surveys, the National Geophysical Data Center (NGDC) Multibeam, and U.S. Army Corps of Engineers (USACE) Dredge Surveys are all assigned uncertainty values according to the appropriate ZOC (Table 2). The NOS Hydrographic Surveys are collected for NOAA charting purposes, and, consequently, have rigorous quality standards. The NOS Hydrographic Surveys in the study area utilized single-beam sonar technology, and the vertical uncertainty is assigned a ZOC of B. The NGDC Multibeam data set is a collection of sonar-derived depths, typically from academic fleets. These data are not held to the same accuracy standards as NOS Hydrographic Surveys. Therefore, even though depths were collected with multibeam sonar, which would typically be assigned a ZOC of A for NOS Hydrographic Surveys, they are assigned a larger uncertainty of ZOC of B. The USACE Dredge Surveys utilized single-beam sonar, and are also assigned a ZOC of B. The ZOCs indicate potential vertical errors at the 95% confidence level. The errors are assumed to be normally distributed and divided by 1.96 to be consistent with the topographic and bathymetric–topographic LIDAR data set uncertainty provided in the data sets' metadata at 1 standard deviation.
Table 2.
Adapted zones of confidence (ZOC) to estimate the vertical uncertainty of bathymetric data sets. The vertical uncertainty is estimated on the basis of the age (and presumed technology and quality standards) of the data set. The variable d in the vertical uncertainty equation represents the water depth. Potential vertical errors are assumed to be normally distributed and are converted to 1 standard deviation by dividing by 1.96. The bathymetric data sets (USACE Dredge Surveys, NGDC Multibeam, and NOS Hydrographic Surveys) in this study area are all assigned a ZOC of B.
Vertical Datum Transformation Uncertainty
Two bathymetric data sets, the NOS Hydrographic Surveys and the USACE Dredge Surveys, are originally referenced vertically to the tidal datum of MLLW. The NOAA vertical datum transformation tool (VDatum) converts elevation data sets to common datums by considering the spatial variability of the relationship between the datums (Parker et al. 2003). VDatum transforms the NOS Hydrographic Surveys and the USACE Dredge Surveys from the tidal datum of MLLW to the orthometric datum of NAVD 88 (Geoid 12b) to be consistent with the topographic and bathymetric–topographic LIDAR data sets. VDatum provides a single, global estimate of vertical datum transformation uncertainty of 0.12 m at 1 standard deviation in this area of SW Florida (NOAA, 2016). The other bathymetric data set, the NGDC Multibeam, was not referenced to a specific tidal datum during collection. The depths are assumed to be referenced to the instantaneous water level, and no datum transformation is performed. The lack of specific tidal datum for the NGDC Multibeam also provides justification for the larger uncertainty designation (i.e. ZOC of B instead of ZOC of A). The measurement uncertainty listed in Tables 1 and 2 (σm) and any vertical datum transformation uncertainty (i.e. 0.12 m) provided by VDatum (σd) are considered independent, and are combined using the root sum of squares (Schmid, Hadley, and Waters, 2013) to calculate the data set source vertical uncertainty (SVU) at 1 standard deviation (Equation 1):
where, SVU = the source vertical uncertainty, σm = the measurement uncertainty, and σd = vertical datum transformation uncertainty.
DEM Cell-Level Source Uncertainty
The DEM cell-level source uncertainty is equivalent to Equation (1) where there is only one measurement in a DEM grid cell. Where there are multiple measurements in a DEM grid cell, the DEM cell-level uncertainty is calculated from the SVU from Equation (1), the subcell terrain variance, and the number of measurements. First, the exact pooled variance is calculated (Rudmin, 2010). “The exact pooled variance is the mean of the variances plus the variance of the means of the component data sets” (Rudmin, 2010, p. 1). The cell-level exact pooled variance (S2; Equation [2]) is equal to the square of the weighted mean of the SVU from Equation (1) for all measurements plus the weighted variance of all measurements around the weighted mean elevation, multiplied by the Bessel small-sample correction factor (Upton and Cook, 2014). The Bessel small-sample correction factor corrects the bias in the estimation of the cell-level variance, especially when the number of measurements in a DEM grid cell is less than 30. The relative data set weighting hierarchy provided in Table 1 is incorporated if there are elevation measurements in a DEM grid cell from more than one data set. Note that Equation (2) also calculates the exact pooled variance where there are multiple measurements in a DEM grid cell from the same data set, as the data set weight is equivalent for all measurements and therefore does not affect the calculation. The cell-level standard error is then calculated from the exact pooled variance in Equation (2), and the number of measurements per grid cell as determined by the spatial resolution of the DEM, to represent the DEM cell-level source uncertainty (Sz̄ ; Equation [3]):
where, S2 = the cell-level exact pooled variance, n = the number of elevation measurements in a given DEM grid cell, SVUi = the measurement and any vertical datum transformation uncertainty calculated from Equation (1) for the ith measurement, wi = the ith measurement data set weight provided in Table 1, zi = the ith measurement elevation value, and z̄ = the weighted mean elevation of all measurements in a given DEM grid cell. where, Sz̄ = the cell-level standard error, S2 = the exact pooled variance calculated from Equation (2), and n = the number of elevation measurements in a given DEM grid cell.
A source uncertainty surface at 1 standard deviation is then derived from the DEM cell-level standard error to reflect that the source uncertainty also propagates into interpolation uncertainty in cells unconstrained by measurements. The interpolation uncertainty (details forthcoming) is calculated by assuming the measured values are the “true” values; however, nearby source uncertainty contributes additional uncertainty in interpolated regions of the DEM. The DEM cell-level standard error calculations are associated with the latitude and longitude of the center of each cell constrained by at least one measurement. The GMT software “surface” tool creates the source uncertainty surface at a spatial resolution of one-ninth arc-second with the cell-level standard error point data using spline interpolation with an adjustable tension value. A tension value of 0.35 is utilized to suppress undesired oscillations and false local maxima or minima (Smith and Wessel, 1990). Furthermore, a lower-limit value of zero is imposed on the output source uncertainty surface to prevent negative uncertainty values in any areas of false local minima created by the spline interpolation.
Interpolation Uncertainty
A split-sample method quantifies interpolation deviations from measured values to derive an interpolation uncertainty equation (Amante and Eakins, 2016). The split-sample method is applied to smaller subgrids within the study area to quantify interpolation deviations for terrains with different slopes and curvatures. The interpolation deviations from the smaller subgrids are associated with the distance to the nearest measurement, and then aggregated to derive a single interpolation uncertainty equation to apply to the entire study area.
Split-Sample Method
A split-sample method consists of randomly omitting a percentage of measurements, applying an interpolation technique, and calculating the differences between the interpolated values and the omitted measurements (Amante and Eakins, 2016). During each split-sample routine, the retained measurements are gridded using spline interpolation with the MB-System tool mbgrid, and the resulting interpolated raster is compared, on a cell-by-cell basis, with the omitted measurement raster to quantify the interpolation deviations. For each cell, the Euclidean distance to the nearest measurement is calculated, measured in raster cell units. Each interpolation deviation is then associated with the distance to the nearest measurement. See Amante and Eakins (2016) for more details on the split-sample method.
DEM Subgrids
The terrain slope and curvature affect the magnitude of spline interpolation deviations from measured values (Amante and Eakins, 2016). Therefore, the DEM is divided into smaller subgrids at the same spatial resolution (i.e. one-ninth arc-second) to perform the split-sample method on different terrains throughout the study area. The number of rows and columns of the subgrids is determined by calculating the distance to the nearest measurement for every cell in the study area. The maximum value of the distance to the nearest measurement, i.e. the maximum interpolation distance, is ∼136 cells and the 95th percentile is ∼59 cells. Subgrid dimensions are automatically generated as four times the 95th percentile of the distance to the nearest measurement for the entire study area, which equates to 236 rows by 236 columns (∼0.5 km2). These dimensions ensure that a statistically significant number of interpolation deviations are quantified using the split-sample procedure for almost all interpolation distances in the DEM.
Split-Sample Subgrids: Criteria and Selection
A stratified, semirandom sampling approach selects the subgrids for split-sample routines, ensuring that the split-sample subgrids are located in areas of relatively dense data and are geographically located throughout the DEM (Figure 4). The cell-sampling density is defined as the proportion of DEM grid cells constrained by measurements, and determines the eligibility of subgrids for split-sample routines. Areas of large cell-sampling densities are preferred as they produce more interpolation deviations, i.e. more samples, to derive the relationship between interpolation deviations from measured values and the distance to the nearest measurement. Furthermore, areas of large cell-sampling densities are typically constrained by higher-accuracy measurement technologies, such as LIDAR. Thus, the measured elevations in these areas depict realistic terrain, which is ideal for quantifying interpolation deviations to subsequently derive the interpolation uncertainty equation.
Figure 4.
Location of the split-sample subgrids for the bathymetry (Bathy), bathymetry–topography (Bathytopo), and topography (Topo) strata to quantify interpolation deviations and derive the interpolation uncertainty equation. The subgrids are in areas of relatively dense data (areas of black), especially for the BathyTopo and Topo strata, and are also geographically located throughout the DEM area to incorporate the effect of terrain slope and curvature on the magnitude of interpolation deviations. Note: The spatial footprints of the sparse bathymetry point data are enlarged by a factor of ∼9× to be visible at the scale of the map. (Color for this figure is available in the online version of this paper.)
The subgrids are initially assigned to three strata on the basis of their elevation values: bathymetry (bathy), bathymetry–topography (bathytopo), and topography (topo). Subgrids with all DEM values below the NAVD 88 zero elevation are bathy, subgrids with all DEM values above zero are topo, and subgrids with DEM values below and above zero are bathytopo. A cell-sampling density percentile threshold is specified for each stratum to ensure that areas of dense data are selected. In this study, all subgrids equal to or greater than the 50th percentile of the cell-sampling density for their respective stratum are eligible for split-sample selection. Twenty-five subgrids for each stratum are then selected for split-sample routines by maximizing the cumulative distance between all subgrids in that stratum, for a total of 75 subgrids (Figure 4).
Derived Interpolation Uncertainty Equation
The split-sample method is implemented on all 75 selected subgrids shown in Figure 4. The split-sample percentage determines the number of measurements retained for interpolation (i.e. training data). The split-sample percentage is automatically determined as the 5th percentile of the cell-sampling density from all DEM subgrids in the study area. The 5th percentile of the cell-sampling density is 0.007%, equating to retaining elevation values for four of the 55,696 raster cells in a subgrid for training in each split-sample routine. Elevations along the outermost edge of the subgrids are also retained to guide interpolation to avoid interpolation edge effects (Amante and Eakins, 2016). Amante and Eakins (2016) indicate that the magnitude of the interpolation deviations decrease at the same distance to the nearest measurement when increasing the cell-sampling density. Thus, deriving an interpolation uncertainty equation using the lower limit of the cell-sampling density (the 5th percentile) ensures a liberal uncertainty estimate, as the cell-sampling density will be larger for most areas of the DEM.
The split-sample method is performed 50 times for each of the 75 subgrids for an aggregated total of 3750 split-sample routines. Fifty million interpolation deviations from the original measurements and their associated distance to the nearest measurement are randomly selected from the aggregated split-sample routines. The interpolation deviations with distances up to the 95th percentile of the distance to the nearest measurement for the entire DEM (i.e. ∼59 cells) are separated into 10 equal-width bins of 5.9 cells. The standard deviation is calculated for each bin to derive an equation representing the interpolation uncertainty as a function of distance to the nearest measurement with a best-fit power law equation (Equation [4]). The interpolation uncertainty equation is then applied to a one-ninth arc-second raster representing the distance to the nearest measurement to derive the interpolation uncertainty surface at 1 standard deviation:
where, I(d) = the interpolation uncertainty at 1 standard deviation, d = the distance to the nearest measurement in raster cells, and A and B are derived coefficients.
Total Vertical Uncertainty
The DEM source uncertainty surface and interpolation uncertainty surface are assumed to be independent, and the total vertical uncertainty (TVU) surface at 1 standard deviation is calculated as the root sum of squares (Equation [5]):
where, TVUsurface = the total vertical uncertainty surface, σs_surface = the source uncertainty surface, and σi_surface = interpolation uncertainty surface.
RESULTS
The primary result of this study is an uncertainty surface that estimates potential DEM vertical errors at the individual cell-level that originate from the (1) elevation measurements, (2) vertical datum transformation, (3) DEM spatial resolution, and (4) interpolation technique.
Source Uncertainty
The DEM cell-level source uncertainty is calculated for every cell constrained by at least one measurement, and is used to derive the source uncertainty surface. The measurement uncertainty, any vertical datum transformation uncertainty, subcell terrain variance, and the number of measurements per grid cell determine the DEM cell-level source uncertainty using Equations (1), (2), and (3). The number of measurements per DEM grid cell is determined by the spatial resolution of the DEM (i.e. one-ninth arc-seconds), and is shown in Figure 5. The number of measurements per grid cell indicates areas of relatively dense and sparse data. In general, areas with LIDAR and multibeam sonar data coverage have many measurements per DEM grid cell, and consequently, small vertical uncertainty. Areas of sparse data include deeper waters that are constrained by a single bathymetric measurement per DEM grid cell, and have large vertical uncertainty. An example of the cell-level source uncertainty for an individual DEM grid cell and its location within the DEM source uncertainty surface is shown in Figure 6. There is small uncertainty (±0.015 m at 1 standard deviation) in the individual DEM grid cell highlighted in Figure 6 due to small source data set uncertainty from topographic LIDAR (Florida Division of Emergency Management [FDEM] 2007 LIDAR), small subcell terrain variance in flat terrain, and 40 measurements in the DEM grid cell due to effective laser pulse penetration in nonvegetated land cover.
Figure 5.
Number of measurements per grid cell affects the magnitude of the cell-level source uncertainty; many measurements per grid cell result in small uncertainty of the DEM value. Areas of dense data are typically found on land where there is topographic LIDAR and offshore within a single swath of multibeam sonar. Areas of dark vertical stripes on land indicate where LIDAR flight paths overlap and effectively double the number of measurements per grid cell, resulting in smaller uncertainty. Areas of sparse data, such as deeper waters, are constrained by a single bathymetric measurement per DEM grid cell and have large uncertainty. Note: The spatial footprints of the sparse bathymetry point data are enlarged by a factor of ∼9× to be visible at the scale of the map.
Figure 6.
An example of the cell-level source uncertainty (standard error) for an individual DEM grid cell (left diagram) and its location within the DEM source uncertainty surface (right diagram). The cell-level source uncertainty for an individual DEM grid cell (left diagram) is calculated from the measurement uncertainty (error bars), the number of measurements (circles), and subcell terrain variance (vertical position of measurements around the dashed line representing the mean elevation). Plus/minus the pooled standard deviation from the mean elevation is shown in light gray, and plus/minus the standard error at 1 standard deviation from the mean elevation is shown in dark gray. The source uncertainty surface (right diagram) is derived from the cell-level source uncertainty (left diagram) and varies spatially, depending on the measurement uncertainty, any vertical datum transformation uncertainty, the number of measurements per DEM grid cell, and subcell terrain variance around the mean elevation of the DEM grid cell. The location of the individual DEM grid cell in the left diagram is indicated by the star in the right diagram, and its relationship to land cover and terrain slope is indicated by the star in Figure 7.
Figure 7.
Effect of land cover and terrain slope on the cell-level source uncertainty (i.e. standard error). Areas of dense vegetation (veg.) have large uncertainty due to few LIDAR ground returns per DEM grid cell due to poor laser pulse penetration. Areas of large terrain slope have large uncertainty due to large subcell terrain variance around the mean elevation. Conversely, flat terrains with sparse vegetation have small uncertainty due to small subcell terrain variance and many LIDAR ground returns per DEM grid cell due to effective laser pulse penetration. The area of low uncertainty indicated by the star is the location of the individual DEM grid cell shown in Figure 6.
The effect of land cover and terrain slope on the magnitude of the cell-level source uncertainty is further illustrated in Figure 7, and statistically quantified in Table 3. Areas of denser vegetation and larger slopes have larger cell-level source uncertainty. Four land cover classes in the study area are identified from the Florida Cooperative Land Cover Map (Florida Fish and Wildlife Conservation Commission and Florida Natural Areas Inventory, 2016) to highlight the effect of land cover on the number of measurements per DEM cell, and consequently, the magnitude of the cell-level source uncertainty (i.e. standard error). Land cover classes representing dense vegetation, “Freshwater Forested Wetlands” and “Mixed Hardwood–Coniferous,” have a relatively small average number of measurements per DEM cell, ∼9 and ∼10 measurements, due to poor laser pulse penetration, and relatively large average standard error, 0.063 m and 0.061 m, respectively (Table 3). Classes representing sparse vegetation, “Improved Pasture” and “Transportation,” have a relatively large average number of measurements per DEM cell, ∼22 and ∼33 measurements, due to effective laser pulse penetration, and relatively small average standard error, 0.030, and 0.028 m, respectively (Table 3). The Spearman's rank correlation coefficient (Spearman, 1904; Table 3) indicates a negative correlation between the standard error and the number of measurements for all four land cover classes, i.e. more measurements per DEM cell results in smaller cell-level standard error, as expected per Equation (3). Furthermore, terrain slope is positively correlated with standard error for all classes, as larger terrain slope results in larger subcell terrain variance, and, consequently, larger standard error, as expected per Equations (2) and (3). The average terrain slope for each of the four classes is less than 2° and the relatively flat terrain results in a smaller-magnitude Spearman's rank correlation coefficient between terrain slope and standard error than the correlation coefficient between the number of measurements and standard error for each land cover class (Table 3). The number of DEM cells, i.e. sample size, for each of the four land cover classes is greater than 250,000, and all Spearman's rank correlation coefficients have p-values < 0.001.
Table 3.
Effect of land cover and terrain slope on the magnitude of the cell-level source uncertainty (standard error). Land cover classes representing dense vegetation, “Freshwater Forested Wetlands” and “Mixed Hardwood–Coniferous,” have a relatively small average number of measurements per DEM cell due to poor laser pulse penetration, and consequently, have relatively large average standard error. Classes representing sparse vegetation, “Improved Pasture” and “Transportation,” have a relatively large average number of measurements per DEM cell due to effective laser pulse penetration, and consequently, have small average standard error. Spearman's rank correlation coefficient indicates a negative correlation between the standard error and the number of measurements for all four land cover classes, i.e. more measurements per DEM cell result in smaller cell-level standard error, as expected. Terrain slope is positively correlated with standard error for all classes, as larger terrain slope results in larger subcell terrain variance, and, consequently, larger standard error, as expected. The number of DEM cells, i.e. sample size, for each of the four land cover classes is greater than 250,000, and all Spearman's rank correlation coefficients have p-values < 0.001.
Interpolation Uncertainty
The interpolation deviations from measured values are plotted as a function of distance to the nearest measurement in panel A of Figure 8. The deviations are not biased, with a mean of approximately 0 m for all distances to the nearest measurement. The interpolation deviations are separated into 10 equal-width bins of 5.9 cells, up to the 95th percentile of the distance to the nearest measurement for the entire DEM (∼59 cells). The interpolation uncertainty equation is derived from the standard deviation of the binned interpolation deviations as a function of distance to the nearest measurement (panel B of Figure 8). The magnitude of the interpolation uncertainty in panel B increases with larger distances to the nearest measurement, and then levels off as the spatial autocorrelation of the terrain decreases. The interpolation uncertainty surface and its relationship to the distance to the nearest measurement is shown in Figure 9. There is large interpolation uncertainty in areas of sparse, bathymetric measurements offshore due to large distances to the nearest measurement. There is small interpolation uncertainty on land and along the coastline in areas of dense LIDAR measurements due to small distances to the nearest measurement.
Figure 8.
Fifty million interpolation deviations from measured values from the split-sample routines (panel A) and the derived interpolation uncertainty equation (panel B). The deviations in panel A are unbiased, with a mean of approximately 0 m. Also, note that the magnitude of the interpolation uncertainty in panel B increases with distance to the nearest measurement, and then levels off as the spatial autocorrelation of the terrain decreases.
Total Vertical Uncertainty
The total uncertainty surface estimates potential DEM vertical errors at the individual DEM cell-level at 1 standard deviation (Figure 10, panel C). The total uncertainty surface is calculated from the root sum of squares of the source uncertainty surface (panel A) and the interpolation uncertainty surface (panel B) in Figure 10. The largest estimations of potential vertical errors are due to a combination of large measurement uncertainty and additional vertical datum transformation uncertainty, few measurements per grid cell due to sparse data and the high spatial resolution of the DEM, and large interpolation uncertainty. Areas of large vertical uncertainty are in deeper waters where there are sparse, less-accurate NOS Hydrographic Survey measurements with transformed vertical datums, and large interpolation distances. Areas of small vertical uncertainty are located on land surfaces with small terrain slopes constrained by dense, relatively accurate FDEM 2007 topographic LIDAR measurements with no vertical datum transformation, and small interpolation distances.
Figure 10.
Source uncertainty surface (panel A), interpolation uncertainty surface (panel B), and the total uncertainty surface (panel C) at 1 standard deviation. The total uncertainty surface represents potential DEM deviations from the actual seabed or land surface. Note that the total uncertainty varies spatially, and reflects the spatial variability of both the source and interpolation uncertainty surfaces. Large vertical uncertainty exists in deeper waters offshore with sparse, echo-sounder measurements. Conversely, small vertical uncertainty exists on flat terrains with dense LIDAR measurements. The total uncertainty surface is the uncertainty product that should be incorporated in coastal process modeling.
DISCUSSION
The methods and results of this study advance previous research by Jakobsson, Calder, and Mayer (2002), Calder (2006), and Hell and Jakobsson (2011) by integrating bathymetric and topographic data sets of disparate age, quality, data density, and vertical datums to create a seamless coastal DEM, and most notably, provide an estimate of the cell-level uncertainty that also incorporates the spatial resolution of the DEM. Uncertainty estimations using the number of measurements per DEM grid cell, as determined by the spatial resolution of the DEM, are remarkably absent from literature on estimating DEM uncertainty, despite being acknowledged by Wechsler (2007) a decade ago.
The TVU surface (Figure 10, panel C) is the uncertainty product that should be incorporated in coastal process modeling. DEM realizations can be created by adding or subtracting the TVU surface at a desired confidence level from the DEM. For example, the TVU surface multiplied by ±1.96 for the 95% confidence level, and then added to the DEM separately, would result in maximum and minimum DEM realizations, respectively. Furthermore, intermediate realizations can be created between the maximum and minimum DEM realizations by multiplying the TVU surface by factors between −1.96 and 1.96, and then adding the resulting uncertainty surfaces to the DEM separately.
Total Vertical Uncertainty
The relative contribution of the measurement, vertical datum transformation, and interpolation uncertainty to the TVU varies throughout the study area and depends on the data set constraining the DEM. The measurement uncertainty is the largest contributor to the TVU for older bathymetric measurements in deeper waters, in accordance with the equation for ZOC B in Table 2. For example, at depths of 18 m, the measurement uncertainty at 1 standard deviation for the NOS Hydrographic Surveys data set is approximately 0.7 m. This measurement uncertainty dominates until distances of approximately 40 cells from the nearest measurement, where the interpolation uncertainty becomes a larger contributor. Conversely, with relatively accurate LIDAR technology, the interpolation uncertainty contribution is larger than the measurement uncertainty when the distance to the nearest measurement is larger than a few cells. The measurement uncertainty is generally much larger than the vertical datum transformation uncertainty (0.12 m), especially in deeper waters, for data sets that use VDatum to transform the data from the tidal datum of MLLW to the orthometric datum of NAVD 88. The vertical datum transformation is also smaller than the interpolation uncertainty for essentially all interpolation distances. Therefore, the vertical datum transformation uncertainty is not a primary contributor to DEM uncertainty in areas of sparse, old, relatively inaccurate bathymetric measurements.
Source Uncertainty
A limitation of the research in this manuscript is that the measurement uncertainty is represented by the global statistic provided in the data sets' metadata. Previous research indicates that the uncertainty of elevation measurements, particularly with LIDAR, is correlated with terrain and land cover, with larger uncertainty in areas of larger terrain slope and dense vegetation (Bater and Coops, 2009; Goulden et al., 2016; Leon, Heuvelink, and Phinn, 2014; Spaete et al., 2011; Su and Bork, 2006). Likewise, horizontal errors can produce large vertical errors in areas of large terrain slope with hydrographic data (Calder, 2006; Marks and Smith, 2008). Future research could improve the estimation of topographic measurement uncertainty by collecting accurate ground control points, correlating measurement error with terrain and land cover, and deriving spatially varying measurement uncertainty estimations (Leon, Heuvelink, and Phinn, 2014). The methodology presented in this manuscript is designed specifically for the development of NOAA NCEI DEMs, and the collection of ground control points for the numerous DEMs developed annually for coastal locations around the United States is not feasible because of limited resources. The current methodology does, however, partially incorporate terrain and land cover effects on the magnitude of the cell-level source uncertainty. The exact pooled variance in Equation (2) is larger in areas of steep slopes because of larger subcell terrain variance around the mean elevation. Furthermore, there are fewer LIDAR ground returns in densely vegetated areas, and thus the standard error calculation using Equation (3) will also be larger when dividing by a smaller number of measurements, n. This partial incorporation of terrain and land cover effects on the magnitude of the cell-level source uncertainty is illustrated in Figure 7 and statistically quantified in Table 3.
Another limitation related to the measurement uncertainty is the assumption of a normal error distribution. Previous research indicates that DEM errors can have a nonnormal error distribution (e.g., Marks and Smith, 2008; Schmid, Hadley, and Waters, 2013). Likewise, limited resources prevent the collection of additional data to determine the specific error distribution and necessitate the normality assumption. The normality assumption may result in an overestimation of the vertical uncertainty in DEMs (Schmid, Hadley, and Waters, 2013).
A final limitation of the source uncertainty estimation is the assumption that the datum transformation uncertainty is uniform across the study area. In areas of complex bathymetry, the relationship between orthometric and tidal datums can vary substantially, and, therefore, have spatially varying uncertainty (NOAA, 2016). The global metric of VDatum uncertainty in this research is, thus, another limitation. There is ongoing research at NOAA to create VDatum uncertainty surfaces, and when completed, this spatially varying vertical datum transformation uncertainty surface will easily be able to be incorporated into the methods described in this manuscript.
Spatial Resolution, Cell-Level Source Uncertainty, and Sample Size
For DEM grid cells constrained by measurements, the cell-level source uncertainty is represented by the standard error. In cells constrained by one measurement, the standard error is simply the data set SVU calculated from Equation (1) because the denominator, the number of measurements, n, in Equation (3) is equal to 1. This is the best approximation of the cell-level source uncertainty, but it is of extremely limited statistical value as there needs to be at least two measurements in a DEM grid cell for any useful metric regarding the uncertainty of an average value. Furthermore, even with a few measurements in a DEM grid cell, the small sample size requires a correction factor from the Student t distribution to convert the standard error, with an assumed normal distribution, to a desired population confidence level (e.g., 95% confidence; Student, 1908). Accordingly, another useful product to be disseminated with NOAA NCEI DEMs is an accompanying grid representing the number of measurements per DEM grid cell, similar to Figure 5. This grid can inform the appropriate Student t correction factor to express the uncertainty at a desired population confidence level. Future research will investigate additional techniques to estimate the uncertainty with small sample sizes, and to properly propagate the DEM uncertainty into the modeling of coastal processes.
One avenue of future research is to develop multiresolution, raster-based DEMs (Hell and Jakobsson, 2011; Wechsler, 2007). Vector-based DEMs, such as triangular irregular networks, allow for spatially varying resolutions on the basis of data density and terrain variance, but the unstructured nature is often not supported or computationally efficient in many modeling algorithms because of the complexity of computational geometry (de Azeredo Freitas et al., 2016; Shingare and Kale, 2013). Future research will investigate developing multiresolution raster grids that allow for different resolutions on the basis of data density, terrain variance, and vertical uncertainty. Hell and Jakobsson (2011) reduce DEM artifacts from spline interpolation introduced by large interpolation distances in areas of sparse measurements, while maintaining terrain details in areas of dense data, by generating a stack of multiple-resolution raster DEMs. Essentially, the higher-resolution grid cells overrule lower-resolution grid cells in the stack where there are sufficient data, and a composite, multiresolution DEM is generated (Hell and Jakobsson, 2011). However, Hell and Jakobsson (2011) do not provide estimates of the vertical uncertainty in the intermediate DEMs or in the final, composite multiresolution DEM.
Hell and Jakobsson (2011) use only the data density to locally adjust the DEM resolution. Future research will utilize the cell-level uncertainty estimates generated with the methods described in this manuscript, which incorporate the data density (i.e. number of measurements per DEM grid cell), in addition to the subcell terrain variance, and source measurement uncertainty, to iteratively adjust the local DEM resolution on the basis of a user-defined uncertainty limit. For example, a limit of 1-m TVU at 1 standard deviation can be established, and several DEMs are generated at progressively coarser resolutions until all DEM grid cells are below the 1-m uncertainty limit. Coarsening the resolution will result in less uncertainty in the average elevation within the footprint of a DEM grid cell for many coastal areas because of more measurements, especially in flat terrains where there is small subcell terrain variance. The resulting stack of DEMs with different spatial resolutions will be compared on a cell-by-cell basis depth-wise, with the highest resolution under the uncertainty limit being represented in the final, composite DEM, using similar methods as Hell and Jakobsson (2011).
This future research will also result in a more statistically robust standard error calculation in areas of sparse measurements, as there will be more measurements in DEM grid cells at coarser spatial resolutions. The specified uncertainty limit and resulting spatial resolution, however, will need to be balanced with the relevant scale of analysis of the DEM application (e.g., coastal inundation modeling, habitat modeling, visualization). For example, a 30-m2 DEM grid cell that contains hundreds of measurements may have a desired low vertical uncertainty of the average elevation within the cell footprint, but this coarse cell size would not be useful for detailed coastal inundation modeling, as significant volumes of water flow within conduits of smaller spatial dimensions.
Interpolation Uncertainty
The interpolation uncertainty equation in panel B of Figure 8, and the resulting interpolation uncertainty surface illustrated in panel B of Figure 9, are global estimates that are derived from numerous terrains throughout the study area to provide an intermediate approximation of interpolation uncertainty. Future research could derive equations to better incorporate the effect of local terrain slope and curvature on the magnitude of interpolation deviations. Interpolation uncertainty equations could be derived for each split-sample subgrid. These separate equations would then be applied to nearby subgrids using a distance-weighted algorithm to produce a continuous, but varying, estimate of interpolation uncertainty across the entire DEM to incorporate the effect of local terrain slope and curvature on the magnitude of interpolation uncertainty.
Morphologic Change
NOAA NCEI DEMs are developed to represent the most-recent data sets, and, therefore, the best approximation of the present-day terrain. Consequently, newer data sets supersede older data sets, and older data sets are removed before DEM generation. Future research will estimate additional vertical uncertainty due to potential morphologic change from the data collection date. For instance, dynamic areas, such as coastal inlets, have additional uncertainty due to morphologic change since the data collection date. An additional uncertainty contribution can be calculated from the cell-level measurement variance among data sets from multiple time periods. The cell-level measurement variance from multiple time periods can identify areas prone to morphologic change from natural sediment transport or storm events. The propensity for morphologic change can be another uncertainty component, especially in between data collections, or when using the DEM to model future coastal processes, such as sea-level rise inundation. For example, researchers modeling future coastal inundation from sea-level rise should incorporate uncertainty in the DEM due to potential, future morphologic change. Morphologic change analysis in dynamic coastal areas is now possible where there are multiple, accurate coastal LIDAR surveys spanning close to a decade, including this study area in Florida. Morphologic change analysis can also indicate areas that are relatively stable vs. dynamic, which aids in prioritizing future data collection areas.
CONCLUSIONS
Integrating several bathymetric and topographic data sets to create a coastal DEM and estimating its vertical uncertainty at the individual cell-level provide numerous challenges (Eakins and Grothe, 2014). The diverse data sets are typically collected with a wide range of technology, at different time periods, and referenced to different vertical datums. Consequently, the data sets have disparate measurement uncertainty, with additional uncertainty introduced by any vertical datum transformation. Furthermore, the incongruent data densities of bathymetric and topographic data sets often necessitate extreme interpolation between sparse bathymetric measurements for the DEM to retain the high spatial resolution of topographic LIDAR. This extreme interpolation adds additional uncertainty into the DEM, especially in areas of large terrain slope and curvature. The DEM spatial resolution, and subsequently the number of measurements per DEM grid cell and the subcell terrain variance, also affect the magnitude of the cell-level vertical uncertainty.
Methods to address the current lack of research on developing coastal DEMs with accompanying uncertainty surfaces that estimate potential DEM vertical errors at the individual cell-level are provided in this manuscript. DEM uncertainty affects the fidelity of coastal process modeling, such as tsunami propagation and coastal inundation. Incorporating the estimated uncertainty by deriving numerous plausible DEM realizations within the uncertainty bounds can produce an ensemble of plausible coastal process modeling results, and in turn, better-informed coastal management decisions. Estimating the spatially varying DEM uncertainty also aids in prioritizing future elevation data collection, which will subsequently also improve the fidelity of coastal process modeling.
ACKNOWLEDGMENTS
This research is part of my doctoral dissertation at the University of Colorado Boulder, and was supported by funding from NOAA NCEI in Stennis, Mississippi. I thank Waleed Abdalati, Barbara Buttenfield, Stefan Leyk, Seth Spielman, and Weiqing Han for their guidance throughout my doctoral studies. I also thank many at NOAA NCEI for their support, and Mike Sutherland and the two anonymous reviewers for their valuable feedback that significantly improved the quality of the manuscript.
