"Interpolation is a black art, and interpolation methods should never be trusted. Careful researchers always interpolate several times using different methods, and choose the result they dislike least." - Roger Bivand

The purpose of this page is to communicate the status of my MSCS Master's Project at RPI Hartford. The project was designed to implement an interpolation algorithm developed by Dr. Wm. Randolph Franklin at RPI Troy. The algorithm is of interest to cartographers and is used to interpolate terrain data from contour lines.

Although Dr. Franklin offers one of the best interpolation algorithms yet developed, the system of linear equations required to solve it for typical cartographic data sets can be quite large, for example one million equations with one million unknowns. The data structure to hold the coefficient matrix for the system can exceed 1000 GB in size.

Fortunately, the matrix is sparse, meaning it is populated with mostly zeros. It is possible to store the data in alternative data structures called sparse matrices whereby only the array indices and values of the non zero elements are stored. This reduces the size of the storage considerably. However, the result is still a very large problem.

My project implemented the algorithm for large input elevation grids of up to 1201 by 1201 elevation postings. A variety of techniques were originally envisioned, including MATLAB sparse solvers, sparse libraries such as SparseLib++, and multigrid iterative methods. Ultimately, a two level iterative approach using the Paige-Saunders LSQR solver was developed. In the process, I discovered the power of iterative solvers, as well as their problems, the details of which may be of some general interest.

I also compared the Frankln approximation algorithm to one commercial and one research-grade contour to DEM interpolator with additional interesting results.

A brief introduction extracted from my project proposal is shown below. My final report is available by clicking on the link above.

1.0 Statement of the problem to be solved.

The goals of this three credit Master’s Project are:

To extend the work described in the Franklin paper entitled ‘Applications of Analytical Cartography’ ^[1]

Demonstrate the feasibility of the algorithm described in that paper when applied to larger datasets approaching commonly used DEM formats (for example input files of 1201 by 1201 elevation grid nodes.)

Compare the results of the algorithm with commercial implementations of contour-to-grid interpolations, for example: RV2, GRASS, PCI, etc.

2.0 Significance of the problem, project justification.

Interpolating information from known data to unknown data points is one of the classical problems of computational cartography and is used widely in many diverse GIS (Geographical Information Systems) applications. One example is the derivation of DEMs (Digital Elevation Models) from topographical contour maps for areas where DEMs are not available. In order to do this, DEMs are sometimes reverse-engineered from contour map data. In this multi step process, contour maps are scanned into digital format if they are not already archived that way. Then the contour line elevation data is separated from the other information on the map. The contour line data is then conditioned (to the degree possible) to correct bleeding, breaks, spurs bridges and other anomalies.

The contours must then be tagged with their elevation data in a machine-readable format. This is typically done by vectorizing the contour lines and tagging the vectors with their elevations either automatically, or more commonly using a semi-automated labor-intensive process. The tagged vector data is then rasterized and transferred to grid using an interpolation algorithm. Finally, the gridded elevation values are written to some type of GIS format that can be used by other applications. The interpolative step in the process is the focus of this proposal.

One algorithm for the interpolation (or rather the approximation) of known elevations to a regular grid is described by Dr. Franklin in the previously cited paper. The algorithm is based on a novel application of the Laplacian PDE (partial differential equation) whereby a system of over determined linear equations is formulated and solved by performing a linear least squares fit to the known data and unknown grid nodes.

The algorithm has several advantages as compared to previous algorithms:

It does not require continuous contours; i.e. it can deal with breaks, such as those that commonly occur when contours are too closely spaced.

It can use isolated point elevations if known, such as mountain tops.

If continuous contours with many points are available, it can use all the points (without needing to sub sample them).

Information flows across the contours, causing the slope to be continuous and the contours not to be visible in the generated surface.

If local maxima, such as mountain tops, are not given, it interpolates a reasonable bulge.  It does not just put a plateau at the highest contour.

Unlike some other methods, it works even near kidney-shaped contours, where the closest contour in all four cardinal directions might be the same contour.

It has the advantage of producing a surface virtually free of the negative artifacts associated with other interpolation algorithms, for example terracing and ringing.

It also offers the advantage of allowing the surface smoothness and fit to be adjusted by selective weighting of the equations corresponding to the known elevations, thereby allowing grid adjustments to be made on a node-by-node basis if desired.

Although the algorithm presents advantages, it presents computational challenges. If the elevation grid is N by N, the solution is formulated in terms of a N² by N² coefficient matrix. Since the number of flops for row reduction for an N by N matrix is approximately 2N³/3, the time complexity for an N² by N² matrix by the direct solution is 2(N²)³ /3 ~ N⁶. Because the proposed system is overdetermined, additional matrix multiplications are required to solve the system by least squares, increasing the computational time and storage requirements even more.

This algorithm has been implemented by Franklin using MATLAB on data sets of the order of 257 by 257 nodes using sparse matrix techniques. Although this was sufficient to demonstrate the feasibility of the algorithm, real data sets are larger. It is proposed to attempt to demonstrate the algorithm on realistic data sets, as large as 1201 by 1201 grid postings. The task is made difficult by the rapid increase in problem size as a function of the number of input elevation grid nodes. It is further proposed to compare the quality of the output grids produced by the algorithm with other algorithms commercially implemented. In doing so, we hope to offer Dr. Franklin’s useful algorithm to the cartographic community.