 |
|
|
||||
| Earth Resources Observation and Science (EROS) |
Significant extensions to an existing finite difference approach to the interpolation of digital elevation models (DEMs) are discussed. The underlying finite difference computational structure is retained, but the minimum curvature interpolation criterion is replaced by a locally adaptive criterion which directly minimises profile curvature, that is curvature of the modelled land surface in the down slope direction. The new method should automatically match landforming processes and hence preserve drainage structure. It should remove biases in the representation of profile curvature which have previously limited applications to the analysis of erosion and soil forming processes. The local nature of the new interpolation condition should also permit greater flexibility in modelling landscape features, including breaklines and cliffs. A statistical analysis of the discretisation errors imposed by representing the landscape with a regular grid DEM leads to a method for optimal smoothing of both point and contour line elevation data. This in turn has led to a simple criterion for matching the spatial resolution of the DEM to the information content of the data, a practical advance toward addressing scale issues in hydrological and environmental modelling.
Topography is a dominant control on earth surface processes. It moderates the spatial distribution of climate, an environmental variable of equal importance in controlling the distribution and productivity of biological systems. Topography also directly moderates the flow of water over and through the earth's surface, in turn moderating soil wetness and soil erosion potential, and influencing soil physical and chemical properties.
The chief limitations of regular grid elevation models (DEMs) appear to lie in not being adaptive to topography with spatially varying complexity and in supporting sometimes over simplistic hydrological analyses. The paper by Gallant and Hutchinson (1996) describes an alternative structural approach to digital elevation modelling. Despite these limitations, regular grid DEMs have become an increasingly popular way of encoding the topography for environmental modelling purposes. They are directly compatible with remotely sensed data sources and complex terrain can be represented provided the DEM has sufficiently fine resolution. This raises the issue of scale which has practical implications for data acquisition and storage and theoretical implications for dependent hydrological and other earth surface models. It also has implications for integration of DEMs with other spatial data sets.
The spatially varying complexity of actual topography and its strongly anisotropic (or directional) behaviour preclude a straightforward analytical approach to interpolation of elevation, especially in view of the very large data sets involved. Input and output data sets may contain over one million points. These difficulties can be addressed with analytical methods, by dividing the region into small adjoining pieces, but problems remain in deciding the size of such pieces, how to join them together and how to address anisotropic behaviour still within, and immediately beyond, each piece.
The finite difference approach developed by Hutchinson (1988,1989) has been shown to address many of these difficulties. Analytic functions are replaced by the regular grid of points and a high quality interpolation criterion related to minimum curvature is used to construct the basic interpolation method. The finite difference formulation can be made computationally efficient, by using a suitable multi-grid iteration method, and facilitates the incorporation of sensible drainage structure by imposing directional constraints on the grid points relating to the direction of surface water flow. The directional constraints can be inferred automatically by a drainage enforcement algorithm and obtained directly from digitised stream line data. The method can also infer directional constraints, defining ridges and stream lines, from the corners in fine scale contour data. With these extension, the method is thus already strongly locally adaptive, although the underlying minimum curvature algorithm is not. Two significant extensions to this finite difference methodology are discussed.
The first of these, which has been developed and already applied to the existing method, involves a statistical analysis of the discretisation errors imposed by representing the landscape with a regular grid DEM. This has led to a locally adaptive method for optimally smoothing both point and contour line elevation data. It has in turn led to the development of a simple criterion for matching the spatial resolution of the interpolated DEM to the information content of the source point and line data. This is intended to be a practical advance toward addressing scale issues in hydrological and environmental modelling.
The second extension, which is still under development, replaces the minimum curvature interpolation criterion by a locally adaptive criterion which directly minimises profile curvature, that is curvature of the modelled land surface in the down slope direction. The new method automatically matches landforming processes and hence preserves drainage structure. It holds the potential to remove biases in the representation of profile curvature which have previously limited applications to the analysis of erosion and soil forming processes. The local nature of the new interpolation condition should also permit greater flexibility in modelling landscape features, including breaklines and cliffs.
The current interpolation method employs a simple multi-grid approach which initially solves the iterative interpolation problem on a very coarse grid resolution. The problem is then solved at succeeding finer resolutions, until the final resolution is obtained. Starting values for each succeeding grid resolution are obtained by bilinear interpolation from the preceding coarser grid.
At each resolution the elevation data are simply allocated to the centre of nearest grid cell. This introduces a small positional error in the data which may be interpreted as a vertical error in the point placed at the centre of the cell. This is illustrated in Figure 1 where the data point A on a sloping terrain surface is allocated to the centre of the grid cell of width h, leading to a horizontal displacement d, which gives rise to a vertical error z.
Figure 1. Vertical error of displaced data point.
The size of this vertical error depends on the slope of the grid cell and the size of the horizontal displacement as shown in Figure 1. If it is assumed that each data point is placed randomly within its associated grid cell then the variance of the corresponding vertical error can be shown to be
var = h*h*slope*slope/12.
The amount of data smoothing can then be determined by weighting the elevation data points, as allocated to the grid cell centres, according to their error variances as estimated by this equation, and determining the smoothing parameter so that the average of the squares of the weighted residuals is 1.0. This matches the expected value of the average of the weighted residuals of the displaced data points from the actual terrain surface.
The amount of smoothing is therefore locally adaptive to the slope of the terrain. Since slope is a function of the terrain surface, it cannot be specified until the terrain surface has been interpolated. The iterative nature of the interpolation procedure is crucial to the imposition of this locally adaptive weighting. It permits the use of steadily improving estimates of the terrain slopes until the process converges for each successive grid resolution. The elevation procedure ANUDEM has been upgraded to apply the locally adaptive data weighting. Care was required in adopting a suitably stable subsidiary iterative process which the smoothing parameter converges to a final value. Further details will be published elsewhere.
As the multi-grid iteration proceeds from coarse DEM resolution to the final specified resolution, the slopes of the fitted grid at the data points steadily increase in magnitude. At coarse resolutions several data points may be allocated single grid cells leading to averaging of the data and resultant smoothing of the fitted DEM in relation to the actual terrain surface. Eventually, the DEM resolution is sufficiently fine for there to be little or no data point averaging, and the slopes of the fitted DEM at the data points stabilise. At this stage, all information has been extracted from the source terrain data, since further refinement of the DEM resolution would not change the DEM slopes.
The African DEM (Hutchinson et al. 1996) was calculated at 10 grid spacings which were successively halved from 6.4 degrees to 0.0125 degrees. Figure 2 shows that the root mean square slope of the DEM at the data points stabilises to a value of 3.5% once the DEM resolution has reached the tenth grid resolution of 0.0125 degrees or about 1.25 km. The point data used in fitting this DEM were digitised from the 1:1M scale air navigation charts. The 1 km resolution used in applications by the USGS of the full 1:1M scale data would appear to be consistent with the information content of this source data set.
Figure 2. Root mean square slope vs resolution of African DEM.Figure 1.
This gives rise to a simple criterion for optimising DEM resolution to the information content of the source data. This has significant practical implications for data inventory, since data storage requirements are very sensitive to grid resolution. The smoothing of the discretisation error can be viewed as providing an optimal representation of the terrain surface at all resolutions up to the final optimised resolution. At each resolution the fitted DEM values are estimates of the average of the actual terrain surface across each grid cell. DEMs produced in this way have facilitated investigations by Gessler et al (1996) of the effects of DEM resolution on derived terrain parameters at various resolutions.
The next step in the development of the locally adaptive method is to develop a locally adaptive roughness penalty to define the interpolation process. This has been achieved until now by applying various localised constraints to an underlying interpolation procedure defined by a general, non-adaptive, roughness penalty. It is intended that these local constraints be incorporated directly into the proposed adaptive roughness penalty. The proposed procedure is to minimise profile curvature of the DEM in the downslope direction.
The definition of this roughness penalty depends on the aspect of the terrain surface being interpolated. As for the adaptive data smoothing criterion, an iterative process is therefore required so that steadily improving estimates of terrain aspect are obtained as the iteration proceeds. It is anticipated that this process will be stable because terrain aspect is less sensitive to DEM resolution than profile slope and curvature. In fact, in the simple but reasonably typical situation of a hill with circular contours, it can be seen that the profile slope and curvature can be prescribed arbitrarily while aspect, as represented by the circular contours, remains constant.
A locally adaptive interpolation criterion also requires substantial additional computation. A prototype implementation has been developed which attempts to compensate for the additional computational burden by making the implementing most of the computations associated with each row of the DEM in a way which can be directly parallelised on a vector processor.
The stability of the overall process is a critical issue. In the current prototype additional stability has been imposed by adding a small factor times the original generalised curvature roughness penalty. An application to the small point elevation data set provided in Table 5.11 of Davis (1986) is shown in Figure 3. The roughness penalty used here was the sum of profile curvature and one tenth of the generalised curvature. The figure shows stronger identification of valley structure than the minimum curvature interpolation shown in figure 1 of Hutchinson (1989). The shape of the DEM is similar to that of the DEM shown in figure 6 of Hutchinson (1989), which was obtained by imposing specific drainage constraints on a modified minimum curvature interpolation procedure.
An preliminary application to contour data has reduced the bias towards the heights of the data contours exhibited by DEMs interpolated using non-adaptive roughness penalties. In particular the application has shown better extrapolation towards peaks defined by surrounding contour line data. This should also reduce biases in representation of profile curvature, as the bias towards the data contour heights using the general roughnes penalty is due to a slight flattening of the fitted DEM in the vicinity of the data contours.
There appears to be trade-off between stability, as currently obtained by adding a generalised curvature roughness penalty, and sensitivity to drainage structure. If the generalised curvature factor is too small then aspect is not stably defined. On the other hand, if it is too large, then drainage structure is too generalised. Further developments in the algorithm are anticipated.
Figure 3. Locally adaptive interpolation of Davis elevation data.
The locally adaptive approach to digital elevation modelling shows promise in matching the spatially variable complexity of actual terrain surfaces. Locally adaptive smoothing of discretisation error has been successfully implemented and has yielded a useful criterion for optimising resolution of fitted DEM to the information content of the source data.
The locally adaptive roughness penalty also shows promise, but problems associated with maintaining stable and accurate definition of terrain aspect during the interpolation process need to be satisfactorily resolved. The parallelisation of the computations helps to reduce the additional computational cost of the locally adaptive technique.
Davis,J.C. 1986. Statistics and Data Analysis in Geology, Second Edition, Wiley, New York.
Gallant,J.C. and Hutchinson,M.F. 1996. Towards an understanding of landscape scale and structure. These proceedings.
Gessler,P.E., McKenzie,N. and Hutchinson,M.F. 1996. Progress in soil-landscape modelling and spatial prediction of soil attributes for environmental models. These proceedings.
Hutchinson,M.F. 1988. Calculation of hydrologically sound digital elevation models. Proceedings, Third International Symposium on Spatial Data Handling, Sydney, Columbus: International Geographical Union, pp. 117-133.
Hutchinson,M.F. 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology 106: 211-232.
|
For further assistance, please contact: Customer Services U.S. Geological Survey Center for Earth Resources Observation and Science (EROS) 47914 252nd Street Sioux Falls, SD 57198-0001 Tel: 800-252-4547 Tel: 605-594-6151 Fax: 605-594-6589 Email: custserv@usgs.gov Business Hours: Monday thru Friday, 8:00 a.m. to 4:00 p.m., central time |
| Accessibility FOIA Privacy Policies and Notices | |
  |
|