# Data interpolation and ListContourPlot

I am fairly new to Mathematica and I have two quick questions on using it for a Hydrology and Hydrogeology class. One is about data interpolation and interpolating without any data defined in an area. First question: I have a set of data:

    data = {
{875, 3375, 632}, {500, 4000, 634}, {2250, 1250, 654.2}, {3000, 875, 646.4},
{2560, 1187, 641.5}, {1000, 750, 650}, {2060, 1560, 634}, {3000, 1750, 643.3},
{2750, 2560, 639.4}, {1125, 2500, 630.1}, {875, 3125, 638}, {1000, 3375, 632.3},
{1060, 3500, 630.8}, {1250, 3625, 635.8}, {750, 3375, 625.6}, {560, 4125, 632},
{185, 3625, 624.2}
}


where $x$ and $y$ are coordinates in space and $z$ is the elevation of bedrock underground. What I want to do is to have a "smooth" and not jagged contour plot. I created the following code which is extremely cumbersome and inefficient, but works:

ListContourPlot[data,
Contours -> Function[{min, max}, Range[min, max, 2]],
ColorFunction -> "TemperatureMap",
Epilog -> {PointSize[.015],
Point[{{875, 3375}, {500, 4000}, {2250, 1250}, {3000, 875}, {2560,
1187}, {1000, 750}, {2060, 1560}, {3000, 1750}, {2750,
2560}, {1125, 2500}, {875, 3125}, {1000, 3375}, {1060,
3500}, {1250, 3625}, {750, 3375}, {560, 4125}, {185, 3625}}],
Text[GB2, {875, 3450}], Text[GB4, {500, 4070}],
Text[220, {2250, 1330}], Text[221, {3000, 940}],
Text[222, {2560, 1250}], Text[223, {1050, 820}],
Text[224, {2060, 1630}], Text[225, {3000, 1810}],
Text[226, {2750, 2630}], Text[227, {1125, 2580}],
Text[10, {875, 3190}], Text[15, {1000, 3430}],
Text[16, {1100, 3550}], Text[17, {1250, 3700}],
Text[18, {750, 3455}], Text[20, {630, 4125}],
Text[21, {185, 3690}],
Text[Style[626.2, Bold, Medium], {390, 3590}],
Text[Style[628.2, Bold, Medium], {600, 3590}],
Text[Style[630.2, Bold, Medium], {800, 3590}],
Text[Style[632.2, Bold, Medium], {1000, 3640}],
Text[Style[634.2, Bold, Medium], {1200, 3375}],
Text[Style[636.2, Bold, Medium], {1500, 3390}],
Text[Style[638.2, Bold, Medium], {1301, 1669}],
Text[Style[640.2, Bold, Medium], {1301, 1469}],
Text[Style[642.2, Bold, Medium], {1301, 1269}],
Text[Style[644.2, Bold, Medium], {1301, 1100}],
Text[Style[646.2, Bold, Medium], {1160, 1000}],
Text[Style[648.2, Bold, Medium], {1470, 1000}],
Text[Style[650.2, Bold, Medium], {1500, 850}],
Text[Style[652.2, Bold, Medium], {2000, 1070}]}]


The ListContourPlot only uses a linear interpolation for the data. I tried adding in

InterpolationOrder -> 5


Any other order does not really change the contour lines. I had a problem using

ContourLabels -> All


; otherwise I would not manually place text for labeling each contour and similarly for labeling each data point.

My second question: The contour plot does not extend contours out past where points are interpolated (I think).

Is there any way to continue these contour lines and extend them to the graph's boundaries? Normally, one would hand draw these contours and extend them.

Thanks so much!

