20
$\begingroup$

I have some data of type {{x1, y1, f1}, {x2, y2, f2}, ...} and want to do ListContourPlot. However, there is a problem that the data f is not smooth enough and has small errors.

Here is a toy example illustrating the problem:

data = Flatten[ 
   Table[{i (1 + RandomReal[0.1]), j (1 + RandomReal[0.1]), 
     i + j}, {i, 0, 100}, {j, 0, 100}], 1];
ListContourPlot[data, Contours -> {80, 120}, ContourShading -> None]

enter image description here

I would need to smooth the lines in the above curve, to get something like

enter image description here

Any ideas to get this kind of smoothed plot? Thanks!

And here is the figure from real data (which is too large to paste here), which I need to smooth:

enter image description here

I find this thread InterpolationOrder for ContourPlot related. But I was unable to get the method work in my case because of (I guess) different input.

Note added:

Thanks a lot for the answers! Those answers works great for the toy example, but for the real data I still cannot get smooth lines so far (edit: so far means before the great guys update their answers ^_^. Now it worked great). Here is my data and plotting function, just in case you would like to give a try. I will also investigate why those interesting methods fail to work when I apply it to data...

https://www.dropbox.com/s/2pgobgrhbo42f4v/contourData.dat

data = << "contourData.dat";
ListContourPlot[data, Contours -> {2.30, 6.18}, 
 ContourShading -> None]
$\endgroup$
2
  • $\begingroup$ I guess the main issue is that your toy example is uniformly sampled, while your real data is very far from it (see ListPlot[Most /@ data, AspectRatio -> 1]: i.stack.imgur.com/DgAAh.png vs. i.stack.imgur.com/zfE5C.png). Also, the aspect ratio of the data is very different from the way you're plotting it. $\endgroup$
    – user484
    Apr 7, 2014 at 17:21
  • $\begingroup$ @RahulNarain : Thanks! Yes it's a very good point. The real data is from Markov chain Monte Carlo (MCMC), which is by purpose not uniformly distributed. $\endgroup$
    – Yi Wang
    Apr 7, 2014 at 19:43

3 Answers 3

24
$\begingroup$

N.B. Your actual data calls for a more sophisticated approach than the quick hack in my original answer, so I've replaced it with a much better and quite general solution.


There are two things that make your actual data harder to work with than the toy example. First, it is highly irregular and nonuniformly distributed:

ListPlot[Most /@ data, AspectRatio -> 1]

enter image description here

and second, it has an awful aspect ratio:

ListPlot[Most /@ data, AspectRatio -> Automatic]

enter image description here

If your $x$ and $y$ axes aren't actually spatial coordinates but represent independent quantities with different units, you would do well to rescale the data so that, say, the variances along both axes are equal:

sd = StandardDeviation[data];
{scalex, scaley} = 1/Most[sd];
scaledData = {scalex, scaley, 1} # & /@ data;

Now to smooth an arbitrary nonlinear function described by noisy samples at irregularly scattered locations, I believe a local regression (LOESS) model is appropriate. In fact, LOESS is a generally useful thing so it's worthwhile to have an implementation available in Mathematica, and since it's pretty simple to implement, I went ahead and did it. My implementation follows Cleveland and Devlin's 1988 paper, i.e. local quadratic regression with tricube weights, except their $q$ is my $k$.

