Search code examples
3dgnuplotdensity-plot

How to plot (x,y,z) points showing their density


My data file is a set of (x,y,z) points located around the axis origin. They represent points where a kind of measure has failed. These points are in this link.

Gnuplot can plot them,

set encoding iso_8859_1
set term postscript eps enhanced color size 4.7in,4in
set xlabel "X"
set ylabel "Y"
set zlabel "Z"
set output "test_gp.eps"
set style line 1 lc rgb '#0060ad' pt 7 ps 0.5 lt 1 lw 0.5 # --- blue
set style fill  transparent solid 0.15 noborder
splot "data.dat" u 1:2:3 w p ls 1 title "P_{error}"

with the resulting plot

enter image description here

The problem is that the figure does not show where the error is more likely to occur. I would like to show the variations in the density of points if it is possible.

I do not know if it is possible to vary the color of the points or give them transparency and try to represent the locations with the most density of points.

If a representation of the density of 3D points is not possible, another way could be to make a projection density of the points in the three planes (x=0,y,z), (x,y=0,z), (x,y,z=0).

Regards

EDIT:

I can plot with different colors the success (0) and error (1) point locations of my experiment. The file in this link has a 4th column with all data samples ( 0's and 1's).

The splot figure

splot "data_all.dat" u 1:2:3:4 w points ls 1 palette title "P_{error}"

is

enter image description here

but this figure is not showing the density of points. For example, in Mathematica the density plot of these data samples is

enter image description here

How could I get the density plot with Gnuplot?. It is likely that Mathematica is doing an interpolation of the points in the middle and give them values between 0 and 1, but I do not know how to achieve it with Gnuplot.


Solution

  • @user1993416, I guess you can do something with gnuplot. You might want to play with the parameter Radius for determining the number of points around a certain point within this radius and calculating the density. With my 8 year old computer 1000 points need approx. 2-3 minutes.

    Update: After revisiting this old post, I noticed that there was a unnecessary loop in test data creation and the script can be made much faster by using only arrays, i.e. no writing to tables which seems to be much slower than arrays. Now, it takes about approx. 2, 7 and 26 seconds for 500, 1000 and 2000 points. So, runtime is O(n^2).

    Script: (works with gnuplot>=5.2.0)

    ### 3D density plot
    reset session
    set term wxt
    TimeStart = time(0.0)
    
    # create some random test data
    set table $Data
        set samples 1000
        plot '+' u (invnorm(rand(0))):(invnorm(rand(0))):(invnorm(rand(0))) with table
    unset table
    
    # put the datafile/dataset into arrays
    stats $Data u 0 nooutput
    N = STATS_records
    array X[N]
    array Y[N]
    array Z[N]
    array C[N]
    stats $Data u (X[$0+1]=$1, Y[$0+1]=$2, Z[$0+1]=$3) nooutput
    
    # look at each datapoint and its sourrounding
    Radius = 2
    V      = pi*4/3*Radius**3   # sphere volume
    p0     = 0;  dp = 10        # progress steps
    do for [i=1:N] {
        if (p=real(i)/N*100, p>=p0 ? p0=(floor(p/dp)+1)*dp : 0) {
            print sprintf("Progress: %.0f%%", p)
        } 
        C[i] = 0
        stats $Data u ( C[i] = C[i] + \
              (sqrt((X[i]-$1)**2 + (Y[i]-$2)**2 + (Z[i]-$3)**2) <= Radius )) nooutput
    }
    print sprintf("Time elapsed: %.3f sec",time(0.0)-TimeStart)
    
    set key noautotitle
    set xyplane relative 0
    set view equal xyz
    set view 65,45,1.2
    set palette rgb 33,13,10
    
    splot X u (X[$0+1]):(Y[$0+1]):(Z[$0+1]):(C[$0+1]/V) w p ps 1 pt 7 lc palette z
    
    set terminal gif animate delay 20
    set output "SO53659762.gif"
        a0     = 40
        a1     = 60
        Frames = 24
        do for [i=0:Frames] {
            set view 65,(sin(2*pi*i/(Frames+1))*(a1-a0)+a0)
            replot
        }
    set output
    ### end of script
    

    Result: (output to gif terminal)

    enter image description here