Digital Terrain Models - Traps for the Unwary

S.M.Wise
Department of Geography and Sheffield Centre for Geographic Information and Spatial Analysis, University of Sheffield
Email: s.wise@shef.ac.uk



1. Introduction

The storage, display and analysis of data about the terrain surface is arguably one of the most widely used areas of GIS functionality. The characteristics of the terrain - its elevation, gradient and aspect - have a fundamental control on many aspects of environmental processes and also play a significant role in controlling Man's activities (Weibel and Heller 1991, Moore et al 1991). The terrain is handled in GIS using the Digital Terrain Model (or Digital Elevation Model - the distinction between these terms is not important here and for simplicity the abbreviation DTM will be used throughout). The commonest form of DTM is the gridded DTM, in which elevation values are estimated on a regular grid, and these are becoming widely available at all scales from the global to the local, and are commercially produced by National Mapping Agencies in many countries.

One of the attractions of the gridded DTM is that it can be handled using any GIS package which can handle raster data - this is in marked contrast to the Triangulated Irregular Network, for example, which although arguably a superior representation of the terrain surface, requires specialist software for its manipulation. However, as Fisher (1997) has pointed out, raster data layers are not necessarily all the same, and may differ in the way that the values stored in the raster pixels should be interpreted. In particular he comments that '... the pixel ... is a delusion which may become a snare for the unwary given the way it is treated in most modern software'. This is particularly true in the case of the gridded DTM, which is often treated as if it were a raster data layer, with the values held in pixels, whereas it is in fact a series of points located on a regularly spaced grid. This apparently minor difference can lead to DTMs being treated in an inappropriate manner by GIS software. The aim of the present paper is to illustrate some of the problems which can arise from the way DTMs are handled in GIS software and suggest how future work by the GIS community might reduce them.

The first section explores ways in which terrain surfaces are handled in GIS software, elaborating upon the difference between point and area-based representations. The next sections explore two of the most important consequences of treating these as if they were the same - ambiguities with the interpretation of the values in the DTM and uncertainty in identifying the location of those values. The paper concludes with some general observations on the role of the data model in representing the terrain surface and some suggestions for future work.

2. Models for terrain data

The terrain is one example of a more general class of geogrphical phenomena in which the values of the phenomena are continuously variable in space, other examples being temperature, pressure and rainfall. Such phenomena are normally handled as surfaces and, as pointed out by Kemp (1993), surfaces in GIS can be represented in 6 ways:

All except the first are commonly used to store data on terrain elevation. Data attached to irregularly spaced points is a very common format for initial observations on elevation, from field survey for example, but cannot truly be considered a model of the surface itself, because there is no information (either explicit or implicit) about the nature of the surface between the data points. The remaining four methods are all commonly derived from irregular point data. The Triangular Irregular Network or TIN, connects the points together to form a set of triangular facets covering the entire surface, and it is these triangles which define the nature of the surface and allow the estimation of elevation and surface derivatives at any point. The other three methods all involve interpolation from the original data - in the case of contours, by keeping the elevation constant and interpolating location, in the other two cases by defining location and interpolating elevation.

The commonest form of surface representation for elevation is the last - the gridded DTM which is a set of point estimates of elevation on a regular lattice. There are probably two reasons why this form of surface storage has become so dominant.

The first is the simplicity and availability of the array data structure used to store the lattice of values. The array is possibly the single most important data structure in computing. The two dimensional array is provided in every major computer language with a syntax which is similar in most cases, and which resembles the syntax of matrix algebra in mathematics. For example in FORTRAN, to declare an array and reference a value in it:

INTEGER IMAGE(10,10)

IMAGE(1,1) = 1

Software to handle gridded DTM's can therefore be written by even a novice programmer using a language such as BASIC. This is in stark contrast to the complex data structures needed to handle elevation data in a TIN or in a contour-based representation.

The second reason for the undoubted popularity of the gridded DTM is the compatibility between this representation and between the storage of other types of geographical phenomena in pixels in a raster GIS. Superficially at least,a gridded DTM appears to be the same as any other raster layer - a set of values of some geographical phenomenon located on a regular grid. However, there are two quite distinct ways in which such a set of data values may be interpreted, corresponding to the last two forms of surface storage in the list above.

