# Plotting large datasets

I have a lot of data (e.g., 64x8192), to plot and it consumes too much memory (more than 2 GB) causing Mathematica to return "no memory available".

For example,

lst = Table[Sin[x^2 + y^2], {x, 2^7}, {y, 2^13}];
ListDensityPlot[lst]


Is there any way to cope with it?

-
Strongly related: mathematica.stackexchange.com/q/140/121 –  Mr.Wizard Jun 28 '12 at 7:02
While this related thread does not deal with downsampling and memory consumption, it allows you to display only a smaller area of a large plot and easily pan through/zoom in/out of the whole. –  István Zachar Jun 28 '12 at 8:05

Scroll to the end for the latest edit.

If you really want all the points in a large data set, I'd use a more basic display method:

Let's say the matrix of data is m, then do this (assuming the entries m[[i, j]] are real numbers):

Image[Rescale[m]]


This will make one pixel for each data point without trying to do any interpolation as is the default in density plots.

Edit

As suggested in the comment, here is a comparison of the timings:

(* A test matrix with random entries between 0 and 1: *)
m = RandomReal[1, {300, 300}];
Timing@Image[m]


Timing@ListDensityPlot[m, MaxPlotPoints -> 1000]


So the ListDensityPlot is several orders of magnitude slower than my Image suggestion. This conclusion seems to be almost independent of the choice for MaxPlotPoints.

With the Image approach, it is no problem to create extremely large images (say $5000\times 5000$). I would not advise using ListDensityPlot for anything like this.

If you want a colored output, just use Colorize as follows:

m = RandomReal[1, {300, 300}];

Timing[plot = Colorize[Image[m], ColorFunction -> "LakeColors"]]


Edit 2

The final question would be, how to re-create the frame and ticks that you get from ListDensityPlot and its relatives. Let's assume that the range of x and y values corresponding to the matrix dimensions of m are $x \in [0, 2]$ and $y\in[0, 3]$, respectively. Then you get the standard display for the plot we just created as follows:

xmin = 0;
ymin = 0;
xmax = 2;
ymax = 3;
Graphics[Inset[
Show[plot, AspectRatio -> Full],
{xmin, ymin},
{0, 0},
{xmax - xmin, ymax - ymin}],
PlotRange -> {{xmin, xmax}, {ymin, ymax}},
Frame -> True]


The trick here is to align the bottom left corner {0,0} of the Image to the lower coordinate limits of the PlotRange I chose for the Graphics, and then (using the fourth argument of Inset) to stretch the image to fit exactly into that PlotRange. The content is only stretchable this way if you show it with the option AspectRatio -> Full.

Edit 3

For completeness, here is the example matrix given in the question, constructed in a slightly modified way to make its entries lie between 0 and 1. Here, Rescale turns out to be very slow when applied to an array of exact numbers, so I speed things up by adding N[] during the construction of the matrix. Thanks to MichaelE2 for pointing this out in the comment. The matrix construction uses Outer and automatically threads Sin over lists in its argument.

In this construction, the x coordinate is represented by Range[2^7] and the y coordinate by Range[2^13] - if you switch these then the plot will have the axes switched. The Reverse below is needed to have the vertical axis pointing up:

s = Sin[Outer[Plus, N[Reverse[Range[2^13]]^2], N[Range[2^7]^2]] ];
plot = Colorize[Image[s], ColorFunction -> "LakeColors"];

xmin = 1;
ymin = 1;
xmax = 2^7;
ymax = 2^13;
Graphics[Inset[
Show[plot, AspectRatio -> Full], {xmin, ymin}, {0, 0}, {xmax - xmin,
ymax - ymin}], PlotRange -> {{xmin, xmax}, {ymin, ymax}},
Frame -> True, AspectRatio -> 1]


The point of this example is to show how much more detail we get with this approach, compared to the slower ListDensityPlot.

Alternative: ArrayPlot

I mentioned using Rescale earlier because the arguments to Image have to be in the unit interval to be mapped onto the grayscale range from white to black. This grayscale can later be colorized too. It's important to keep this fact in mind: matrix entries that fall outside the interval $[0, 1]$ will not lie in the range where the grayscale varies and will therefore appear at the same extreme grayscale (or color afte using Colorize).

Once we add the rescaling and the frame ticks, the Image approach de-facto converges to the built-in function ArrayPlot. This is another easy-to-use way of getting a fast plot of a large matrix. You would use it like this: ArrayPlot[m, ColorFunction -> "LakeColors", DataReversed -> True]. I just tried this with the last example matrix, and it works nicely.

One more observation

According to this linked answer, when using ColorFunction, ArrayPlot isn't as fast as the direct creation of a Raster image. The solution proposed there is not Image as I used here, but instead something like this: Graphics[Raster[Rescale[...], ColorFunction -> "LakeColors"]]

