Kansas Geological Survey, Computer Contributions 23, originally published in 1968
1Pomona College, 2Stanford University, and 3Oberlin College
Originally published in 1968 as Kansas Geological Survey Computer Contributions 23.
Geologists are interested in contouring because it is a traditional and highly efficient technique for presenting information visually. Hand contouring, however, is a subjective procedure, susceptible to all pitfalls and hidden biases of other methods of subjective description. Empirically defined mechanical procedures, such as those used in the accompanying programs, are not necessarily better than other approaches, but they have the advantages of reproducibility and consistency. They also have the very real advantage of extreme rapidity of preparation.
The purpose of the programs is to plot contoured maps with irregularly spaced data, trend surfaces, and residuals. The contour program listed here is based on a square grid, the points of which are determined by second-degree trend surfaces fitted to the nearest stations. The package has been thoroughly tested and used for production work for several years.
Other publications of interest include those by Batcha and Reese (1964); Krumbein (1959); Ojakangas and Basham (1964); and Surkan, Denny, and Batcha (1964).
David D. Pollard and Roger Smith contributed to development of the programs while holding NSF Undergraduate Research Participation grants; some of the programs were developed as part of a project supported by NSF grant GA-686. We are grateful to Dianne Willis of IBM for helping with the coding of subroutine DRAFT.
Computer time was provided generously by Western Data Processing Center, University of California at Los Angeles. From 1965 Pomona College provided time on its IBM System/360 and Calcomp plotter. Part of the off-line plotting was paid for by a grant from the Shell Assists Fund, Shell Companies Foundation.
The first difficulty encountered in automatic contouring is that a plotter can move the pen only in either the up or down position along a straight line. Subroutine PLOT uses the coordinates of the new point to interpolate control points on the most direct path possible. A curved line must therefore be subdivided into straight-line segments, each short enough to give a satisfactory approximation. It is true that rounded contours can be drawn with an analog plotter (such as EAI model 3440) by taking advantage of the momentum of the pen and ignoring some of the computed points. Although this sort of filtering is subjective, the results are pleasing. But this is a technique employed by the operator of an off-line plotter, and the task of the computer programmer is unchanged.
Now, if the contours are to be subdivided into straight-line segments, the area must be dissected into triangles within which the surface is taken as planar. Ideally the data points should be at the corners of equilateral triangles forming a hexagonal pattern; but we have chosen instead a square grid, with four isosceles triangles meeting in the center of each square. Subroutine DRAFT works alternately down and up the columns of squares, determining as it goes the coordinates of the end points of each straight line segment. Because these points lie on the sides and diagonals of the squares, discontinuities in the contours are not possible. The output can be conceived of as a matrix with as many rows as there are straight-line segments, and whose four columns specify the end points (x1, y1) and (x2, y2).
After plotting begins, the matrix can be scanned to see if there is another line to be drawn without raising the pen; if there is not, the pen is lifted to the nearest remaining segment. This logic is incorporated in subroutine SIFT. Because the pen moves at the same speed whether it is up or down, this procedure saves considerable time, and is, of course, essential if the smoothing possibilities of an analog plotter are of interest. If the plotter is off-line, the size of the matrix to be sifted is determined by the available memory and the relative costs of computer and plotter time. But it is a disadvantage to have an on-line digital plotter idle while the sifting process is going on, so the relative sizes of the plotter and sift buffers should be adjusted to keep the plotter usefully employed. The buffer sizes given in the listings are for an IBM System 360/40 with 32K bytes and an on-line Calcomp model 565 plotter.
Data rarely are distributed on a grid, but a great simplification is effected if we can construct a grid that is an acceptable substitute. Provided that points which are close together or coincident have similar values attached to them, a very fine grid can come close to honoring every data point. If this condition does not exist, the data are described as noisy, for there is significant within-stations variance. The trouble is, of course, that the amount of computer time increases fourfold when the side of a grid-square is halved, and there will be many more individual line-segments to be sifted. Moreover, in order that the fine grid may honor every data point, spurious jogs will appear in the contour lines. Geological data usually are noisy, so a certain amount of smoothing is an advantage. If we are trying to contour mountainous terrain on the basis of a small number of first-order triangulation stations, our object being to portray the area as a whole, it is perhaps not necessary or even advisable to insist that the contours should exactly fit each data point. If the distribution had been better, data of lower individual quality might have been preferable. This is a sampling problem that has often been ignored. It should be remembered that each data point has to support the surface over a polygonal area extending outwards towards the adjoining points, and it is a waste of effort to determine elevation to a fraction of an inch at one point when there is an unrecorded change in elevation of a hundred feet only a few steps away.
In constructing the grid we use a least-squares criterion and fit a polynomial surface to the points closest to each grid intersection. This is subroutine SGRID (phase DBMCI2). A plane can be fitted to three points; but then no degrees of freedom remain, the surface cannot rise above or fall below the control points, and meaningless results are obtained if the points are approximately collinear. To avoid this problem a quadratic surface is fitted. Because six coefficients are to be determined we must use at least six points, and because these may be distributed poorly we insist on a minimum of eight! A circle is constructed such that the probable number of points contained within it is ten, and if fewer than eight points are found, the radius is incremented and a second count is made. The number and magnitude of the increments can be changed, but experience shows that if more than one increment is allowed the result is likely to be unsatisfactory. Points are weighted inversely as the square of their distance from the grid intersection, and the surface is fitted. On this basis a value is assigned to the grid point. This value may be unsatisfactory because of the distribution of the data points, but instead of a time consuming test of the distribution we use the arbitrary criterion that no predicted value may lie outside the range of the values used to compute it by more than 20 percent of that range. If for any reason the program is not able to assign a satisfactory value to a grid point, it is flagged 999.9 and is ignored during plotting.
If too few data points are used in the definition of a grid point, the loss of degrees of freedom will obviously result in a poor prediction. If too many are used the result also may be bad. A second-degree surface may be an unsatisfactory fit over the larger area containing more numerous data points. The necessary compromise is not easily achieved, because it depends in part on the density of data points relative to the complexity of the actual surface. If the density or distribution of data points is poor, or the data are noisy, it is best to try for a smoothed surface. In this instance, the points should not be weighted (in DBMCI2 remove cards DBM2-4,5 and DBM2-152 through 164, and add "WT(I) = 1.0" after statement #613, DBM2-174). The process then becomes that of fitting a moving average trend surface of the second degree.
It is sometimes convenient to consider the distribution of a variable as the superposition of local effects on a regional trend, and such data can be analyzed by studying the residuals from a least-squares surface. Areas of geologic interest are sometimes outlined more clearly by residuals than by the data in their original form. If the number of polynomial coefficients is increased, the sum of squares of the residuals will diminish until the deviation from the surface approaches the within-stations variance. The surface then represents the data with the noise eliminated. Many measurements thus are condensed into a small number of coefficients, and the noise level may be estimated even if the within-stations variance has not been measured independently. In 1959, Krumbein described a method of analysis applicable if the spacing of data points is irregular. His program was divided into parts, the output from one being the input for the next, and the largest matrix inverted was the lOx 10 for a third-degree surface. We extended Krumbein's method to include eighth-degree surfaces, and published the program in 1963 (see McIntyre, 1963). It was written as a single unit for operation on the 7090/7094 system at Western Data Processing Center (WDPC). Being in FORTRAN, it used an inefficient method of building the large matrix. We later replaced this portion by a subroutine written in FAP, with considerable savings of space and time. The contouring program was developed to supplement the trend-surface program, but the two were kept separate because it proved advantageous to review the output of the first before using computer and plotter time for the second. For several years both programs were operated at WDPC, and the tapes were sent to service bureaus for off-line plotting. To prevent wasteful processing of bad tapes, status reports were printed while the plotter tapes were being written. This feature is retained in the version given here, although in its present form the program is for an on-line plotter. Pomona College took delivery of an IBM System 360/40 in 1965, and a Calcomp 565 plotter was attached on-line. Our original programs were modified for this configuration. Because the memory of the System 360 was only about one quarter that of the 7094, the program had to be broken into phases that are successively fetched for execution and overlay one another in core.
Read the PDF version (6.4 MB)
Kansas Geological Survey
Placed on web Sept. 5, 2019; originally published 1968.
Comments to webadmin@kgs.ku.edu
The URL for this page is http://www.kgs.ku.edu/Publications/Bulletins/CC/23/index.html