Figure 1: The lattice (left) and pixel models in raster GIS.

Firstly the values may be interpreted as point estimates of the phenomenon in question, as shown in the model on the left of Figure 1. Secondly they may be interpreted as values relating to the whole area of a pixel, as shown on the right of Figure 1. Since both are loosely referred to as raster representations in the GIS literature, there is no agreed terminology to distinguish them. In this paper the terms lattice and pixel respectively will be used to refer to them.

In the lattice model, the location of the values is known, and strictly speaking the layer is not a raster layer at all, but a point layer in which all the points happen to lie in a regular lattice arrangement. However it is far more efficient to store such a set of points using an array data structure than as a series of explicit point objects in a vector data structure. In the pixel model, the value actually relates to the entire pixel area, but for convenience its location is often treated as being in the centre of this pixel. This makes it simple to estimate distances between features, such as roads and buildings, when these are stored in a raster GIS, and it is generally recognised that such distances have an element of imprecision arising from the nature of the raster data.

The problem is that very few GIS packages make any distinction at all between the lattice and pixel models, and those that do, do not always do so in a very consistent or thorough manner. Treating the two as if they were the same leads to two areas of ambiguity:

Interpreting the attribute values. Given a raster layer, how does the user know whether the values in the cells relate to a true point at the centre of the cell or to the whole pixel area? Perhaps a more important question is whether this matters or not (Fisher 1997).

One situation where the distinction may be important is in the case of DTMs at varying scales. It is very clear, that with a DTM produced from large scale maps or air photos, the values are estimates of elevation at a point (i.e. a lattice). Although the same interpolation procedure may be used to produce global scale DTMs, from much smaller scale maps (e.g. Hutchinson 1993), it may be more appropriate to treat the resulting DTM as a pixel model due to the degree of uncertainty associated with the elevation estimates. In most cases the original source data are themselves not accurate measurements of elevation at known locations but have been generalized to show broad trends. For example, a map at 1:1000000 scale will certainly show areas of upland and their general configuration, but the contours cannot be treated as accurate estimates of height. The same is true of maps at larger scales of course, but I would suggest that the relative error is far less in such cases, because less generalisation will have taken place. The nature of the source data means there will be far more uncertainty attached to the elevation estimates, and one way of interpreting this is as a degree of uncertainty in the location of the estimate.

This has implications for the interpretation of surface derivatives (e.g. gradient and aspect) calculated from global scale DTMs. Again, these are calculated using the same algorithms as with local scale DTMs, but it is not clear that a gradient calculated on the basis of 9 lattice points spaced 1km apart can be interpreted in the same way as one calculated from points spaced 100m apart. Nyquist sampling theory tells us that to sample any phenomenon which varies in space or time, it is necessary to use a sampling interval less than half the wavelength representing this variation. With many natural terrains, points at a spacing of 1 km may well fail this criterion, and will certainly fail to detect all but the broadest trends in the terrain. The calculation of aspect is less problematical - as Hutchinson (1993) points out aspect is likely to remain fairly constant whatever resolution is used for its calculation.

There is thus an argument which says that at global and local scales, the values in a DTM should be regarded as pixel rather than point based - i.e. broad estimates of elevation within the pixel, rather than estimates at the central point.

Identifying the location of the value. In the pixel model, the geographical extent of the layer is very clearly defined as the outer boundary of the set of pixels. Therefore, as long as the coordinates of this boundary are known, then the position of each pixel is known - both its areal extent and its centre. However, in the lattice case, the relationship between the point locations and the geographical extent of the layer is more ambiguous - strictly speaking, since the lattice model consists of points rather than areas, the two need not be related at all. As with a vector point model, it would be possible to have a lattice of points located within an irregular area boundary. This is a purely fictitious example, but illustrates the fact that the relationship between the point locations and the layer extent is one of convention rather then being predetermined. This is important because there are in fact two distinct conventions in common use as illustrated below:

Figure 2: Location of points in the lattice model

