Data interpolation and ListContourPlot

Question

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!

Answer 1

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]

Answer 2

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.

Answer 3

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

Answer 4

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},
                           ContourShading -> None];
Show[colors, Graphics[{contours[[1]] /. {Line[x___] -> BSplineCurve[x]},
                       AbsolutePointSize @ 5,Point /@ data[[All, 1 ;; 2]]}]]

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.