Contour interpolation tool - Documentation

This text describes the use of the contour_interpolate tool that implements an interpolation algorithm for calculating a value grid fully filled with values from gridded contour files.

Compiling the program

This tool requires the CImg image processing library and the GDAL library to be built and should be compilable on all platforms supported by these libraries.

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:

1. 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.
2. Based on this first approximation the unknown height contours are drawn into the contour map using the average height along their path.
3. 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.
4. 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.
5. 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.
6. The last three steps are repeated (normally once, could use more for better results) with a decreased amount of smoothing.
7. From this elevation model the local minima and maxima are determined and their heights set to avoid flat hill/mountain tops.
8. 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.
9. 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.
10. Finally the water mask image is applied and the results written to a file.

Program options

the program offers the following command line options:

-i contour-file:water-file:undefined-height-contour-file ...

Specifies an gridded contour file data source. This option has to be specified last and several input sources may be listed. Each consists of three files: The normal contours (pixel values denote height), the water surfaces (pixel values denote height as well) and the undefined height contours (integer pixel values denote the identifying index of the contour, all pixels with the same value get the same height)

-o File:StartX:EndX:StartY:EndY:SizeX:SizeY

Output file (with corner coordinates and resolution specified)

-s Height

Height value to be interpreted as sea level in the input data. Default: 0.0.

-so Height

Height level to assign to all sea surface pixels in the output file. Default: 0.0.

-z

Use this to clip all height values below sea level. Default: off.

-sz 0|1|2

Setting this to values different from 0 activates an enlarging of the sea surface near the coast line. Default: 0.

-e Value

Error bound (in height units) for the interpolation algorithm. Note the technique does not guarantee the error to stay below this value. Too low values can lead to terracing artefacts. Default: 1000.0.

-g Value

Grid stepping for the gradient line calculation step. Default: 4.

-n

Dry run, only draw the input data, no interpolation. Default: off.

-h

show available options

various progress information is printed to stderr during operation. Several additional output images are written (test_*.png) containing debugging information.

Legal stuff

This program is licensed under the GNU GPL version 3. A copy of the license should be included in the program package. If you can not find it you can read it on http://www.gnu.org/licenses/gpl-3.0.html

Christoph Hormann, August 2007