In the first example, the points are considered to lie at the centre of the 'pixels', so that the location of the lower left point is offset by half a pixel width from the lower left corner of the geographical boundary. In the second case, the geographical boundary of the layer is taken to be the same as the geographical extent of the points themselves, so that the lower left point is located on the lower left corner of the geographical boundary.

The difficulty is that in most GIS systems the raster model has been developed to deal with pixels rather than points, and therefore no distinction is made between the geographical extent of a layer, and the location of the values stored in that layer - it is assumed that one is defined by the other.

3. DTM error assessment

With these general points in mind, let us now consider some examples of how these issues can lead to genuine difficulties in creating and using DTMs. These are all drawn from the author's own work on the effect of interpolation errors on the use of DTMs (Wise 1998) and while this is not the focus of this paper, it is helpful to briefly describe this research to provide a context for the remainder of the paper.

A good deal of work has been done to consider the accuracy of DTM creation methods, and most agencies publish figures reporting the accuracy of their DTM products, However, the focus of attention is normally the accuracy of the elevation values when compared with some independently measured values. The normal practice is to use a set of spot height measurements, which have not been used in the interpolation process, and to calculate the rmse of elevation:

[Equation 1]

However, this is a very incomplete measure of accuracy for a number of reasons

  1. It is normally based on a small sample of points.
  2. Small errors in elevation can produce large errors in surface derivatives, such as gradient and aspect estimated from the DTM.

What is far more important is the effect of interpolation errors on the final results of any analytical operations which use the DTM. Again, some work has begun to look at this issues, mostly with regard to the effect of differences in the algorithms used for the calculation of surface derivatives (Skidmore 1989) or viewsheds (Fisher 1993) although some work has also considered the effects of errors in the elevation values (Lee et al 1992)

Wise (1998) describes the first stages in a project to consider the role both of interpolation artefacts algorithmic differences on the results derived from DTMs. A single data source (a set of digitised contours in the first instance) was used as the input to a series of different interpolation procedures to produce a set of DTMs. These were then compared in terms of the accuracy of their elevation values, by comparison with spot heights in the normal manner, and also in terms of the accuracy of surface derivatives calculated from them. They were also used to undertake a series of analytical operations, such as identfying drainage lines, and determining the watershed of different streams. The main findings were that DTMs with very similar rmse values could produce widely differing results (hence conirming the poor performance of the rmse as a measure of DTM quality) and that very slight differences between DTMs could produce very large differences in the results of analytical operations in certain cases.

In the original work (Wise 1998) a small test area was used for which only paper maps were available to the author, but in later work an area has been chosen for which existing digital datasets are also available, allowing a comparison between these and DTMs created using standard GIS packages. In the course of this work, a number of problems have come to light which have arisen solely from the ambiguities in the way that the gridded DTMs are being handled in the GIS software, and which illustrate the sort of snares which await the unwary DTM user. Before describing these, let me emphasis that the purpose here is not to gripe about particular software packages. The two which will be discussed - IDRISI and ARC/INFO - are probably no worse than many others in the way they handle raster data. The root of the problems described here is that the models of spatial data on which they are based do not distinguish sufficiently clearly between the lattice and pixel representations of geographical data, and the solution to these difficulties lies in the development of more sophisticated models of spatial data.

3.1 Errors in interpolation

The original work (Wise 1998) deliberately used methods of DTM creation which would be easily available to academic researchers needing to create DTMs for their own study areas, rather than using specialist software which might be used by agencies such as the Ordnance Survey. Therefore one of the methods tested was the INTERCON module in IDRISI (Eastman 1992) There are two stages to producing a DTM with INTERCON:

INTERCON is known to be a relatively poor interpolator because of defficiencies in the second stage - it has a tendency to produce elevation values which are highly correlated along the search directions, producing a characteristic striping along these directions in the DTM (Wood and Fisher 1993). However, the first stage in the process is also a source of error because LINERAS assumes that the raster layer being produced consists of pixels. This means that any pixel which a contour passes through is assigned the value of the height of that contour. At the second stage of the process, this rasterised contour layer is implicitly assumed to be a lattice, because in calculating the gradients around each DTM point, the estimation process must assume that the heights are located at the centres of the pixels. Thus even if a contour passes through the corner of a pixel, that height will effectively be assigned to the centre point of that pixel. It is very simple to demonstrate the sensitivity of this to the location of the contour lines. Figure 3 shows the result of rasterising two straight lines using LINERAS in IDRISI. In the lower case, the upper end of the line has been displaced to the right by a distance which is 0.05 of the pixel width, so that it just passes through the corner of the upper right hand pixel, causing the dislocation of five of the pixels representing the line.

