Interpolating elevation data grids from contour lines
This page describes a technique for interpolating terrain elevation data grids
(DEMs) or other 2D scalar data grids from contour lines. The technique as well as the implementation
should be considered as unfinished and experimental - the reason i outline it here is that it has
been used for generating elevation data files that have been
published and i want to
document what steps were taken to generate these files.
The task of generating a continuous model of the shape of the earth surface (or at least a regular grid representation of it) from just a few constant altitude contour lines (and maybe other information like spot heights, river lines etc.) has been a subject of considerable research. If you introduce additional criteria as for example the demand for maximum smoothness it is possible to find something like the best possible interpolation concerning these demands and the boundary conditions defined by the contours. Such an optimization is done by various programs available.
For the specific task of approximating the earth surface however criteria like maximum smoothness of the surface are not valid globally. In general the idea of finding the best interpolation technique for any kind of contours representing any kind of terrain is hopeless. One can only try to create a technique that works reasonably well in a large number of situations.
Most interpolation techniques commonly used only work reasonably well in some situations but fail quite miserably in cases they were not designed for. The system i want to describe probably is not really different from this but it has meanwhile been used to process contours of most of Scandinavia. This means it has been tested on a large variety of terrain forms (mountain regions and narrow fjords in Norway as well as flat regions in Finland for example). It scales reasonably well to be used on large regions (DEM sizes of 7500×2500 pixels have been routinely processed but significantly larger sizes should pose no real problem either). I have the impression the resulting elevation models compare quite well to the contour derived DEMs available elsewhere.
As already said above the technique is unfinished and highly experimental. There are surely a lot of ways the system can be improved. Various 'magic' constants are present in the code that should either be chosen automatically my some analysis of the input data or be user-selectable. A really good interpolation would actually require the program to analyze and 'understand' the terrain structure. The algorithm outlined here contains some elements that could help designing such a 'terrain aware' interpolation technique but the way it currently works is fairly straight away and dump.
The algorithm
The input data for this interpolation system are the contour lines drawn into the DEM grid supplemental grid files contain the water surfaces (which are supposed to be flattened) and contours of constant but unknown heights (usually small lakes). It is fundamentally important that the grid coordinate system does not have a strong distortion in the covered area (i.e. the scale is similar in both grid directions and angles are approximately correct).
The following steps are followed to calculate the heights of the pixels between the contours:
- A first approximation interpolation is calculated using directional search in 32 directions from every pixel and weighting all found contour heights in those directions with the inverse distance.
- Based on this first approximation the unknown height contours are drawn into the contour map using the average height along their path.
- This approximation is slightly smoothed using a gaussian blur filter and then at the contour line pixels the difference to the known heights is calculated. This is a simple approximation of the surface curvature at those points.
- This curvature map forms the basis of a new approximation which is now done using iterative blur filter runs and before each resetting the values on the contour pixels. The result is an interpolated version of the error/curvature map.
- The resulting interpolated correction map is added (slightly diminished by a factor of 0.85) to the smoothed (and therefore slightly incorrect at the contour pixels) initial approximation.
- The last three steps are repeated (normally once, could use more for better results) with a decreased amount of smoothing.
- From this elevation model the local minima and maxima are determined and their heights set to avoid flat hill/mountain tops.
- Furthermore from this approximation gradient lines are calculated following the gradient directions of the elevation model. Along those lines heights are linearly interpolated between the know contour pixels. Density of these gradient lines is selectable by the user. The original contour pixels together with the set minima/maxima altitudes and the gradient lines form the new basis of the final approximation.
- The final interpolation is done the same way as the first (steps 1 to 4) just with different settings and 10 iterations instead of 2.
- Finally the water mask image is applied and the results written to a file.
The implementation
The implementation is like the technique itself - experimental and much improvable. You can download it on a separate page.
Short bibliography
Here is a short list of literature on the matter of generating elevation grids from contour lines. This is no way exhaustive but it should get you started if you are interesting in looking more deeply into that matter.
GMT surface tool.closed source Windows program.