-
Normally you would smooth the contours with InterpolationOrder set to, say 3. According to the documentation, "With a limited number of points, irregular data is linearly interpolated". So, to get smoothing, you need more data points, it would seem. Interpolating is for between points: you cannot "extrapolate" beyond the given data range by interpolating. – DavidC Oct 9 '12 at 3:13
The "natural" extrapolation isn't going to give the results you may expect ContourPlot[Interpolation[data][x, y], {x, 400, 3000}, {y, 900, 4000}, ColorFunction -> "TemperatureMap", Contours -> 30] results in i.stack.imgur.com/2YqPb.png – Dr. belisarius Oct 9 '12 at 3:43
@belisarius One could try to use the Ingolf Dahl's Obtuse Package for smooth 2D non-grid interpolation (I have not tried). – Alexey Popkov Oct 9 '12 at 6:12
@AlexeyPopkov I knew I've seen that one, but I couldn't find it for linking a reference. Thanks!. Anyway, methinks there are too few points in the the input data to get something good enough. Who knows :) – Dr. belisarius Oct 9 '12 at 6:14
@AlexeyPopkov I can't access the link. I may have DNS problems, or perhaps the link is really broken – Dr. belisarius Oct 9 '12 at 6:16

## 3 Answers

The approach taken by Rahul is very nice, I think. I attempted to use this approach with both Interpolation and FindFit (using a sum of scaled Gaussians). Both of these attempts failed; so I'm certain that it was pretty tricky. Ultimately, though, I think the paucity and irregularity of the data dooms this type of approach.

Another approach that I'd suggest is to use ListContourPlot to get a linear approximation (literally containing piecewise-straight contours) and then to approximate those contours with smooth splines. As we see below, we can do this quite easily, if we're willing to sacrifice color. If you do want color, then we need to lay polygons on top of one another in the correct order, which is a bit of a hassle. In addition, it would be nice if their boundaries didn't intersect, which becomes more and more problematic as the number of contours increases.

Assuming that your data has been defined, here's a code that takes all this into account. Note that it is not entirely automated. Relayering the polygons and getting the colors right took a bit of experimentation.

max = Max[Last /@ data];
min = Min[Last /@ data];

(* Initial approximation and then some re-ordering. *)
initPic = ListContourPlot[data, Contours -> Range[min,max,2]];
initPicLines = Cases[Normal[initPic], Line[pts_] -> pts, Infinity];
initPicLines = Join[Reverse[initPicLines[[1;;3]]],initPicLines[[4;;9]],
{initPicLines[[11]]},initPicLines[[13;;]],{initPicLines[[12]],initPicLines[[10]]}];

(* Set the color for the 17 polygons. *)
Clear[col];
col[3] = ColorData["TemperatureMap"][1];
col[2] = ColorData["TemperatureMap"][0.95];
col[1] = ColorData["TemperatureMap"][0.9];
Do[col[k] = ColorData["TemperatureMap"][1-0.05k],{k,4,15}];
col[16] = ColorData["TemperatureMap"][0.35];
col[17] = ColorData["TemperatureMap"][0.42];

(* A function that smooths and extends the piecewise-straight contours *)
smoothAndExtendWithColor[pts_, {i_}] := Module[
{splineFunction, line, start, fin},
splineFunction = BSplineFunction[pts, SplineDegree -> 2];
line = First[Cases[
ParametricPlot[splineFunction[t], {t, 0, 1}],
Line[lp_] :> lp, Infinity]];
If[First[pts] == Last[pts],
{col[i], Tooltip[Polygon[line],i]},
start = First[pts] - splineFunction'[0];
fin = Last[pts] + splineFunction'[1];
{col[i],
Tooltip[Polygon[Join[{start}, line, {fin}]],i]}]];

(* Put it all together. *)
Graphics[{Opacity[1],EdgeForm[Black],
MapIndexed[smoothAndExtendWithColor, initPicLines],
PointSize[Medium], Point[Most /@ data]},
PlotRange -> {{0, 3000}, {500, 4000}}, Frame -> True,
FrameTicks -> False, PlotRangeClipping -> True,
Background -> ColorData["TemperatureMap"][0.7]]


Again, if you're willing to sacrifice color, then things are easier. Here is the simplest version of such code.

initPic = ListContourPlot[data, Contours -> 8];
Graphics[{Cases[Normal[initPic], Line[pts_] :>
BSplineCurve[pts], Infinity],
PointSize[Medium], Point[Most /@ data]}]


It's a bit more work to extend them, since we need to work directly with the spline functions, rather than spline primitives, but it's still not too bad.