Figure 3: Line rasterisation using LINERAS.

The effect is to produce an error in the estimate of height for pixels which are apparently on the line. The top right hand pixel in Figure 3 shows the extreme case in which the error in height will be:

[Equation 2]

Since the minimum height error possible is zero, the mean error will be approximately half of the value in equation 2. Thus in the sample DTM used in the initial studies with 10m pixels and an average slope of 12%, the mean error in the height estimates from the rasterisation process alone would be 0.4 m. For comparison the RMSE for the final DTM was 1.3m.

This illustrates one possible consquence of interpreting the same set of values firstly as a pixel and secondly as a lattice. A series of other problems all had their roots in the uncertainty about the location of the lattice points.

3.2 Location of lattice points

The comparison between the different DTMs was undertaken using GRID, the raster GIS module of ARC/INFO. The documentation for GRID does in fact point out the distinction between pixel layers (which it calls grids) and lattices (ESRI 1995) but comments that exactly the same data storage is used for both.

Like most GIS packages, GRID stores information about the geographical boundaries of each layer (called the extent of the coverage in ARC/INFO terminology) and the size of the pixels (or equivalently the spacing between lattice points). However it does not explicitly report the location of the pixel centres for grids, or of the points for lattices, and this can lead to a good deal of uncertainty when using the software.

Two DTMs were available from outside organisations, both produced from contours digitised from the Ordnance Survey 1:50000 map series.The first was a section of the Ordnance Survey PANORAMA dataset, the second a section from a UK DTM produced by the Institute of Hydrology (Morris and Flavin 1990). Both covered an area of 10x10km, with a lattice spacing of 50m, but with the OS data the lattice points were positioned on the corners of the 50m pixels, giving 201 points in X and Y whereas with the IH DTM they were positioned in the centres giving a 200x200 lattice. This made it impossible to make a direct comparison in the height values between the two DTMs since the lattice points were 35m apart. However, out of interest it was decided to see what action ARC/INFO would take if a direct comparison were attempted. In GRID the difference between two raster layers can be calculated using a simple expression:

diff = os - ih

where:

diff
name of layer containing height differences.
os
name of layer containing OS DTM.
ih
name of layer containg IH DTM.

ARC/INFO performed this operation without a murmur! The resulting layer, diff had 201 points in both X and Y, but the last row and column were completely filled with missing value indicators. Close inspection revealed that each value in the 201x201 DTM had been subtracted from the value for the point above and to the right from the 200x200 layer - in the rightmost column and northernmost row there were no suitable values from the 200x200 layer and so a missing value was reported.

The action which ARC/INFO has taken is not unreasonable in itself, but what is worrying is that the user is given no indication that anything is amiss. Essentially the points are being treated as if they were pixels - if two pixels overlap they can be regarded as being co-located. This treatment of the location of lattice points has lead to a persitent problem of ensuring that the various raster layers being used in the error assessment work have lattice points/pixels which are located in the same place. In trying to compare DTM interpolation methods it was clearly important that the DTMs were all produced with points which were in the same positions - otherwise, this itself would be a source of difference between the DTMs.

As has been shown, INTERCON uses a pixel-based approach, and so the elevation estimates are at the centres of the pixels (i.e. giving a 200x200 lattice). In ARC/INFO two methods were used produce a gridded DTM from the original contour data (Wise 1998):

Both TINLATTICE and TOPOGRID allow the user to specify the origin and spacing of the lattice points, and both provide defaults, but they differ considerably in how this is done. In all cases the origin of the lattice is considered to be the lower left point.

