Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
    
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. –  David Carraher 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 –  belisarius Oct 9 '12 at 3:43
1  
@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 :) –  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 –  belisarius Oct 9 '12 at 6:16
show 2 more comments

3 Answers 3

up vote 4 down vote accepted

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

share|improve this answer
    
"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 Narain 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
add comment

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:

enter image description here

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.

share|improve this answer
    
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 Narain Oct 10 '12 at 3:59
1  
@rcollyer "most people don't know that there are other norms besides 2" — true :D –  rm -rf 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 Narain Oct 10 '12 at 19:08
show 1 more comment

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, Apply[Abs[#1 + I #2] &, point - tx]] /@ 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

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.