smoothAndExtend[pts_] := Module[{},
If[First[pts] == Last[pts],
BSplineCurve[pts],
splineFunction = BSplineFunction[pts];
start = First[pts] - splineFunction'[0];
fin = Last[pts] + splineFunction'[1];
line =
First[Cases[pic = ParametricPlot[splineFunction[t], {t, 0, 1}],
Line[lp_] :> lp, Infinity]];
Line[Join[{start}, line, {fin}]]]]
Graphics[{smoothAndExtend /@
Cases[Normal[initPic], Line[pts_] -> pts, Infinity],
PointSize[Medium], Point[Most /@ data]},
PlotRange -> {{0, 3000}, {500, 4000}}, Frame -> True,
FrameTicks -> False, PlotRangeClipping -> True]


-
"I think the paucity and irregularity of the data dooms this type of approach." I don't understand this sentence at all. Thin plate splines are a standard technique for scattered data interpolation, and can clearly be seen to work in my answer. Perhaps I should edit my answer to make it clearer that this method is not something I pulled out of a hat and fiddled with until it worked. – Rahul Oct 10 '12 at 21:13
@RahulNarain I mean that, looking at the picture, it seems that there are details that might not be totally justifiable from the data. That certainly doesn't nullify the approach or make it wrong. I upvoted it and, as I indicated in my answer, I think that your approach is quite impressive indeed! It's more global and easily automatable than mine. As such, I see the two approaches as complementary. – Mark McClure Oct 10 '12 at 21:32

One thing you can do is fit a smooth function to the data, and draw the contour plot of that instead. Using the thin plate case of polyharmonic splines (see also this nice article by David Eberly), I get the following plot:

Here is my code. Being fairly new to Mathematica, I am open to suggestions for improvement.

data = {{875, 3375, 632}, {500, 4000, 634}, {2250, 1250,
654.2}, {3000, 875, 646.4}, {2560, 1187, 641.5}, {1000, 750,
650}, {2060, 1560, 634}, {3000, 1750, 643.3}, {2750, 2560,
639.4}, {1125, 2500, 630.1}, {875, 3125, 638}, {1000, 3375,
632.3}, {1060, 3500, 630.8}, {1250, 3625, 635.8}, {750, 3375,
625.6}, {560, 4125, 632}, {185, 3625, 624.2}};

{xs, ys, zs} = Transpose[data];

phi[r_] := Which[r == 0, 0, r < 1, r Log[r^r], True, r^2 Log[r]];

n = Length[data];
f[p_] := Sum[
a[j] phi@EuclideanDistance[p, {xs[[j]], ys[[j]]}], {j, n}] +
b[0] + {b[1], b[2]}.p;
sol = Solve[
Table[zs[[i]] == f[{xs[[i]], ys[[i]]}], {i, n}]~
Join~{Sum[a[i], {i, n}] == 0, Sum[a[i] xs[[i]], {i, n}] == 0,
Sum[a[i] ys[[i]], {i, n}] == 0},
Table[a[i], {i, n}]~Join~{b[0], b[1], b[2]}];

ContourPlot[f[{x, y}] /. sol, {x, 0, 3500}, {y, 500, 4500},
AspectRatio -> Automatic, Contours -> 20,
ColorFunction -> "TemperatureMap",
Epilog -> {PointSize[.015],
Point[{{875, 3375}, {500, 4000}, {2250, 1250}, {3000, 875}, {2560,
1187}, {1000, 750}, {2060, 1560}, {3000, 1750}, {2750,
2560}, {1125, 2500}, {875, 3125}, {1000, 3375}, {1060,
3500}, {1250, 3625}, {750, 3375}, {560, 4125}, {185, 3625}}]}]


Note: My previous implementation was incorrect as it omitted the polynomial terms parametrized by b.

