26
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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.. $\endgroup$
    – Silvia
    Aug 8, 2012 at 18:06

4 Answers 4

14
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ theData = RandomReal[NormalDistribution[], {10000, 2}]; i = Image@ ListPlot[theData, PlotStyle -> {Black, Opacity[.3]}, Axes -> False] Lighter@Lighter@i $\endgroup$ Aug 7, 2012 at 21:18
  • 1
    $\begingroup$ 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 $\endgroup$
    – acl
    Aug 7, 2012 at 22:41
  • $\begingroup$ @verde Did you forget a semicolon after the list plot? Without it it just exceeds the recursion limit... $\endgroup$
    – Ajasja
    Aug 8, 2012 at 7:48
  • $\begingroup$ Nope, but there is a <return> after the first plot that is not showed in the comment $\endgroup$ Aug 8, 2012 at 14:14
19
$\begingroup$

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]]]
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Aug 8, 2012 at 9:05
  • $\begingroup$ @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 $\endgroup$ Aug 8, 2012 at 12:20
  • $\begingroup$ 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. $\endgroup$ Aug 8, 2012 at 12:34
  • $\begingroup$ @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 $\endgroup$ Aug 8, 2012 at 14:10
  • $\begingroup$ @OleksandrR. see last edit using Interpolating Function $\endgroup$ Aug 8, 2012 at 14:26
8
$\begingroup$

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

$\endgroup$
4
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – user484
    Jul 9, 2014 at 20:15

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.