(* precompute spatial data structure because we'll be needing nearest neighbours a lot *)
nearest = Nearest[scaledData /. {x_, y_, z_} :> ({x, y} -> {x, y, z})];
(* local quadratic regression with k neighbours around point (x, y) *)
loess[nearest_, k_][x_, y_] := Module[{n, d, w, f},
  n = nearest[{x, y}, k];
  d = EuclideanDistance[{x, y}, Most[#]] & /@ n;
  d = d/Max[d];
  w = (1 - d^3)^3;
  f = LinearModelFit[n, {u, v, u^2, v^2, u v}, {u, v}, Weights -> w];
  f[x, y]]

Now we can plot the contours of the regression function, scaling the coordinates back to the original data:

fit = loess[nearest, 100];
{{xmin, xmax}, {ymin, ymax}, {zmin, zmax}} = {Min[#], Max[#]} & /@ Transpose[data];
ContourPlot[fit[scalex x, scaley y], {x, xmin, xmax}, {y, ymin, ymax}, 
 Contours -> {2.30, 6.18}, ContourShading -> None]

enter image description here

Seems to work pretty well.

enter image description here

$\endgroup$
3
  • $\begingroup$ Isn't this sort of MovingAverage? $\endgroup$ Apr 7, 2014 at 16:33
  • 2
    $\begingroup$ @RahulNarain I would suggest to add smouthData[n_] and this gif for fanciness purposes ;o) $\endgroup$
    – Öskå
    Apr 7, 2014 at 16:43
  • $\begingroup$ The result is amazing! I will read it carefully and try to understand the algorithm! $\endgroup$
    – Yi Wang
    Apr 8, 2014 at 7:24
10
$\begingroup$

This is likely a consequence of taking the toy example too seriously, but LinearModelFit seems like a good choice:

lmf=LinearModelFit[data,{x,y},{x,y}]

ContourPlot[lmf[x,y],{x,0,100},{y,0,100},Contours->{80,120},ContourShading->None]

enter image description here

For a the data provided you might get some use from:

basis[n_,m_] := Flatten[Table[x^i y^j,{i,0,n},{j,0,m}],1]

lmf = LinearModelFit[data,basis[4,6],{x,y}];

Now use this linear model:

cp = ContourPlot[lmf[x, y], {x, 0, 0.5}, {y, -4, 4}, 
  Contours -> {2.30, 6.18}, ContourShading -> None, 
  ContourStyle -> Red]

Show[ListContourPlot[data, Contours -> {2.30, 6.18}, 
  ContourShading -> None], cp]

enter image description here

Not too bad, of course you can adjust the basis used as appropriate. At some point, if you have the information to create a nonlinear model that will be better of course.

$\endgroup$
3
  • $\begingroup$ If I understand correctly, his real data are not linear and the model probably is unknown. $\endgroup$ Apr 7, 2014 at 16:31
  • $\begingroup$ NonlinearModelFit might also be appropriate as well. In either case, I suspect they don't want to "smooth" as much as fit. The picture they show doesn't look to terrible to model. $\endgroup$
    – chuy
    Apr 7, 2014 at 16:34
  • $\begingroup$ Sorry, I'm facing a similar problem. How can I use NonlinearModelFit in this algorithm? $\endgroup$
    – Orion
    Jun 23, 2020 at 8:32
7
$\begingroup$

I was thinking - how can we average but without loss of points? Well we can randomly sample, interpolate and average - as many times as we want.

Let's take a look at more complicated data:

data = Flatten[Table[{i (1 + RandomReal[0.1]), j (1 + RandomReal[0.1]), 
       i^2 + j^2}, {i, 0, 100}, {j, 0, 100}], 1];

This is like ~ 10^4 points. Grab samples by 1000 - and many of those - and Interpolate - ListContourPlot anyways does that:

Do[f[k] = Interpolation[RandomSample[data, 1000], InterpolationOrder -> 1], {k, 100}]

Average them:

F[x_, y_] = Sum[f[k][x, y], {k, 100}]/100;

Now let's see:

Show[
  ContourPlot[F[x, y], {x, 0, 100}, {y, 0, 100}, 
   ContourShading -> None, ContourStyle -> Directive[Red, Thick]],
  ListContourPlot[data, ContourShading -> None, Contours -> 9]
  ] // Quiet

enter image description here

Anyway - something along these lines.

============ OLD VERSION ============

Well, what I will suggest is brutally simple. You want something better. But for the sake of thought entertainment... You get too much detailed interpolation because you got so many data points. Reduce them?

Manipulate[
 Show[
  ListPlot3D[data[[1 ;; -1 ;; n]], PlotStyle -> Opacity[.5], 
   MeshFunctions -> {#3 &}, Mesh -> 10, 
   MeshStyle -> Directive[Opacity[.5], Thick], 
   PerformanceGoal -> "Quality"],
  ListPointPlot3D[data[[1 ;; -1 ;; n]], 
   PlotStyle -> Directive[Opacity[.5], Red]]
  , ImageSize -> 500]
 , {{n, 1}, 1, 147, 1, Appearance -> "Labeled"}]

enter image description here

ListContourPlot[data[[1 ;; -1 ;; 139]], Contours -> {80, 120}, ContourShading -> None]

enter image description here

You can also use moving average, but it also removes points:

ListContourPlot[MovingAverage[data, 50], ContourShading -> None, Contours -> {80, 120}]
$\endgroup$
5
  • 2
    $\begingroup$ Hey Vitaly, you really are into the screencast thing (I like it) ;D $\endgroup$
    – Yves Klett
    Apr 7, 2014 at 16:22
  • $\begingroup$ @YvesKlett haha yes I love that ;-) $\endgroup$ Apr 7, 2014 at 16:26
  • $\begingroup$ @VitaliyKaurov : Thanks a lot for the answer! The data[[a ;; b ;; c] method indeed improves the figure. But as you said, some data are lost and in my real example the quality is still not good enough. The moving average way is lovely. However, my {{x1, y1, f1}, {x2, y2, f2}, ...} list is actually not ordered in x1 < x2 < x3... but instead in a random way. Thus MovingAverage does not directly apply for my real need... Also, the screencasting is amazing :) $\endgroup$
    – Yi Wang
    Apr 7, 2014 at 16:33
  • $\begingroup$ @YiWang Then RahulNarain moving average method via nearest should work. $\endgroup$ Apr 7, 2014 at 16:35
  • $\begingroup$ @VitaliyKaurov : Thanks a lot for the update! Yes, random sampling, re-sampling and averaging is a very interesting idea :) $\endgroup$
    – Yi Wang
    Apr 8, 2014 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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