Data interpolation and ListContourPlot

27

22

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).

contour plot

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!

John Lombardi

Posted 2012-10-09T02:52:52.297

Reputation: 373

Isn't this what Kriging is all about? How do the solutions relate to that concept?

– gwr – 2015-07-04T10:24:02.583

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 – 2012-10-09T03:13:10.850

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 http://i.stack.imgur.com/2YqPb.png

– Dr. belisarius – 2012-10-09T03:43:31.970

2

@belisarius One could try to use the Ingolf Dahl's Obtuse Package for smooth 2D non-grid interpolation (I have not tried).

– Alexey Popkov – 2012-10-09T06:12:28.760

@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 – 2012-10-09T06:14:43.137

@AlexeyPopkov I can't access the link. I may have DNS problems, or perhaps the link is really broken – Dr. belisarius – 2012-10-09T06:16:23.790

@belisarius I just checked again - the link works for me. If you wish, I could upload it somewhere. – Alexey Popkov – 2012-10-09T06:18:22.093

@AlexeyPopkov nah, thanks! I'll just switch to another DNS – Dr. belisarius – 2012-10-09T06:19:54.333

Answers

8

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]]

enter image description here

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]}]

enter image description here

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]

enter image description here

Mark McClure

Posted 2012-10-09T02:52:52.297

Reputation: 31 084

"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. – None – 2012-10-10T21:13:40.883

@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 – 2012-10-10T21:32:21.370

20

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:

contour 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.

user484

Posted 2012-10-09T02:52:52.297

Reputation:

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 – 2012-10-10T03:36:04.870

Thanks for the suggestions. I'll keep MapIndexed in mind for future reference and try to use @ and /@ more often. – None – 2012-10-10T03:59:36.737

1@rcollyer "most people don't know that there are other norms besides 2" — true :D – rm -rf – 2012-10-10T04:17:56.013

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.'s ennui – 2012-10-10T11:43:11.127

@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? – None – 2012-10-10T19:08:25.893

It actually does detect the system is linear (and acts appropriately), but the good thing about directly using LinearSolve[] is that you don't have to generate and carry around the a[k] and b[k]. – J. M.'s ennui – 2012-10-10T23:48:37.470

14

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]}]

contour plot of polyharmonic spline


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 = Internal`PrecAccur[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.

J. M.'s ennui

Posted 2012-10-09T02:52:52.297

Reputation: 115 520

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. – Aaron Bramson – 2016-08-15T09:09:37.170

@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[]? – J. M.'s ennui – 2016-08-15T09:50:26.457

@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. – Aaron Bramson – 2016-08-16T04:14:17.437

@Aaron, yes, that's a problem I've found with radial basis functions in general; for a sufficiently perverse point configuration, the underlying matrix becomes quite ill-conditioned. I've seen algorithms that are more stable than the implementation given here, but they are more intricate, and I have yet to find the time to implement them. – J. M.'s ennui – 2016-08-17T16:59:16.097

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.'s ennui – 2012-10-11T04:17:10.413

I have list of data (in .dat file). it has "three columns" and I want to plot it by contour plot. however, the contours are very sharp (the same as fig.1 in this page) . I really want to use the code written by the author of "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)..... Unfortunately it doesn't plot anything!! would you please help me to solve the related problem?? – user65564 – 2019-10-08T07:38:22.793

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 – 2013-03-28T13:39:46.193

3

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]]}]

contour plot of irregular data

My implementation of Mark's answer is below;

colors = ListDensityPlot[data, ColorFunction -> "Rainbow"];
contours = ListContourPlot[data, Contours -> 20,
                           ContourStyle -> {AbsoluteThickness @ 3, Thin, Thin, Thin},
                           ContourShading -> None];
Show[colors, Graphics[{contours[[1]] /. {Line[x___] -> BSplineCurve[x]},
                       AbsolutePointSize @ 5,Point /@ data[[All, 1 ;; 2]]}]]

smooth contour plot of irregular data

Note that this does not preserve the crisp color bands, but it makes a nice contour plot of irregular data while retaining any styling applied to the original plot.

atlasgeo

Posted 2012-10-09T02:52:52.297

Reputation: 141

1This 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] – J. M.'s ennui – 2017-08-24T18:00:45.330

Just curious, but why did you change the Epilog to use Drop instead of the Part brackets? – atlasgeo – 2017-08-24T18:46:00.937

It was just a different way of doing the same thing. You extract the first two, and I remove the last one. – J. M.'s ennui – 2017-08-24T18:53:16.347