In TINLATTICE, the default is to place the origin on the lower left corner of the geographical extent of the TIN layer, giving a 201x201 DTM. To produce a 200x200 DTM it was necessary to calculate the positions of points displaced from the lower left corner and the upper right corner by 25m and specify these to TINLATTICE.

TOPOGRID also uses the boundary of the input layer (in this case the set of digitised contours), placing the origin of the lattice on the lower left corner. However the output layer will have a geographical extent which is larger by half a pixel width on all four sides than the input layer - this makes it seem different from the DTM produced by TINLATTICE although it is in fact also a 201x201 lattice. To produce a 200x200 DTM, the XYZLIMITS command must be used within TOPOGRID. Although described in the documentation as defining the limits of input data which will be used, its effect is to control the placing of the lattice points on the output layer - therefore the same positions as specified to TINLATTICE are supplied to XZYLIMITS.

The methods for ensuring that lattice points are where the user wishes them to be are therefore cumbersome and differ from command to command. More seriously, there is no simple method to determine just where the points have been placed, since the relationship between the geographical extent of a layer and the position of the lattice points within that layer depends on the commands used to produce the lattice.

It is also possible to illustrate the effect these apparently small differences can have. Part of the test of the different DTMs involved the delineation of the watershed of various streams in the area using the different DTMs. In one of the tests, the stream for which the watershed was calculated was taken from a vector layer digitised from the 1:50000 map. (This is not necessarily good practice, but made it very simple to extract particular tributaries and test each separately).

The procedure used was to extract the vector line defining a tributary, convert this to raster using LINEGRID in ARC/INFO and then specify this as the target for the WATERSHED function in ARC/INFO GRID. Because the DTMs differed in the location of the lattice points, care had to be taken once again to ensure that the stream was rasterised to the correct lattice location and resolution.

Figure 4: Watershed of Harley Brook (highlighted on the left) determined using the same DTM, but slightly different rasterised versions of the brook.

The consequences of slightly altering the rasterisation of the stream can be seen in Figure 4. which shows two of the watersheds predicted for the small tributary stream in the south east corner of the map. It can be seen that in one case a fairly plausible watershed is predicted (albeit with a gap produced by a sink in the DTM), while in the other case the tributary is wrongly predicted to 'catch' flow coming down the main stream, with a huge addition to the predicted area drained. The two maps have been produced using the same DTM and the same digitised stream - the only difference is the way the stream has been rasterised.

In Figure 4a, the LINEGRID function within GRID has been used. This allows the user to specify the pixel size of the raster layer, but determines the origin of the raster layer from the geographical extent of the vector layer. Since the stream had been extracted into a separate layer, its geographical extent was not aligned with the lattice locations of either the 200x200 DTM or the 201x201 DTM - the pixel centres for the rasterised stream are thus displaced with respect to both DTMs. In Figure 4b, the LINEGRID command in ARC/INFO has been used which does allow the user to specify the origin of the raster layer - in this case the stream pixel centres are co-located with the DTM lattice points.

4. Discussion

The previous section has illustrated how the problems with the way that GIS software packages handle DTMs can have serious, and largely undetectable consequences for users. It is worth re-iterating that the problems described above are not evidence of weaknesses in particular software packages, but of a weakness in the raster model itself.

In an early discussion of data models in GIS, Peuquet (1984) identified three stages in converting features in the real world to a digital representation in the computer. In the first of these, a data model is constructed, and Peuquet (1984) distinguished the two main data models still recognised in most GIS packages - vector and raster. Once the features of interest have been modelled using vector and raster, then an appropriate data structure is chosen which can be implemented on a computer. This is the choice of the software developer rather than the user, with some systems using a quadtree structure and others opting for some form of simple array for example. The final level of the scheme is the file structure, in which the data structure is implemented on a particular type of computer hardware.

The problem with this scheme is that what is called the raster data model is not really a model at all. The problems described above arose because the software assumed that anything represented by a regular tesselation of values has the same characteristics, when in fact this is not true. In effect a data structure - the array - is being treated as if it were a data model. The paucity of current data models in GIS has been noted by several other workers (Goodchild 1991, Burrough 1992, Jones 1997) although most of the discussion has centred on the failure of current data models to represent certain classes of phenomena, such as temporal change or true three-dimensional forms. The example of the DTM shows that current data models are also failing to properly represent those features which are normally regarded as unproblematical.

