Fusing Digital Elevation Models to Improve Hydrological Interpretations ()
1. Introduction
There is a growing demand for high-resolution datasets representing the Earth’s surface as technology and applications of spatial modeling advance. Among these datasets are digital elevation models (DEMs), i.e., gridded datasets representing continuous elevation changes across landscapes at both local and global extents [1] [2] . In turn, DEMs are used for spatial visualization and modeling of topographic, geomorphologic, and hydrological properties (e.g., slopes, soil erosion, basin borders, stream and river networks, and stream discharge) [3] [4] [5] [6] . In this regard, DEM elevation accuracies and resolutions vary based on differences in mode of elevation capture, processing, point-to-point interpolations, and timing of capture [7] [8] [9] [10] [11] . Hence, slope and hydrological DEM interpretations such as DEM-delineated basin boundaries and flow channels can vary in considerable detail [12] [13] [14] .
In terms of acquisition and availability, DEMs of varying origins are becoming freely accessible [8] [10] [15] - [23] . Some of these are listed in Table 1 [19] [20] [21] [22] [23] , i.e., the Shuttle Radar Topography Mission (SRTM90 and SRTM30 DEMs) and the Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER DEM). These are satellite derived and global. Also listed are: the nationally available DEM for Canada, i.e., the Canadian Digital Elevation Model (CDED, generated from elevation contours and spot heights), and a provincial example, e.g., the New Brunswick DEM (NB-DEM) derived through stereo-photogrammetry. Also increasingly available are Light Detection and Ranging (LiDAR) DEMs generated from point-cloud data through air-borne laser-light scanning and pulse return classification [24] .
Elevation differences between LiDAR and other DEMs (hereby referred to as
![]()
Table 1. Overview of open-sourced DEMs and LiDAR-DEMs utilized in this study. Note variation in coverage, resolution, and vertical error.
1) http://www.snb.ca/gdam-igec/e/2900e_1c.asp , 2) https://earthexplorer.usgs.gov/ 3) http://geogratis.gc.ca/site/eng/extraction/ , 4) http://www.snb.ca/geonb1/e/DC/catalogue-E.asp
non-LiDAR DEMs) vary by methodology, accuracy and resolution [21] [24] [25] [26] [27] . Non-LiDAR DEMs such as SRTM30 and STRM90 are based on surface reflections, while ASTER, CDED, NB-DEM represent bare-earth elevations. In this regard, the former refer to Digital Surface Models (DSMs), while the latter refer to Digital Elevation Models (DEMs). Both are subsets of Digital Terrain Models (DTMs) [5] . LIDAR point cloud data can be used to generate DSMs as well as DEMs through determining elevation differences between the first and last laser pulse returns, generally at a vertical accuracy of ± 0.15 m [24] .
All DEMs can be used to generate hydro-topographic interpretations pertaining to, e.g., flow directions, flow accumulation, stream channel networks, upslope basin areas and borders, location of depressions, and extent of areas subject to flooding and well-to-poor soil drainage. These interpretations are, however, influenced by DEM accuracy [28] - [34] . For example, 10 to 20 m DEMs are preferable to determine seamless flow-channel connectivities across roads, but lead to underestimating slope steepness especially along shorelines, road cuts and deeply incised stream valleys. High-resolution LiDAR DEMs can serve both purposes, by using them at their finest resolution to determine where roads potentially block channel flow and therefore need to be breached, or where the DEMs need to be re-sampled towards coarser resolutions to approximately connect flow channels across barriers without breaching. In contrast, attempts to increase the resolution of photogrammetrically-derived non-LiDAR DEMs through re-interpolation alone do not generate more information, but introduce DEM artifacts such as “ridging” [35] .
Several approaches are available to reduce vertical and lateral DEM errors through fusion. These involve DEM re-projecting, re-sampling, re-interpolation, amalgamation, and hydrological enforcement (e.g., [36] [37] [38] [39] [40] ). For example, the extent of “ridging” can be reduced through DEM editing by way of Tin Random Densification [22] and/or spatial filtering [41] [42] . Luedeling et al. [43] applied a DEM fill technique to fill SRTM voids with ASTER data. Fusing DEMs for the purpose of DEM-error reduction was reviewed by [39] . Apart from DEM fusion, hydrological enforcement ( [44] [45] [46] [47] ) is needed to guide DEM-based flow-direction and flow-accumulation algorithms towards already delineated streams, rivers, lakes and shorelines. This enforcement is done by (i) lowering the elevation of all stream and open-water referenced pixels to their immediate lowest neighborhood elevations, and (ii) ensuring that all streams and rivers drop monotonously in elevation towards their nearest shorelines. This enforcement is generally applied to non-LiDAR DEMs, but hydro- topographical adjustments are also needed for high-resolution LiDAR-DEMs to remove artifacts due to low laser pulse returns from open water surfaces such as lakes, rivers and shores.
The objective of this article is to demonstrate how a systematic non-LiDAR- to-LiDAR reduction in elevation differences by way of non-LIDAR DEM fusion at 10 m resolution leads to DEM-generated hydrological interpretation improvements. This was done comprehensively for the Province of New Brunswick as a case study. The extent of the DEM-based interpretation improvements are analyzed below in terms of:
1) non-LIDAR DEM fusion;
2) using bare-earth 1m LiDAR elevation data as reference DEM;
3) nearest-distance conformance testing of non-LiDAR versus LiDAR-derived flow channels;
4) evaluating false positives and false negatives within DEM-derived cartographic depth-to-water indices (DTW [32] ) for each DEM.
2. Methodology
Located in Eastern Canada and spanning an area of 7,282,014 ha, the Province of New Brunswick encompasses an array of diverse geomorphologies and landforms, and has province-wide non-LiDAR DEM and hydrographic network coverages available from Service New Brunswick’s GeoNB data catalogue, and LiDAR coverages for parts of the province (Figure 1, Table 1). The SRTM and ASTER DEMs were acquired through the United States Geological Survey Branch’s Earth Explorer application. The CDED data were obtained from Natural Resources Canada. A semi-automated ArcGIS 10.1 five-stage process (Figure 2) was developed as follows:
![]()
Figure 1. LiDAR DEM coverages (~11%) used for improving non-LiDAR DEMs across all of New Brunswick by way of a 5-stage process. Also shown: latest LiDAR-DEM acquisition for New Brunswick, LiDAR coverage for optimization and process validation, and scan line (green), used to represent elevation differences associated with each non- LiDAR DEM. Background: hill-shaded NB-DEM. LiDAR-DEM source: http://www.snb.ca/geonb1/e/DC/DTM.asp.
![]()
Figure 2. A 5-stage workflow developed for the purpose of unifying spatial referencing and resolution, and reducing non-LiDAR to LiDAR elevation differences in bare-earth elevation. The preferred technique for each stage is highlighted in red.
Stage 1. Homogeneity between the non-LiDAR DEMs in terms of the spatial reference and resolution was ensured through re-projection. All downloaded open-sourced DEM tiles were combined into seamless coverages and re-pro- jected (CSRS NAD 1983 New Brunswick Stereographic) in three ways (bilinear, cubic, nearest) for three cell sizes (10 m, default, 90 m).
Stage 2. The Stage 1 data layers resampled at 8, 10, 12, 14, 16, 18, 20 m cell size were randomly point sampled to generate a point shape file containing (n = 86,915) for each layer.
Stage 3. The points of the Stage 2 shape file were re-interpolated using three interpolation methods (Inverse Distance Weighting, Kriging, Nearest Neighbor) using three cell sizes (5, 10, 20 m).
Stage 4. The best-fitted non-LiDAR to LiDAR match was generated through
1) determining which of the Stage 3 generated data layers had the least non- LiDAR to LiDAR elevation differences, by DEM source (done through basic statistics and boxplot comparisons); for this purpose, the bare-earth LiDAR-DEM, with the last-return points originally interpolated to 1 m, was re-sampled at 10 m resolution;
2) determining which combination and weighting of the best Stage-3 non- LiDAR DEMs should be used for the NB-DEM fusion process;
3) obtaining the best-fitted non-LiDAR combination through linear regression analysis, with the best-elevation generating DEM-layers including the fused NB-DEM layer as independent LiDAR-DEM predictor variables. This analysis was done by systematic selection of evenly-spaced points across the non-LiDAR DEMs within the LiDAR coverages, at 1 point per 1 km2, yielding 7255 points for each layer. The output of the regression analysis included the best-fitted coefficients, their standard error of estimate, Student’s t-value (= best-fitted regression coefficient estimate/stand error of estimate), and Pvalue (= the probability that the best-fitted regression coefficient is significantly different from zero).
Stage 5. The best-fitted regression result (hereby referred to as “NBDEM-Op- timized”) was applied province-wide. With this, flow channel and associated cartographic DTW datasets were derived across the province [32] . The results so produced were validated using a 40 by 40 km tile within a newly acquired LiDAR coverage (Figure 1).
The progressive non-LiDAR to LiDAR elevation difference reduction was achieved by selecting the best-fitting Stage 1 and 2 permutations based on layer-to-layer differences in terms of averages, standard deviations, minimum and maximum values, root mean square differences (RMSD), and elevation difference percentages falling within the ±2 m and ± 4 m LiDAR elevation ranges. Also obtained were percentile box plots for the elevation differences, by layer and by stage. The sequence of best choices so generated is highighted in Figure 2. In detail, using default cell sizes and the cubic resampling technique for re- projection generally led to better results than using the bilinear, nearest, and majority re-projection methods (Stage 1). Similarly, randomly extracting re- projected DEM elevation points at 12 points per hectare (Stage 2), followed by IDW re-interpolation (Stage 3) at 10 m instead of 5 or 20 m resolution generally led to better re-interpolation results than using the kriging and natural neighbor methods.
Each of the non-LiDAR DEMs resulting from the 5-stage process and the LiDAR DEMs were used to derive continuous flow channel networks and associated DTW index across New Brunswick. These processes involved using:
1) the D8 algorithm for deriving flow direction and flow accumulation datasets [29] , with stream channels defined to have a minimum upslope flow- accumulation area of 4 ha;
2) a least-cost algorithm to determine the extent of the least upward elevation change away from each flow channel and shore line [32] .
The resulting DTW pattern generally reflects soil drainage patterns, with DTW ≤ 10 cm very poorly drained, 10 < DTW ≤ 25 cm poorly drained, 25 < DTW ≤ 50 cm imperfectly drained, 50 < DTW ≤ 100 cm moderately well drained, and DTW > 100 cm well and DTW > 20 m generally excessively well drained. The distances between the nearest non-LiDAR to LiDAR-derived flow channels for each DEM were conformance tested by plotting their cumulative frequency distributions. The differences in each non-LiDAR versus LiDAR- derived DTW < 1 m coverages were evaluated by determining the extent of false negative and false positive DTW < 1 m areas, respectively.
3. Results
Elevation differences, by DEM Layer. The layer-by-layer non-LiDAR to LIDAR elevation differences are overlain on the corresponding hill-shaded DEMs in Figure 3. Among these, the contour-derived CDED differences are smoothest, the NB-DEM differences are ridged, and the ASTER differences are the most variable. Across the LiDAR-DEM coverages, the differences tend to be largest along forested ridges and valleys (where steep elevation changes occur), and smallest on open areas. Hence, negative SRTM and ASTER elevation differences relative to the LiDAR-DEM occur along valleys, shores, and forest edges, undoubtedly due to the differences in resolution (1m versus 30 and 90 m). Owing to the general terrain and intensive forest cutting patterns across New Brunswick, there are almost as many negative SRTM and ASTER to LiDAR-DEM differences as there are positive differences. As a result, mean differences only have a small positive bias at 1.62 and 1.47 m for ASTER and SRTM (90 m) respectively. The extent of the elevation differences within and across all of the non-LiDAR generated data layers from the 5-stage process are summarized in Table 2. In comparison, the standard deviation and root mean square
![]()
Figure 3. Example for non-LiDAR to LiDAR elevation differences (color-coded) draped over DEM-derived hill-shade for each non-LiDAR DEM. Bottom right: hill-shaded bare-earth LiDAR DEM.
errors are largest for the ASTER DEM, with only 23.2% of its elevation range remaining within the ±2 m LiDAR elevation range, and least so for the best-fitted NBDEM-optimized regression result. For the latter, the remaining elevation differences with respect to the LiDAR DEM are within ± 2 m at 70.1%, and within ± 4 m at 92.9%.
Stage 1 to 4 non-LiDAR to LiDAR elevation difference reductions. The entries in Table 2 and the box plots in Figure 4 concerning the non-LiDAR to LiDAR-DEM elevation differences display a narrowing towards the 0 m difference. The fused and regression optimized DEMs for NB were obtained as follows:
(1)
(weights were chosen to eliminate elevation “ridging” across NBDEM-Fused), and
(2)
(a, b, c, d, e refer to intercept and regression coefficients). The best-fitted regression results so produced are compiled in Table 3. In combination, the influence
![]()
Figure 4. Box plots for the non-LiDAR to LiDAR elevation differences for SRTM90, SRTM30, CDED, NB-DEM, NBDEM-Fused, and NBDEM-Optimized across the DEM optimization extent, showing the 10th, 25th, 50th, 75th and 90th percentiles, the points above and below the 10th and 90th percentiles, and a reference line at 0 m.
![]()
Table 2. Non-LiDAR and LiDAR DEM elevation difference comparison across the DEM optimization extent: statistical summary, in m, including increase in ≤ ±2 to ±4 m elevation difference percentages.
![]()
Table 3. Stepwise regression results with LiDAR-DEM as dependent variable and SRTM90, SRTM30m, CDED, and NBDEM-Fused as predictor variables: regression coefficients and related errors in m; R2 values near 1 across DEM optimization extent.
of stepwise non-LiDAR layer inclusion into the regression analysis generated the following improvement sequence:
.
Hence, the SRTM, CDED and NBDEM-Fused DEMs account for most of the non-LiDAR to LiDAR elevation difference reductions. Entering the selected CDED and SRTM DEMs individually into the stepwise regression process is of further benefit, but only marginally so. Using the ASTER DEM produced no improvements. In detail, the NBDEM-Fused and NBDEM-Optimized layers matched the LiDAR-DEM with the least bias, least minimum, and maximum differences, and with about 70% (i.e., 90th percentiles) of the elevation differences falling within the ±2 m LiDAR elevation range (Figure 4).
Hydrological interpretations (flow channels, wet areas), by DEM layer. Examples of the DEM-produced flow-channel and DTW patterns are presented in Figure 5, by data layer. As shown, the patterns obtained are narrow and widely pixelated with the SRTM90 and SRTM30 DEMs, are somewhat terraced across the ASTER DEM, and are ridged across the NB-DEM. In contrast, the CDED, NBDEM-Fused, NBDEM-Optimized and LiDAR DEMs layers follow fairly consistent patterns. The main hydrological benefit of the Stage 1 to 5 process stems from partially cancelling the elevation differences among the NB-DEM, CDED, and SRTM DEMs, and from essentially eliminating the NB-DEM ridging error.
The 3000 m line scan in Figure 6 (see Figure 1 for location) provides an example of how the non-LiDAR to LiDAR elevation differences change across a forest management area in relation to a recent surface image and the corresponding hill-shaded LiDAR DEM. Also shown are: (i) the 0 to 1 m cartographic
![]()
Figure 5. Flow-channel and cartographic depth-to-water patterns overlaid on the hill- shaded non-LiDAR and LiDAR DEMs.
DTW pattern, shaded dark to light blue from, 0 to 1 m, respectively, and (ii) the DTW and LiDAR-derived canopy height profiles along the scan line. These profiles correspond closely with the image-recognized and hill-shaded DEM features, e.g., low tree height across recent cuts and wetlands, incised flow channel locations, elevated or incised roadbeds, and transmission corridors. The same, however, does not apply across the non-LiDAR elevation difference profiles. For these, the ASTER and CDED profiles have the widest excursions while the SRTM and NB-DEM profiles reflect―in part―the forest height variations, except along recent image-revealed forest cuts. The SRTM to LiDAR elevation differences in Figure 6, however, only amount to about one-half to one third of the LiDAR-generated tree height maxima.
Validation results regarding flow-channel and wet-area mapping. The DEM validation results, done outside the Stage 3 and 4 processing areas
![]()
Figure 6. Non-LiDAR to LiDAR elevation differences along a 3000 m scan line (Figure 1) related to the DTW pattern overlaid on recent surface image (bottom) and hill-shaded LiDAR DEM (top). Also shown, scan line for LiDAR-generated DTW and vegetation heights.
(Figure 1), are compiled in Table 4 and Table 5. For this tile, all layers are upwardly biased, but least so for NBDEM-Optimized at 0.58 m, with standard deviation and root mean square differences also being least at 2.84 and 2.90 m, respectively. The ≤ ±2 and ±4 m percentages of the elevation differences remained about the same overall across the DEM optimization and validation extents, but dropped somewhat from 70.1% to 60.5% for ≤ ±2 m, and increased from 92.9% to 95.4% for ≤ ±4 m (compare Table 2 with Table 4).
The conformance results for the nearest flow channel distances between the non-LiDAR and LiDAR-derived DEMs in Figure 7 improved by layer as follows
ASTER < NB-DEM< CDED < SRTM90 < SRTM30
There is a similar performance increase from ASTER to the NBDEM-Opti- mized layers in terms of false positive and false negative DTW < 1 m reductions, with the optimized NBNB-DEM layer being closest to the corresponding 10 m re-sampled LiDAR DEM derivation (Table 5).
In summary, the 5-stage process of combining to the original NB DEM with the CDED and SRTM DEMs and subsequently calibrating the result with available LiDAR DEM pieces not only improved but also extended the flow-channel delineation across New Brunswick in a systematic and comprehensive manner. This is illustrated in Figure 8 by way of an example within the validation tile (Figure 1).
![]()
Figure 7. Cumulative frequency of nearest distances between the non-LiDAR- and LiDAR-derived flow-channel networks, across validation extent.
![]()
Table 4. Non-LiDAR and LiDAR DEM elevation difference comparison: statistical summary for validation site, in m, including percentage of elevation differences ≤ ±2 m and ≤ ±4 m.
![]()
Table 5. Reduction in area (as percentage) of false positive and false negative DTW < 1 m areas relative to the 10m LiDAR-derived DTW pattern across validation extent.
4. Discussion
The above analysis demonstrates that the fusion of non-LiDAR DEMs can lead to systematic elevation difference reductions, which can be further enhanced through regression calibration with available LiDAR datasets. The DEM so optimized can then be used to generate province or region-wide hydrographic DEM interpretations that are similar to what can be derived from LiDAR-gen- erated DEMs, at 10 m resolution at least.
The non-LiDAR-DEM uncertainty reductions by way of the 5-stage process appear to be small numerically at a 20% reduction of false positives, and a 10% reduction of false negatives (Table 4). There is, however, a substantial 8 times out of 10 distance-to-flow-channel improvement from ±25 m (NB-DEM) to ±5 m (NBDEM-Optimized) (Figure 7). This, by itself and in view of the illustrations
![]()
Figure 8. NBDEM-Optimized (bottom left) and LiDAR-DEM (bottom right) derived flow-channel, and DTW < 1 m delineations, versus the water courses of the New Brunswick Hydrographic Network, overlain on aerial imagery and hill-shaded LiDARDEM (top right only). Note the correspondence between the DEM-derived flow channels and the riparian vegetation buffers. Location: part of the validation tile in Figure 1.
in Figures 5-7 has led to considerable improvements in comprehensively mapping flow channels and wet areas across New Brunswick. In turn, this has had a positive effect on local and regional operations planning in, e.g., forest management (harvest block layout and access), transportation and trail routing, environmental impact assessment, land transactions, and emergency responses. In part, some of the province-wide benefits accrue from improved DEM delineations and determinations of upslope basin areas. These areas are used to estimate, e.g., potential stream discharge rates and required culvert and bridge dimensions at any road-stream crossing, by extreme weather events (e.g., 100 mm or stream discharge per day, roughly equivalent to a 100-year storm event). In this regard, DEM error reductions are especially important across flat terrain, where minor elevations changes (artificial, natural, or man-made) can lead to substantial errors in determining storm water flow directions.
The optimized NB-DEM layer is also useful for checking province-wide wetland coverages. In general, already delineated wetlands fall into the DTW< 0.5 m range [30] . Extending the conformance checking between DEM-delineated and actual wetland border using the original NB-DEMlayer amounted to ±80 m, 8 times out of 10. Air-photo delineated wetland borders generally conform to distances within ±40 m, 8 times out of 10. In contrast, LiDAR-DEM derived DTW = 0.5 m contours conform to actual wetland borders within ± 10 m, 8 times out of 10 (details not shown). The optimized DEM, while somewhat less accurate than the LiDAR-DEM, can nevertheless be used to trace and survey-check wetland-to-wetland connections across the province.
The fact that DEM fusion leads to DEM improvements in terms of vertical and lateral error reductions has been demonstrated repeatedly (e.g., [37] [38] [46] [47] ), but the fusion processes pertaining to cell re-sizing, re-projection, re- sampling re-interpolation, DEM weighting and noise filtering all vary. For example, [40] produced “a nearly-global, void-free, multi-scale smoothed, 90 m digital elevation model” called EarthEnv-DEM90 (http://geomorphometry.org/content/earthenv-dem90). Tran et al. [48] fused ASTER with SRTM30 data through (i) DEM quality assessment and pre- processing, (ii) hydrologic DEM enforcement, (iii) void filling and projection shifting, (iv) DSM versus DTM bias elimination by landform, and (v) DEM de- noising. The 5-stage process in this study varies from [48] by way of systematic cell size and DEM re-processing (re-projecting and re-interpolation), and using the regression process for bias removal. The reference DEMs also differ: using a DEM generated from the 1:10,000 topographic map 5 m intervals and spot- height elevation data (i.e., similar to CDED) versus using 1 m LiDAR DEM re- sampled at 10 m.
In terms of applying the above approach to other areas with similar and/or different DEMs including LiDAR DEM coverages, it is important to examine each re-projected, re-sampled and re-interpolated DEM in reference to artifacts and hydrographic correctness. In this regard, all open water surfaces need to be rendered flat, and streams and rivers need to drop monotonously through their surrounding terrain towards their receiving shores.
5. Conclusions
The systematic fusion of currently available DEM layers for all of New Brunswick not only led to considerable non-LiDAR to LiDAR DEM elevation difference reductions, but also produced a closer and verifiable correspondence between the resulting flow-channels and wet-area derivations. In summary, the 5-stage process as described in this article has shown that:
1) Non-LiDAR elevation differences relative to LiDAR DEMs can be reduced through careful analysis requiring re-projecting, re-sampling, and re-interpolation, followed by selective non-LiDAR DEM amalgamation.
2) The fusion process of the SRTM, CDED and NB-DEM layers, each used at about equal weight, was effective in generating a much improved non-LiDAR- DEM coverage across New Brunswick. Using the ASTER DEM did not improve the best-fitted fusion result so obtained.
3) The fusion process removed many layer-specific artifacts and large layer- to-layer elevation differences, while the non-LiDAR to LiDAR-DEM regression process reduced the NBDEM-Optimized elevation bias to less than 1 m.
4) While the results are specific to and can be applied comprehensively across New Brunswick, similar non-LiDAR DEM improvements could be incurred elsewhere.
Acknowledgements
Financial support for this work was received from NSERC’s AWARE Research Network, the Canadian Wood Fibre Research Centre, Natural Resources Canada, and from the Forest Watershed Research Center at UNB. We thank Bernie Connors (Service New Brunswick) for discussion and comments.