# 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. Oct 9, 2012 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 Oct 9, 2012 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). Oct 9, 2012 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 :) Oct 9, 2012 at 6:14
• @AlexeyPopkov I can't access the link. I may have DNS problems, or perhaps the link is really broken Oct 9, 2012 at 6:16

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.
– user484
Oct 10, 2012 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. Oct 10, 2012 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. Oct 10, 2012 at 3:36
• Thanks for the suggestions. I'll keep MapIndexed in mind for future reference and try to use @ and /@ more often.
– user484
Oct 10, 2012 at 3:59
• @rcollyer "most people don't know that there are other norms besides 2" — true :D
– rm -rf
Oct 10, 2012 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[]... Oct 10, 2012 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?
– user484
Oct 10, 2012 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. Oct 11, 2012 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. Mar 28, 2013 at 13:39
• I tried to implement both versions, but with just the code here and the data, the points show up but the contour plot doesn't (all white) in v10.4. Rahul's solution works perfectly. Aug 15, 2016 at 9:09
• @Aaron, for diagnostics: after starting on a fresh kernel, right after the definition f[x_, y_] = polyharmonicSpline[data, x, y];, what do you get when you execute something like f[0, 0]? Did you make sure you had both definitions for polyharmonicSpline[]? Aug 15, 2016 at 9:50
• @J.M. I don't know if it was the kernel or if I had missed something, but I copied and pasted everything into a fresh notebook and your code above worked. On my data it took 10 minutes and then reported "Result for LinearSolve of badly conditioned matrix ... may contain significant numerical errors". It also gave a result that was a very poor match of the data for my purposes, but that's probably due to my data's characteristics. Aug 16, 2016 at 4:14

I recently stumbled on this post because I was wanting to contour irregular data but wasn't happy with the 'linear' contours. I wanted to add a comment to Mark McClure's answer, but since this is my first response I'm not allowed to do so. I merely wish to expand on his solution.

His method of replacing the Line[] with a BSplineCurve[] is clever, but there is an easy way to preserve the color. Consider the following;

data = {{1772721, 582282, -3547}, {1781139, 585845, -3663}, {1761209,
581803, -3469}, {1761897, 586146, -3511}, {1757824,
586542, -3474}, {1759248, 593855, -3513}, {1751962,
595979, -3488}, {1748562, 600461, -3495}, {1749475,
601824, -3545}, {1748429, 612332, -3656}, {1747542,
610708, -3631}, {1752576, 610150, -3650}, {1749236,
605604, -3612}, {1777262, 614320, -3984}, {1783097,
614590, -3928}, {1788724, 614569, -3922}, {1788779,
602482, -3928}, {1783525, 602816, -3827}, {1782876,
595479, -3805}, {1790263, 601620, -3956}, {1786390,
587821, -3748}, {1772472, 591331, -3549}, {1774055,
585498, -3580}, {1771047, 582144, -3528}, {1769765,
592200, -3586}, {1784676, 602478, -3866}, {1769118,
593814, -3606}, {1774711, 589327, -3632}, {1762207,
601476, -3666}, {1767705, 611207, -3781}, {1760792,
601961, -3653}, {1768391, 602228, -3758}, {1760453,
592626, -3441}, {1786913, 605529, -3748}, {1746521,
614853, -3654}};


and then plotting data with some nice contour weights and coloring;

ListContourPlot[data, Contours -> 20,
ContourStyle -> {AbsoluteThickness @ 3, Thin, Thin, Thin},
ColorFunction -> "Rainbow",
Epilog -> {AbsolutePointSize @ 5, Point /@ data[[All, 1 ;; 2]]}]


My implementation of Mark's answer is below;

colors = ListDensityPlot[data, ColorFunction -> "Rainbow"];
contours = ListContourPlot[data, Contours -> 20,
ContourStyle -> {AbsoluteThickness @ 3, Thin, Thin, Thin},

• This should work: ListContourPlot[data, Contours -> 20, ContourStyle -> {AbsoluteThickness[3], Thin, Thin, Thin}, ColorFunction -> "Rainbow", Epilog -> {AbsolutePointSize[5], Point[Drop[data, None, -1]]}] /. Line[l_] :> BSplineCurve[l]` Aug 24, 2017 at 18:00