It is clear that the long term solution to these problems lie in the development of more comprehensive data models, and one possibility is the concept of field and object models which has been suggested by a number of authors (Goodchild 1991, Jones 1997, Worboys 1995). In a more recent review of data modelling in GIS, Worboys (1995) identifies four stage in the process, distinguishing two types of data model - the application domain model and the conceptual computational model and it is certainly tempting to consider the field/object dichotomy as representing the former, with the vector/raster models as the latter. This is certainly feasible in the sense that vector systems can certainly be used to represent fields or objects. The six representations of fields identified by Kemp (1993)are a mixture of models which could be reprsented in raster (pixels), vector (lattice, contours, TIN, random points) or both (polygons) but there will certainly be differences in the operations which are possible using each representation - the question then becomes one of deciding at which level in the data modelling these differences should manifest themselves.

The short term solution to the difficulties identified in this paper is a better handling of lattice and pixel models in current GIS systems, or rather an explicit recognition that the two are different and need to be treated differently. Until this happens, then all those using DTMs in GIS need to take care that the results they are getting mean what they think they do!

References

Burrough P.A. (1992) Are GIS data structures too simple minded? Computers and Geosciences 18(4), 395-400.

Eastman J.R. (1992) IDRISI version 4.0 : Technical reference. Clark University, Worcester MA.

ESRI (1995) ARC/INFO online documentation. Environmental Systems Research Incorporated, Redlands CA.

Fisher P.F. (1993) Algorithm and implementation uncertainty in viewshed analysis. Int.J.Geographical Information Systems 7(4),331-347.

Fisher P.F. (1997) The pixel - a snare and a delusion. Int.J.Remote Sensing. 18(3), 679-685.

Goodchild M.F.(1991) Geographical data modelling. Computers and Geosciences 18(4), 401-8.

Hutchinson M. F. (1993) Development of a Continent-Wide DEM with Applications to Terrain and Climate Analysis. in Goodchild M.F., Parks B.O. and Steyaert L.T. (eds) Environmental modeling with GIS. 392-399. Oxford, Oxford University Press.

Jones C. B. (1997) Geographical Information Systems and Computer Cartography. Longman, Harlow

Kemp K. K. (1993) Environmental modelling and GIS: dealing with spatial continuity. in Kovar K. and Nachtnebel H.P. (eds) Applications of Geographic Information Systems in Hydrology and Water Resources Management, IAHS Publication 211, Wallingford, 107-115.

Lee J., Snyder P.K. and Fisher P.F. (1992) Modeling the effect of data errors on feature extraction from digital elevation models. Phot.Eng.Rremote Sensing, 58(10),1461-7.

Moore I.D., Grayson R.B. and Ladson A.R. (1991) Digital terrain modelling : a review of hydrological, geomorphological and biological applications. Hydrological Processes 5(1),3-30.

Morris and Flavin (1990) A Digital Terrain Model for Hydrology. Proc. 4th International Symposium on Spatial Data Handling, University of Zurich. 250-262.

Peuquet (1984) A conceptual framework and comparison of spatial data models. Cartographica 21(4), 66-113.

Skidmore A.K. (1989) A comparison of techniques for calculating gradient and aspect from a gridded digital elevation model. Int. J. Geographical Information Systems 3(4),323-334.

Weibel R. and Heller M. (1991) Digital Terrain Modelling. in Maguire D.J., Goodchild M.F. and Rhind D.W. (eds) Geographical Information Systems, Volume 1, 269-297. Longman, Harlow.

Wise S.M. (1998) The effect of GIS interpolation errors on the use of DEMs in geomorphology In S.N.Lane, K.S.Richards and J.H.Chandler (eds) Landform Monitoring, Modelling and Analysis. Wiley, Chichester.

Wood J.D. and Fisher P.F. (1993) Assessing interpolation accuracy in elevation models. IEEE Comp. Graphics and Applications 13(2),48-56.

Worboys M. (1995) GIS - A computing perspective. Taylor and Francis, London.