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.

The problem

In a dense data set plotted with opaque points, data points can pile on top of each other. This can cause serious problems with the interpretation of the plots. I'll use as an example a scatterplot of normally distributed points.

theData = RandomReal[NormalDistribution[], {100000, 2}];

Unless one is careful, these plots end up showing a solid disk around the origin, which obscures the distribution since point densities above a certain level appear clipped. This also gives inappropriate emphasis to outliers, by plotting them with the same intensity as what might be a pile of a million points close together at the origin.

Here is the same plot shown at different sizes.

Row[ListPlot[
 theData,
 AspectRatio -> Automatic, ImageSize -> #, PlotStyle -> Directive[Black]
] & /@ {90, 180, 360}]

enter image description here

An artificial black disk appears at the origin and its apparent size depends on the plot size, even though the data are identical.

Here is the same plot with different settings for PointSize.

Row[ListPlot[
 theData,
 AspectRatio -> Automatic, ImageSize -> Medium, PlotStyle -> Directive[Black, PointSize[#]]
] & /@ {0.01, 0.005, 0.0025}]

enter image description here

The disk appears a different size depending on point size, even though the data are identical.

Here are plots of successively smaller subsets of the data.

Row[ListPlot[
 theData[[;; #]],
 AspectRatio -> Automatic, ImageSize -> Medium, PlotStyle -> Directive[Black]
] & /@ {100000, 50000, 25000}]

enter image description here

The distributions are obviously all the same, but they look different in the plot because of the pileup of points.

Solving pileup by changing opacity

For a given point size, one solution is to choose the minimum opacity where the darkest pixel in the plot is black. Here is an ugly, slow, procedural, brute-force way of choosing the correct opacity:

Block[
 {darkestPixel},
 darkestPixel=1;
 idealOpacity=0;
 While[
  darkestPixel>0,
  idealOpacity=idealOpacity+0.001;
  darkestPixel=Min@Flatten@ImageData@ColorConvert[
   Rasterize@ListPlot[
    theData,
    AspectRatio->Automatic,ImageSize->Medium,
    PlotStyle->Directive[Black,PointSize[0.01],Opacity[idealOpacity]],Axes->False
   ],
  "Grayscale"
  ]
 ]
]

Here are the resulting plots for different subsets of the data:

Row[ListPlot[
 theData[[;; #[[1]]]],
 AspectRatio -> Automatic, ImageSize -> Medium, 
 PlotStyle -> Directive[Black, PointSize[0.01], Opacity[#[[2]]]]
] & /@ {{100000, 0.018}, {50000, 0.042}, {25000, 0.1}}]

enter image description here

Solving pileup by changing point size

Another way is to choose a fixed opacity and increase the size of the points until the darkest pixel is black.

Block[
 {darkestPixel},
 darkestPixel = 1;
 idealPointSize = 0;
 While[
  darkestPixel > 0,
  idealPointSize = idealPointSize + 0.001;
  darkestPixel = Min@Flatten@ImageData@ColorConvert[
   Rasterize@
    ListPlot[theData,
     AspectRatio->Automatic, ImageSize->Medium, 
     PlotStyle->Directive[Black,PointSize[idealPointSize],Opacity[0.05]],Axes->False
    ],
   "Grayscale"
  ]
 ]
]

The results for different subsets of the data:

Row[
 ListPlot[
  theData[[;;#[[1]]]],
  AspectRatio->Automatic,ImageSize->Medium,
  PlotStyle->Directive[Black,PointSize[#[[2]]],Opacity[0.05]]
 ]&/@{{100000,0.003},{50000,0.01},{25000,0.015}}
]

enter image description here

The question

What's the best way to quickly and automatically choose the correct opacity, point size or combination of the two so that there is no "clipping" in a dense plot? Is it possible to do this inside of PlotStyle instead of as a pre-processing step? For bonus points, explain the best way to do this for multiple data sets on the same plot.

share|improve this question
    
Adjusting Opacity/PointSize to fit a plot with all data is nice, but there might be limitations. I myself encountered one when dealing with a similar data visualization problem, where the data is so dense, I used a very small Opacity and end with a discontinuity of the color.. –  Silvia Aug 8 '12 at 18:06

4 Answers 4

up vote 12 down vote accepted

At the risk of stating the obvious (and also not directly answering the question):

A ListPlot with opacity less than 1 is very very similar to a histogram, so why not just use that?

theData = RandomReal[NormalDistribution[], {10000, 2}];
opts = {ColorFunction -> Function[c, GrayLevel[1 - c]], 
   PlotRange -> {{-4, 4}, {-4, 4}}, ImageSize -> Medium};
Row@{
  DensityHistogram[theData, opts],
  SmoothDensityHistogram[theData, opts]
}

Examples of histograms

(When doing histogram of large data sets in mma 8 please be aware, that there is a memory leak in HistogramList. See this answer)

EDIT So to also answer the original question. From the histogram it should be possible to estimate the opacity and point size (although I would not recommend this as the rendering is much slower then a normal histogram). Here is my take on this

{bins, counts} = HistogramList[theData];
maxcounts = Max@counts;
binarea = Times @@ Part[Differences /@ bins, All, 1];
r = N@Sqrt[binarea/Pi];
s = 10;
ListPlot[theData, 
   PlotStyle -> 
    Directive[Black, Opacity[1/maxcounts*s], PointSize[r/s]], 
   ImageSize -> #] & /@ {90, 180, 360}

10 000 points: Listplots of different sizes 10 000 points

100 000 points: Listplots of different sizes 100 000 points

share|improve this answer
    
theData = RandomReal[NormalDistribution[], {10000, 2}]; i = Image@ ListPlot[theData, PlotStyle -> {Black, Opacity[.3]}, Axes -> False] Lighter@Lighter@i –  belisarius Aug 7 '12 at 21:18
1  
I was in the middle of writing my own routine to count number of points per "sqaure" (ie to coarse-grain, or bin, the data) when I realized mma must have this built-in. then I saw this! +1 –  acl Aug 7 '12 at 22:41
    
@verde Did you forget a semicolon after the list plot? Without it it just exceeds the recursion limit... –  Ajasja Aug 8 '12 at 7:48
    
Nope, but there is a <return> after the first plot that is not showed in the comment –  belisarius Aug 8 '12 at 14:14

tl;dr

Final results first:

(*Function Definition*)
ClearAll[opaFun];
Options[opaFun] = Options[ListPlot];
opaFun[points_, opts : OptionsPattern[]] := Module[{f, steps = 10 },
   f[x_] := Min[Norm /@ Flatten[ImageData@
              ListPlot[points, opts, Axes -> False, PlotStyle -> {Black, Opacity[x]}],1]]/Sqrt@3;
   Return@NestWhileList[{#, f[#]} &@(#[[1]] + 1/steps) &, {0, 1}, 
                         UnsameQ @@ ({##}[[All, 2]]) &, 2, steps][[-2, 1]]
];
(*Usage*)
opts = {AspectRatio -> 1/2, ImageSize -> 300, Axes -> False};
Timing@ListPlot[#, opts, PlotStyle -> {Black, Opacity[opaFun[#]]}] &@
                  RandomReal[NormalDistribution[], {10000, 2}]

Mathematica graphics

Now the full answer:

A little but time consuming experiment first. You'll not need to do it, as it is only for finding a model.

Let's see how the Min value of the plot intensity varies with Opacity:

theData = RandomReal[NormalDistribution[], {2000, 2}];
f[x_] := f[x] = 
   Min[Norm /@ Flatten[ImageData@Rasterize[
        ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
                   PlotStyle -> {Black, Opacity[x]}, Axes -> False]], 1]];
Plot[f[x] , {x, 0, .4}, PlotRange -> Full]

Mathematica graphics

So, it is an exponential. (Note: in the last edit to this post I got rid of the exponential model by fitting a few points with an Interpolation, and it works pretty nice)

Let's fit it:

data = Table[{i, f[i]}, {i, 0, 1, .1}]

model = a Exp[b x];
fit = FindFit[data, model, {a, b}, x];
modelf = Function[{t}, Evaluate[model /. fit]]

Show[ListPlot@data, Plot[modelf[x], {x, 0, 1}]]

Mathematica graphics

Now you are ready to set the min value of the brightness of the plot to whatever you want:

(The Sqrt@3 is a normalization factor for the intensity of the {1,1,1} RGB pixel.)

Let's use it:

opac = x /. Solve[# == a E^(b x)/Sqrt@3, x] /. fit & /@ {1/2, 1/4, 1/20, 1/200}

ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200, 
   PlotStyle -> {Black, Opacity[#[[1]]]}, Axes -> False] & /@ opac

Mathematica graphics

Edit

Let's pack that all together in a function and plot two very different point sets with the same maximum darkness.

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] :=
 Module[{f, a, b, model, fit, modelf, x},
  f[x_] := 
   f[x] = Min[
     Norm /@ Flatten[
       ImageData@
        Rasterize[
         ListPlot[points, Axes -> False, 
          PlotStyle -> {Black, Opacity[x]}]], 1]];

  model = a Exp[b x];
  fit = FindFit[Table[{i, f[i]}, {i, 0, 1, .1}], model, {a, b}, x];
  modelf = Function[{t}, Evaluate[model /. fit]];
  Return[x /. Quiet@Solve[# == modelf[x]/Sqrt@3, x][[1]] /. fit &@
    desiredOpacity];
  ]

theData  = RandomReal[NormalDistribution[], {2000, 2}];
theData1 = RandomReal[NormalDistribution[], {10000, 2}];

opad = opa[.5,  theData,  AspectRatio -> Automatic, ImageSize -> 200];
opad1 = opa[.5, theData1, AspectRatio -> Automatic, ImageSize -> 200];

Grid[{{ListPlot[theData, Axes -> False, PlotStyle -> {Black, Opacity[opad]},
                         AspectRatio -> Automatic, ImageSize -> 200], 
       ListPlot[theData1,Axes -> False, PlotStyle -> {Black, Opacity[opad1]}, 
                         AspectRatio -> Automatic, ImageSize -> 200]}}, 
   Frame -> All]

Mathematica graphics

The same, but darker

Mathematica graphics

Remember that the default plot is:

Mathematica graphics

Edit

Answering @Oleksandr comments, the following does not assume an exponential model:

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] := 
 Module[{f, model, x}, 
  f[x_] := f[x] = 
    Min[Norm /@ Flatten[ImageData@ Rasterize[
         ListPlot[points, Axes -> False, PlotStyle -> {Black, Opacity[x]}]], 1]];
  model = Interpolation[Table[{i, f[i]}, {i, 0, 1, .1}]];
  x /. FindRoot[desiredOpacity == model[x]/Sqrt@3, {x, 0, 1}][[1]]]
share|improve this answer
    
Don't the model parameters depend on the distribution of the data? It seems to me that, as you'd have to fit the model for every distribution that you want to plot, the time-intensive step can't be factored out in general. –  Oleksandr R. Aug 8 '12 at 9:05
    
@OleksandrR. The time intensive step was an experiment. For each distribution you need an appropriate model, yes but you may "guess" it from an interpolation dome over data = Table[{i, f[i]}, {i, 0, 1, .1}] instead of fitting the data like I did here. Posting that generalization later –  belisarius Aug 8 '12 at 12:20
    
I think I've now convinced myself that it'll always be exponential, but the parameters will still need to be determined each time. You could also use your current Plot approach with PlotPoints/MaxRecursion rather than a fixed step. –  Oleksandr R. Aug 8 '12 at 12:34
    
@OleksandrR. Getting the parms for a specific distribution is fast. Unless you have to do that a thousand times, I guess it is not a problem –  belisarius Aug 8 '12 at 14:10
    
@OleksandrR. see last edit using Interpolating Function –  belisarius Aug 8 '12 at 14:26

Maybe the clipping can be avoided by not using ListPlot. Here is my variation using Point that you might find helpful (or maybe not), no time for extensive testing of all your conditions:

dist = Rescale[1 - EuclideanDistance[{0, 0}, #] & /@ theData];

Graphics[{Thick, 
  Point[theData, 
  VertexColors -> (Blend[{{0, {Black, Opacity[.01]}}, {1, {Black, 
  Opacity[.1]}}}, #1] & /@ dist)]}, AspectRatio -> 1, 
  Axes -> True, PlotRange -> {{-4, 4}, {-4, 4}}]

enter image description here

share|improve this answer

In version 10, ListPlot chooses point size automatically. They have used this exact example to demonstrate it here.

share|improve this answer
    
That's pretty cool! I wonder what happens if you have two different datasets in the same plot -- are their point sizes adapted independently or are they the same? I can see arguments in favour of either. –  Rahul Jul 9 at 20:15

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.