How to make a topographic map from sparse sample data?

Question

I need to make a topographic map of the area for which I have only fairly rare samples (x, y, height) of data. Obviously, I cannot make an absolutely accurate map, but I would like it to be “smooth” in a sense. I need to quantify "smoothness" (possibly the inverse of the average square squared surface curvature), and I want to minimize the objective function, which is the sum of two values:

• Surface roughness
• The average square distance between the height of the surface at the reference point and the actual measured height at this point

Since I really want this to be a topographic map, I'm really looking for a way to build contour lines with a constant height, and there might be some kind of smart geometric way to do this, not to mention surfaces. Of course, I want the contour lines to be smooth as well.

Any suggestions are welcome. I hope this is a well-known numerical problem. I am very comfortable in C and have a good knowledge of FORTRAN. About Matlab and R I'm pretty clueless.

---

Regarding where our samples are located: we plan to approximately equal spacing, but we will take more samples where the topography is more interesting. So, for example, we will consider the mountains more densely than the plain. But we definitely have a choice to sample, and we can even take samples if this simplifies the situation. The only problems are

• We do not know how many places we will need to match in order to find the functions we are looking for.

• Sampling is moderately expensive, on the order of 10 minutes. Thus, sampling a 100x100 grid can take a lot of time.

Answer 1

Kriging interpretation can be useful for smoothly interpolating your sparse samples.

Answer 2

R has many different related tools. In particular, look at the spatial view . A similar question was asked in R-Help before, so you can look at it .

Take a look at the contour functions. Here are some details:

 x <- seq(-3,3) y <- seq(-3,3) z <- outer(x,y, function(x,y,...) x^2 + y^2 ) 

The initial schedule is somewhat rough:

 contour(x,y,z, lty=1) 

Bill Dunlap suggested an improvement: “It often works better to fit a smooth surface to data, evaluate that surface on a finer grid and pass the result to the path. This ensures that the path lines do not intersect and tends to avoid false loops that you could get by smoothing the contour lines themselves. Thin plate splines (Tps from the library ("fields")) and loess (among others) can match the surface. "

 library("fields") contour(predict.surface(Tps(as.matrix(expand.grid(x=x,y=y)),as.vector(z)))) 

This leads to a very smooth graph because it uses Tps() to pick data first and then calls contour . It looks like this (you can also use fill.contour if you want it to be shaded):

For the graph, you can use either lattice (as in the above example), or in the ggplot2 package. Use the geom_contour() function in this case. An example can be found here (ht Thierry) :

 ds <- matrix(rnorm(100), nrow = 10) library(reshape) molten <- melt(data = ds) library(ggplot2) ggplot(molten, aes(x = X1, y = X2, z = value)) + geom_contour() 

Answer 3

For an excellent overview of the algorithm, you may need to stitch the surface first to interpolate it onto the grid.

Answer 4

maybe you can use:

• Geomap
• geomapdata
• gtm

from

• Matrix
• Sparsem
• slam

in R