Copyright (c) Susan Laflin. August 1999.

Probably the most familiar example of a contour plot is the method used by the Ordnance Survey to indicate the height of the ground on one of their large scale maps. However, the idea of drawing contours to indicate the shape of a surface is not restricted to this application, and the next figure shows an example of a contour plot of magnetic intensity over an area on an archaeological site.

Example of Contour Plot

In the methods described in this chapter, the data is assumed to be an array of spot-heights giving the values at the mesh points of a rectangular grid. If the "squares" between each group of four values are not in fact square, then you will have to scale the resulting output to extend either the x or the y values by the required amount.

In generating the curves drawn to represent the contour, there are two main approaches. The first one chooses a particular contour and draws it over the whole area. It repeats this process for each contour value in turn. The second approach draws all contours in the first square, and then moves to the next and draws all its contours, and continues in this way until all the squares have been covered. Both methods have advantages and drawbacks and will be discussed here.

Drawing One Contour over the Whole Area.

There are several requirements which must be satisfied by any contouring algorithm to make it suitable for the first approach:

Completeness

Contours over an area

Since the surface will probably have a number of peaks and dips, there will be a number of contours surrounding the peaks which will take the form of several closed loops rather than a single curve. Any method used to follow contours across an area must make sure that all these loops are drawn. This problem does not arise in the "Square-by-square" method since any contour loops lying entirely within a square will not be drawn.

Contour over Spot-height Array

The figure shows an example of one particular contour value c. The dots indicate the positions at which spot-heights have been recorded, and in this case dx and dy are equal, so that the `squares' are in fact square. To ensure that all branches of the contour c have been located, it is necessary to scan the entire array, row by row, and record all the cases where the contour crosses the dotted line joining the two spot-heights. These cases will be indicated by one of the following conditions:

either		h(i,j) less than c  less than h(i,j+1)  
or		h(i,j) greater than c greater than h(i,j+1) 

Form of Interpolation Function

In order to calculate values between the spot-heights, some method of interpolation must be used. The formulae applied to the contouring methods discussed here, make certain reasonable assumptions about the surface, and these influence the shape of the resulting contours. For example, the surface is assumed to be continuous and single-valued , and this is usually correct although it is easy to imagine cases where this is not true. If you have a quarry or mine-shaft with a vertical face, then the surface will have a jump or discontinuity in the values of its height as you move across the edge of the quarry. (If you are carrying a theodolite at the time, you may also have a broken leg and a large repair bill so discontinuities are best avoided. ) Worse still, if the quarry is undercut, there will be a few places where there are three values of the height of the surface for a given (x,y) value. However, these anomalies do not occur very often and are usually ignored in the contouring packages.

In comparing different contouring methods, it is necessary to have a fast and accurate method of estimating the surface height between the data points, in order to calculate the sequence of coordinate pairs which define the contour. It is then necessary to decide whether Polyline output is acceptable or whether some form of smooth curve is necessary.

Different Interpolation Formulae.

If the Polyline output leads to the effect shown in figure a, then some other form of output is needed. However, not all forms of smooth curve are suitable. Figure b shows an example where the curves are fitted to each set of contour values independently, and this leads to contours of different value crossing one another. This is not acceptable if the surface is supposed to be single valued. Figure c shows the case where the interpolation formula is derived from a surface patch. This means that the curves, although smooth are not independent of each other and so cannot cross.

Method to Draw One Contour.

Assume the contour of value c must be followed across the whole area, and that the array has been scanned to locate all the cases where this contour crosses between h(i,j) and h(i,j+1). First trace the contour across the square with h(i,j) at its bottom left-hand corner. Evaluate the mid-point of the square and assume it is either greater or less than c. The six possibilities are shown in the next figure. Cases a and b have three values on one side of the contour and one on the other, and the value of the mid-point determines whether the contour loops round the mid-point or not. Cases c and d have the two left-hand values on one side of the contour and the two right-hand values on the other side of it. The value of the mid-point determines on which side the contour passes. The last two cases occur when the surface has a saddle-point and imply another branch of the contour, shown dotted, in the same square.

Tracing a Contour

The fastest form of method is shown here, where the coordinates of the intersections of the contours with the straight lines are calculated by linear interpolation, and these points are joined by straight lines to approximate the contour.

The contour enters one side of the square. Joining the corners of the square to the mid-point divides the square into four triangles. Each time the contour enters a triangle, two vertices will lie on one side of the contour and the third will lie on the other. Thus you can calculate through which side the contour leaves the triangle and hence which is the next triangle to be entered. This is true for all triangles over the whole area and so the whole contour may be drawn.

If the contour has started at one edge of the area, then the drawing of this branch of the contour will finish when it crosses an edge of the area, either the same edge or a different one. If, on the the other hand, the contour started somewhere in the middle of the area, then if it returns to this point, it will have completed the drawing. The most difficult case to program occurs when the contour starts in the middle of the area and finishes on an edge, with the early part of the contour undrawn. This can be avoided by drawing all contours which cross the edges first, and then filling in any remaining loops.

Different Interpolation Formulae.

If the linear method leads to output such as that in figure a, some other technique will be needed. The problem can be detected by checking the change in gradient from one linear segment to the next, just as in the `smooth' curves discussed earlier. However, in this case a circular or cubic segment could lead to problems as shown in figure b. The appropriate solution to this problem is to use the equation of a surface patch and draw very much smaller linear segments in this area. Possible patches are the bi-linear surface patch or the bi-cubic surface patch.

Square-by-Square Method

The other approach to contouring is to draw all the contours in the first square, then all in the second and continue until all the squares have been filled in. For this method to result in successful output, you must be careful that the ends of the contours match up exactly, and this means that you must avoid excessive roundoff errors in calculating the coordinate values. It has the advantage that the number of points in any Polyline remains small.

This approach is also useful if you intend to use the GKS Fill-area function to shade between the contours, since the number of vertices per square will be considerably less than the total number for two complete contours. Also using the edges of the square in the definition of the Fill-area polygon avoids problems which arise because an annulus does not have a single polygonal boundary.

If the Polyline style of output is not acceptable, then you will need to fit a surface patch to the area and use this equation to generate a larger number of linear segments, so that the resulting output appears smooth. The next figure shows the area between contours c and c+1, ready to be coloured by the Fill-area command, and when this is completed for all the squares the resulting picture resembles that produced by the most detailed form of mosaic output.

Outline of Contours

The two surface patches most commonly used are the bi-linear patch and the bi-cubic patch.

Example of Contour Output.

The bi-cubic surface patch was used to calculate the contours in the above figure. For highly curved surfaces, the high quality of the output justifies the additional calculation needed to produce it.