# ListDensityPlot behaviour and gnuplot's splot analogue

I would like to reproduce the following plot in Mathematica: This was produced with Gnuplot's splot from a data file with the structure (a text file, each line containing):

x y f(x,y)

The function I found in Mathematica to do this job was ListDensityPlot but there are a few problems:

1. The data file is rather large (~40 MB): 947376 lines and this is too large for Mathematica to handle (it takes Gnuplot about 3 seconds to produce the image above but Mathematica cannot handle this). Note, I do not want a solution that selects only few lines from the result, actually, I would like to double the size of the list containing the data to show the periodicity.
2. Even when I restrict myself to a very small data slice, I get strange behaviour, such as:

The peak on the y-axis is not supposed to be there (or more precisly, there are some data point missing for which $y=1.5$), the y data axis is always $(0,1.5)$ with step $0.1$ and thus there are 16 data points for each x axis point, therefore data[[1:96]] should display full range of y axis for each x point.

But things can get even more out of hand:

As you can see, the last data point is correctly the maximum y point but also, the x point is 2.96546*10^7 but on the plot, this is just an empty region.

I am probably doing something very wrong (this behaviour may be related to the interpolation function of ListDensityPlot but I am not sure about this because InterpolationOrder->1 does not change a thing) but I have not been able to figure out what. I can live with gnuplot for now but I would really like to figure out how to reproduce the same plots in Mathematica --- first to get rid of the problem 2 and then to make Mathematica handle large amount of data with the same swiftness as gnuplot.

For those interested, the data file is located here. The 7z archive contains data.dat which is used in Mathematica, data.dat2 which is the same data reformatted for use with gnuplot and plot1.plt to generate the gnuplot plot and finally the plot shown above.

Mathematica seems to handle arrays much better than lists of {{x,y,f[x,y]},..} tuples for functions like ListContourPlot and ListDensityPlot. So, when making the data files, if you can, just keep it stored as an array of z-values.

I ran into the same trouble as the OP when trying to plot the data as-is,

temp = Import["data.dat"];
ListDensityPlot[temp]


and it just hung for more minutes than I was willing to wait. But since the data is on a rectangular grid, just reshape it as an array and drop the x and y values

Dimensions@temp
data = GatherBy[temp, {#[[2]] &, #[[1]] &}];
data = MapAt[Last@*First, data, {All, All}];
Dimensions@data
{xrange, yrange, zrange} = MinMax /@ Thread[temp]
(* {947376, 3} *)
(* {16, 59211} *)
(* {{2.96051*10^7, 5.92101*10^7}, {0., 1.5}, {274.887, 325.645}} *)


Even with the reshaping, this plot command takes about a minute

ListDensityPlot[data, DataRange -> {xrange, yrange}]


You can speed things up further by resampling the data to use fewer than 59k gridpoints in the horizontal direction

data2 = ArrayResample[data, {16, 500}];

ListDensityPlot[data2, DataRange -> {xrange, yrange}]


That last plot is downright speedy.

I think Gnuplot's splot is closer to Mathematica's ArrayPlot:

data = ReadList["data.dat", {Number, Number, Number}]; // AbsoluteTiming
(* {1.325872, Null} *)

{xarray, yarray, zarray} = DeveloperToPackedArray /@ Transpose@Sort@data;

zmatrix = Transpose@Partition[zarray, yarray // Union // Length];

{xmin, xmax} = xarray[[{1, -1}]]
(* {2.96051*10^7, 5.92101*10^7} *)
{ymin, ymax} = yarray[[{1, -1}]]
(* {0., 1.5} *)
{zmin, zmax} = {Min@#, Max@#} &@zarray;

AbsoluteTime[];
Legended[DensityPlot[False, {x, xmin, xmax}, {y, ymin, ymax},
AspectRatio -> 1/GoldenRatio, PlotRangePadding -> None]~Show~
ArrayPlot[zmatrix, DataRange -> {{xmin, xmax}, {ymin, ymax}}, PlotRange -> All,
DataReversed -> True, ColorFunction -> "SunsetColors"],
BarLegend[{"SunsetColors", {zmin, zmax}}]]
AbsoluteTime[] - %%
(*1.8237344*)
`