To try this with my frame construction, you would have to write it like this (where s is the same as defined in Edit 3):

xmin = 1;
ymin = 1;
xmax = 2^7;
ymax = 2^13;
Graphics[Inset[
Show[
Graphics[Raster[s, ColorFunction -> "LakeColors"]],
AspectRatio -> Full,
{xmin, ymin},
{0, 0},
{xmax - xmin, ymax - ymin}
],
PlotRange -> {{xmin, xmax}, {ymin, ymax}},
Frame -> True, AspectRatio -> 1
]


When I wrap this in Timing or AbsoluteTiming, I indeed get even lower values. However, I don't really believe those values because the actual time it takes for the plot to appear in the notebook is noticeably slower with this last approach. Nevertheless, I think it's worth including this here.

More timing results

To test these display methods with the Mandelbrot set of this linked answer, I did a different timing measurement based on recording the absolute time before and after the image/graphic has been displayed. Here is the code to prepare the pictures:

mandelComp =
Compile[{{c, _Complex}},
Module[{num = 1},
FixedPoint[(num++; #^2 + c) &, 0, 99,
SameTest -> (Re[#]^2 + Im[#]^2 >= 4 &)]; num],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];

data = ParallelTable[
a + I b, {a, -.715, -.61, .0001}, {b, -.5, -.4, .0001}];

rc = Rescale[mandelComp[data]];


Now the timing tests (not showing the images that are generated):

a1 = AbsoluteTime[]; Graphics[
Raster[rc, ColorFunction -> "TemperatureMap"], ImageSize -> 800,

a2 = AbsoluteTime[]

a2 - a1

(* ==> 6.878487 *)


This is pretty slow. Compare this to the Image approach:

a1 = AbsoluteTime[]; Colorize[Image[rc, ImageSize -> 800],
ColorFunction -> "TemperatureMap"]

a2 = AbsoluteTime[]

a2 - a1

(* ==> 1.119959 *)


This is much faster. The ArrayPlot method is twice as slow as the last result, but still faster than the Raster method:

a1 = AbsoluteTime[];
ArrayPlot[rc, ImageSize -> 800, ColorFunction -> "TemperatureMap"]

a2 = AbsoluteTime[]

a2 - a1

(* ==> 2.089126 *)


Therefore, my conclusion is: when in doubt, go with Image.

Edit (19 June 2014)

It seems useful to wrap the conclusion of this answer in a general-purpose function that works similarly to ListDensityPlot but is based on directly creating an Image from the data:

Clear[listPixelPlot] ;
Options[listPixelPlot] = {Frame -> True, FrameLabel -> {},
ColorFunction -> "LakeColors" ,
AspectRatio -> 1, DataRange -> All , ImageSize -> Automatic,
PlotRange -> All, PlotRangeClipping -> True};

listPixelPlot[array_, OptionsPattern[]] :=
Module [{field, plot, xmin, xmax, ymin, ymax,
dataRange = OptionValue[DataRange],
plotRange =
Replace[OptionValue[PlotRange],
Full | All | Automatic -> {Automatic, Automatic}], range},
field = Reverse@Transpose[array] ;
{{xmin, xmax}, {ymin, ymax}} = range = If[
MatrixQ[dataRange] ,
Reverse[dataRange],
Map [ {0 , #} &, Dimensions[array] ]
] ;
plot = Colorize[Image[Rescale[N[field] ] ] ,
ColorFunction -> OptionValue[ColorFunction]] ;
Graphics[
Inset[
Show [plot, AspectRatio -> Full] , {xmin, ymin} , {0,
0} , {xmax - xmin, ymax - ymin}],
PlotRange ->
Table[plotRange[[i]] /. Full | All | Automatic -> range[[i]], {i,
2}],
Sequence @@
Map[# ->  OptionValue[#] & , {Frame, FrameLabel, ImageSize,
AspectRatio, PlotRangeClipping}]]
]


This function takes some of the same options that ListDensityPlot or ArrayPlot accept (see Options[). With the Mandelbrot test data rc from above, you would call this function (and test its timing) as follows:

a1 = AbsoluteTime[];
listPixelPlot[rc, ColorFunction -> "TemperatureMap", ImageSize -> 800]

a2 = AbsoluteTime[];

a2 - a1

(* ==> 1.679611 *)


This is only slightly slower than the fastest result above. The reason is simply that in this function the argument array is always processed by Rescale before converting it to an Image. To make sure Rescale sees a packed array of machine numbers, the core of the function reads Image[Rescale[N[field] ] ]. A lot of the code deals merely with option handling. It can be extended to provide additional options.

-
You've solve my problem more than I can expect, and you go further than kguler. I really enjoy your advice although I have accepted kguler's answer. –  yulinlinyu Jun 28 '12 at 5:14
@yulinlinyu You can accept another answer at any time if you wish. –  Alexey Popkov Jun 28 '12 at 6:12
@Jens, Not only speed and efficiency but also the looks are uncomparably nicer with Image; you may consider posting an image to highlight this additional advantage using, say, a variation of the OP's example data m= Table[Sin[x^2 + y^2], {x, 2^10}, {y, 2^10}]; Image[m]//Colorize. (+++1 needless to say) –  kglr Jun 29 '12 at 0:52
I really enjoy your solution. But I've met some other problem. Since I cannot paste images in the comment, please compare dt = Table[Log10[Abs@Sin[0.2 x] Exp[-y^2]], {x, 100}, {y, 10}];ListDensityPlot[dt, PlotRange -> All, ColorFunction -> ColorData["GreenPinkTones"], ColorFunctionScaling -> True, MaxPlotPoints -> 1000, ImageMargins -> 5] and img = Colorize[Image[dt], ColorFunction -> ColorData["GreenPinkTones"]];Graphics[Inset[ Show[img, AspectRatio -> Full], {0, 0}, {0, 0}, {1, 1}], PlotRange -> {{0, 1}, {0, 1}}, Frame -> True]. Any advice? –  yulinlinyu Jun 29 '12 at 0:54
@yulinlinyu Yes - the data in td have to be rescaled first before they can be used in Image. The numbers always have to be between 0 and 1 for this to work. So you could in this case replace td by Rescale[td]. Or even better, prepare the data so they are rescaled to begin with, as I did in the last example I posted. Then the image will show the color range as expected. Since the x axis has only 10 data, the colors won't look smooth in Image, but the point is that for large dimensions you won't notice this discreteness. –  Jens Jun 29 '12 at 1:02

ListDensityPlot uses all the points in the dataset. You can set MaxPlotPoints option to control downsampling:

  ListDensityPlot[lst, MaxPlotPoints -> 1000]


Memory consumption using MaxMemoryUsed:

 miu0 = MaxMemoryUsed[]
lst = Table[Sin[x^2 + y^2], {x, 2^7}, {y, 2^13}];
usg0 = MaxMemoryUsed[] - miu0
miu1 = MaxMemoryUsed[]
ListDensityPlot[lst, MaxPlotPoints -> 1000];
usg1 = MaxMemoryUsed[] - miu1
miu2 = MaxMemoryUsed[]
ListDensityPlot[lst];
usg2 = MaxMemoryUsed[] - miu2


Incremental memory consumption:

 {usg0, usg1, usg2}
(* ==> {96 881 240, 447 831 576, 3 243 923 016} *)


That is, the input data and plotting with MaxPlotPoints->1000 should not trigger Out Of Memory issues, while without the explicit MaxPlotPoints producing the plot consumes more than 3G of additional memory beyond what is needed for holding the input data.

-
I think this is just I want! –  yulinlinyu Jun 28 '12 at 3:13
@yulinlinyu, thank you for the accept. –  kglr Jun 28 '12 at 3:45
@yulinlinyu, you can change your mind about which answer to accept any time; and if this were my question i would choose Jens's answer. Try, for example, lst = Table[Sin[x^2 + y^2], {x, 2^10}, {y, 2^10}]; Image[lst]//Colorize to see why. –  kglr Jun 28 '12 at 10:33
Thank you for your kindness. The code in your advice is so amazing that I cannot imagine it before. –  yulinlinyu Jun 28 '12 at 11:56

You can use a Span to only use some of the points. E.g.

lst[[;; ;; 10, ;; ;; 10]]


will take every 10th point (in each dimension). This is of course a rather crude way of down sampling.

-
This is really a easiest way! –  yulinlinyu Jun 28 '12 at 12:01

The correct way to deal with this issue is to reduce the number of points in your plot to a manageable level.

Computer monitors have only so many pixels and human eyes have only so much resolution. Having this many points in a graphic will not convey more information or make it look clearer.

The example function you gave, Sin[x^2+y^2] cannot be reasonably made into a graphic over the range you provided. The structure of the function wouldn't be noticeable.

-
I use Sin[x^2+y^2] here is just to generate a large list. Actually, the list comes from very complicated numerical calculations. –  yulinlinyu Jun 28 '12 at 2:08
If your function is smooth enough to make a sensible plot with an enormous amount of data points, it will be smooth enough to make a good looking plot with a much smaller number of points. Try re-sampling the data set so that it is manageable. –  Searke Jun 28 '12 at 2:12
Not quite true. Resolution of 1600x1200 gives of order 10^6 points. If your are to make a plot with about 5x10^5 points having 1 point per pixel, it will cause you plotting problems with information being not oversampled. –  Alexey Bobrick May 12 '14 at 11:25