-
I must disagree with @rm-rf on one or two things: Norm[p-q] is pretty clear, and most people don't know that there are other norms besides 2, so I wouldn't worry about that. So, instead of tps@EuclideanDistance[p, #] use tps@Norm[p-#]. As to using Total@Table instead of Sum, use Sum it is faster (by an order of magnitude on my system) because Total@Table constructs the list first and then sums it. It is likely that it is more memory efficient, also. For small lists, it won't matter much, though. – rcollyer Oct 10 '12 at 3:36
Thanks for the suggestions. I'll keep MapIndexed in mind for future reference and try to use @ and /@ more often. – Rahul Oct 10 '12 at 3:59
@rcollyer "most people don't know that there are other norms besides 2" — true :D – R. M. Oct 10 '12 at 4:17
Since the a[j] and b[j] enter linearly, it might be expedient to rewrite the system in matrix format, and then use LinearSolve[]... – J. M. Oct 10 '12 at 11:43
@J.M. Very true. For simplicity I chose to enter the equations in the most direct way possible, and Solve worked fine for this example. Perhaps Mathematica can automatically detect that the system is linear and use a linear solver internally? – Rahul Oct 10 '12 at 19:08

I might as well... here's an implementation of the thin plate polyharmonic splines in Rahul's answer that uses LinearSolve[] under the hood, as well as exploiting the block structure of the underlying coefficient matrix (thus reducing the computational burden):

polyharmonicSpline[data_List, vars : {__}] /; MatrixQ[data, NumericQ] :=
Module[{bb, cofs, ls, lx, n, p, tx, xa, xap, wa, ws, Φ},
{n, p} = Dimensions[data];
If[Length[vars] + 1 != p, Return[$Failed]]; wa = data[[All, -1]]; xa = Drop[data, None, -1]; tx = Transpose[xa]; xap = PadRight[xa, {n, p}, 1]; Φ = Function[r, Piecewise[{{r Log[r^r], 0 < r < 1}, {r^2 Log[r], 1 < r}}, 0], Listable]; ls = LinearSolve[Φ[N[Function[point, Sqrt[Total[(point - tx)^2]]] /@ xa, Precision[data]]]]; ws = ls[wa]; lx = ls[xap]; xap = Transpose[xap]; bb = LinearSolve[xap.lx, xap.ws]; (ws - lx.bb).Φ[EuclideanDistance[vars, #] & /@ xa] + bb.Append[vars, 1]] polyharmonicSpline[data_List, vars__] /; MatrixQ[data, NumericQ] := polyharmonicSpline[data, {vars}]  The routine is designed to work for any$n$-dimensional data. Try it out: f[x_, y_] = polyharmonicSpline[data, x, y]; ContourPlot[f[x, y], {x, 0, 3500}, {y, 500, 4500}, AspectRatio -> Automatic, Contours -> 20, ColorFunction -> "ThermometerColors", Epilog -> {AbsolutePointSize[4], Point[Most /@ data]}]  Using the new experimental function DistanceMatrix[] in version 10.3, here is an update of the polyharmonic spline routine: polyharmonicSpline[data_List, vars : {__}] /; MatrixQ[data, NumericQ] := Module[{bb, cofs, ls, lx, n, p, prec, xa, xap, wa, ws, Φ}, {n, p} = Dimensions[data]; prec = InternalPrecAccur[data]; If[Length[vars] + 1 != p, Return[$Failed]];
wa = N[data[[All, -1]], prec]; xa = N[Drop[data, None, -1], prec];
xap = PadRight[xa, {n, p}, N[1, prec]];
Φ = Function[r,
Piecewise[{{r Log[r^r], 0 < r < 1}, {r^2 Log[r], 1 < r}}, 0],
Listable];
ls = LinearSolve[Φ[DistanceMatrix[xa]]];
ws = ls[wa]; lx = ls[xap];
xap = Transpose[xap]; bb = LinearSolve[xap.lx, xap.ws];
(ws - lx.bb).Φ[EuclideanDistance[vars, #] & /@ xa] + bb.Append[vars, 1]]


In theory, one should be able to implement other radial basis function methods with a similar approach.

-
If you're a bit concerned about the redundant evaluations of EuclideanDistance[], see this for a different way to construct the underlying matrix. – J. M. Oct 11 '12 at 4:17
I actually just looked at this for the first time. The fastest method for constructing the underlying matrix is likely Leonid's method which takes advantage of the vector processing capabilities even though it performs the same number of calculations as Outer`. – rcollyer Mar 28 '13 at 13:39