Gnuplotting

Create scientific plots using gnuplot

September 29th, 2014 | No Comments

Some time ago I introduced already a waterfall plot, which I named a pseudo-3D-plot. In the meantime, I have been asked several times for a colored version of such a plot. In this post we will revisit the waterfall plot and add some color to it.

Colored waterfall plot

Fig. 1 Waterfall plot of head related impulse responses. (code to produce this figure, color palette, data)

In Fig. 1 the same head related impulse responses we animated already are displayed in a slightly different way. They describe the transmission of sound from a source to a receiver placed in the ear canal dependent on the position of the source. Here, we show the responses for all incident angles of the sound at once. At 0° the source was placed at the same side of the head as the receiver.

The color is added by applying the Moreland color palette, which we discussed earlier. The palette is defined in an extra file and loaded, this enables easy reuse of defined palettes. In the plotting command the palette is enabled with the lc palette command, that tells gnuplot to use the palette as line color depending on the value of the third column, which is given by color(angle).

load 'moreland.pal'
set style fill solid 0.0 border
limit360(x) = x<1?x+360:x
color(x) = x>180?360-x:x
amplitude_scaling = 200
plot for [angle=360:0:-2.5] 'head_related_impulse_responses.txt' \
    u 1:(amplitude_scaling*column(limit360(angle)+1)+angle):(color(angle)) \
    w filledcu y1=-360 lc palette lw 0.5

To achieve the waterfall plot, we start with the largest angle of 360° and loop through all angles until we reach 0°. The column command gives us the corresponding column the data is stored in the data file, amplitude_scaling modifies the amplitude of the single responses, and +angle shifts the data of the single responses along the y-axis to achieve the waterfall.

Even though the changing color in the waterfall plot looks nice you should always think if it really adds some additional information to the plot. If not, a single color should be used. In the following the same plot is repeated, but only with black lines and different angle resolutions which also have a big influence on the final appearance of the plot.

Colored waterfall plot

Fig. 2 Waterfall plot of head related impulse responses with a resolution of 5°. (code to produce this figure, data)

Colored waterfall plot

Fig. 3 Waterfall plot of head related impulse responses with a resolution of 2.5°. (code to produce this figure, data)

Colored waterfall plot

Fig. 4 Waterfall plot of head related impulse responses with a resolution of 1°. (code to produce this figure, data)

June 23rd, 2014 | 7 Comments

Occasionally it is a good idea to create a zoom of some part of your main plot, especially if you have a small part of your plot where the data is hiding each other.

Including a zoom with multiplot

Fig. 1 Including a zoom into your figure to emphasize some data. (code to produce this figure, data)

In Fig. 1 the interaural time difference between a sound signal reaching the two ears of a listener is plotted with different colors for different frequencies. The data is very dense around 0°, so we include a zoom into this region in the same figure at a free place.

This can be done via multiplot and the plotting of the same data in a smaller figure.

set origin 0.12,0.17
set size 0.45,0.4
set xrange [-10:0]
set yrange [0:0.1]
plot for [n=2:13] 'itd.txt' u 1:(column(n)*1000) w lines ls n

The tricky part is that we have a grid in our main figure and if we do nothing the grid will also be visible in the zoomed in version as shown in Fig. 2.

Including a zoom with multiplot without grid correction

Fig. 2 Including a zoom into your figure, without correcting the grid. (code to produce this figure, data)

To avoid this we have to hide the grid in the background of the zoomed graph. This is done with the trick of placing an empty white rectangle at the place the zoom plot should appear in the figure.

set object 1 rect from -88,0.03 to -49,0.41
set object 1 rect fc rgb 'white' fillstyle solid 0.0 noborder

This will then finally lead to the desired result presented in Fig. 1.

June 3rd, 2014 | No Comments

Suppose we have a file containing the following data from the plotting data tutorial:

# absolute_data.txt
# X   Y
1   2
2   3
3   2
4   1
Plotting absolute data

Fig. 1 Plotting absolute data points. (code to produce this figure, data)

This can be plotted in a straightforward manner and will result in Fig. 1. Now suppose we have the same data points stored as relative coordinates in our data file, resulting in:

# relative_data.txt
# deltaX deltaY
1   2
1   1
1   -1
1   -1

If we want to plot that data in gnuplot we have to keep track of the current position manually by storing its (x,y) value as variables by

x=0.; y=0.
plot 'relative_data.txt' u (x=x+$1):(y=y+$2) w p ls 1

Here, we define the starting point to be (0,0) and add to it the values from the first and second column for every line of the data file. By doing so, this results again in Fig. 1. Note, that the addition is always performed first, before the resulting point is plotted which means we get no point at (0,0). Now assume that we also want to add steps going from point to point as shown in Fig. 2. Gnuplot has the steps plotting style to achieve this, but we have to be carefully regarding our (x,y) variables.

Plotting relative data

Fig. 2 Plotting relative data points. (code to produce this figure, data)

Every single line of a plotting command is executed after each other which means our (x,y) pair will not be set to (0,0), but to (4,1) at the beginning of the second line of the plotting command. To avoid this we introduce another (a,b) pair for the second line and get finally.

x=0.; y=0.
a=0.; b=0.
plot 'relative_data.txt' u (x=x+$1):(y=y+$2) w steps ls 2,\
     ''                  u (a=a+$1):(b=b+$2) w points ls 1

April 16th, 2014 | 4 Comments

Gnuplot comes with the possibility of plotting histograms, but this requires that the data in the individual bins was already calculated. Here, we start with an one dimensional set of data that we want to count and plot as an histogram, similar to the hist() function we find in Octave.

Histogram of angle data

Fig. 1 Two different distributions of measured angles. (code to produce this figure, hist.fct, data)

In Fig. 1 you see two different distributions of measured angles. They were both given as one dimensional data and plotted with a defined macro that is doing the histogram calculation. The macro is defined in an additional file hist.fct and loaded before the plotting command.

binwidth = 4
binstart = -100
load 'hist.fct'
plot 'histogram.txt' i 0 @hist ls 1,\
     ''              i 1 @hist ls 2

The content of hist.fct, including the definition of @hist looks like this

# set width of single bins in histogram
set boxwidth 0.9*binwidth
# set fill style of bins
set style fill solid 0.5
# define macro for plotting the histogram
hist = 'u (binwidth*(floor(($1-binstart)/binwidth)+0.5)+binstart):(1.0) smooth freq w boxes'

For a detailed discussion on why @hist calculates a histogram you should have a look at this discussion and the documentation about the smooth freq which basically counts points with the same x-value. The other settings in the file define the width of a single bin plotted as a box and its fill style.

Histogram of angle data

Fig. 2 Two different distributions of measured angles. The bins of the histograms are shifted to be centered around 0°. (code to produce this figure, hist.fct, data)

It is important that the two values binwidth and binstart are defined before loading the hist.fct file. These define the width of the single bins and at what position the left border of a single bin should be positioned. For example, let us assume that we want to have the bins centered around 0° as shown in Fig. 2. This can be achieved by settings the binstart to half the binwidth:

binwidth = 4
binstart = -2
load 'hist.fct'
plot 'histogram.txt' i 0 @hist ls 1,\
     ''              i 1 @hist ls 2

February 21st, 2014 | 4 Comments

For the measurement of distances T-shaped arrows are often used to highlight the length. In gnuplot there is an easy way to achieve this.

Diffraction on a slit

Fig. 1 Diffraction of light for a slit with the length d. (code to produce this figure)

Have a look at Fig. 1 which tries to explain the diffraction phenomenon of a slit with the length d. At a distance a the diffraction pattern is drawn. The different distances, the distance between the first minima of the diffraction pattern, and the wave length are indicated by T-shaped arrows. This kind of arrays can be achieved in gnuplot with the following code.

Theads = 'heads size 0.5,90 front ls 201'
set arrow from -24,-2 to -24, 2       @Theads
set arrow from -22, 2 to -21.44,1.92  @Theads
set arrow from 1.5,-pi to 1.5,pi      @Theads
set arrow from -22,2.5*pi to 0,2.5*pi @Theads

Here, 90 is the relevant entry after size as it describes the opening angle of the arrow head.
In addition, an open circle is drawn to indicate the angle θ. This is achieved by specifying the opening angle for the circle object.

set object circle at -22,0 size 6 arc [-8:0]

January 31st, 2014 | 2 Comments

And another plot of the world. This time we are dealing with the raster data from Natural Earth. This data is normally available as tif-files. To use them in gnuplot we have to convert them first, then we can create a plot as shown in Fig. 1.

Cross-blended Hypsometric Tints

Fig. 1 Cross-blended hypsometric tint plot of the world. (code to produce this figure, data)

The conversion is done by the convert_natural_earth script. There the tif-file is first scaled down to the desired resolution using imagemagick. Afterwards it is converted to a text file and reordered for the splot command of gnuplot. The text file includes the longitude, latitude and three rgb color values.
You have to invoke the script in the following way.

$ ./convert_natural_earth $RES $FILE

where $RES is the desired resolution in pixel of your gnuplot plot and $FILE the input tif-file.
After finished we can plot the resulting text file simply by

set datafile separator ','
set size ratio -1
plot 'HYP_50M_SR_W_350px.txt' w rgbimage

which results in Fig. 1.

Natural Earth II Shaded Relief and Water

Fig. 2 Natural Earth II shaded relief and water plot of the world. (code to produce this figure, data)

The image can also be projected on a 3D figure of the world as shown in Fig. 2. To achieve this the three rgb values have to be summarized in one value and the rgb variable line color option has to be chosen together with pm3d.

rgb(r,g,b) = 65536 * int(r) + 256 * int(g) + int(b)
set mapping spherical
set angles degrees
splot 'NE2_50M_SR_W_700px.txt' u 1:2:(1):(rgb($3,$4,$5)) w pm3d lc rgb variable

December 21st, 2013 | 2 Comments

After plotting the world several times we will concentrate on a smaller level this time. Ben Johnson was so kind to convert the part dealing with the USA of the 10m states and provinces data set from natural earth to something useful for gnuplot. The result is stored in the file usa.txt.

USA election

Fig. 1 Election results of single U.S. states. (code to produce this figure, USA data, election data)

Two double lines divide the single states. This allows us to plot a single state with the help of the index command. At the end of this post the corresponding index numbers for every state are listed.
In addition to the state border data we have another file that includes results from an example election and strings with the names of the states. The election result can be 1 or 2 – corresponding to blue and red. With the help of these two data sets we are able to create Fig. 1 and Fig. 2.
For drawing a single state in red or blue we first collect the results for every single state in the string variable ELEC. The stats command is suitable for this, because it parses all the data but doesn’t try to plot any of them. During the parsing of every line the election result stored in the second column will be added at the end of the ELEC variable.

ELEC=''
stats 'election.txt' u 1:(ELEC = ELEC.sprintf('%i',$2))

In a second step we plot the state borders and color the states with the help of the ELECstring. ELEC[1:1] will return the election result for the state with the index 0.

plot for [idx=0:48] 'usa.txt' i idx u 2:1 w filledcurves ls ELEC[idx+1:idx+1],\
                    ''              u 2:1 w l ls 3

Alaska and Hawaii are then added with additional plot commands and the help of multiplot.

The data file with the election results includes also the names of the single states and a coordinates to place them. This allows us to put them in the map as well, as you can see in Fig. 2.

USA election

Fig. 2 Names and election results of single U.S. states. (code to produce this figure, USA data, election data)

The plotting of the state names is easily achieved by the labels plotting style:

plot for [idx=0:48] 'usa.txt' i idx u 2:1 w filledcurves ls ELEC[idx+1:idx+1],\
                    ''              u 2:1 w l ls 3,\
                    'election.txt'  u 6:5:3 w labels tc ls 3

At the end we provide the list with the index numbers and the corresponding states. If you want to plot a subset of states – as in Fig. 2 – you should adjust the xrange and yrange values accordingly.

0  Massachusetts
1  Minnesota
2  Montana
3  North Dakota
4  Idaho
5  Washington
6  Arizona
7  California
8  Colorado
9  Nevada
10 New Mexico
11 Oregon
12 Utah
13 Wyoming
14 Arkansas
15 Iowa
16 Kansas
17 Missouri
18 Nebraska
19 Oklahoma
20 South Dakota
21 Louisiana
22 Texas
23 Connecticut
24 New Hampshire
25 Rhode Island
26 Vermont
27 Alabama
28 Florida
29 Georgia
30 Mississippi
31 South Carolina
32 Illinois
33 Indiana
34 Kentucky
35 North Carolina
36 Ohio
37 Tennessee
38 Virginia
39 Wisconsin
40 West Virginia
41 Delaware
42 District of Columbia
43 Maryland
44 New Jersey
45 New York
46 Pennsylvania
47 Maine
48 Michigan
49 Hawaii
50 Alaska

October 24th, 2013 | 2 Comments

If you want to put labels into a graph using the epslatex terminal you are probably interested in using a smaller font for these labes than for the rest of the figure. An example is presented in Fig. 1.

Photo density flux

Fig. 1 Photon flux density for different characteristic tail state energies E0 dependent on the photon energy. (code to produce this figure, data)

Figure 1 shows again the photon flux density from one of the last posts, but this time plotted with the epslatex terminal. The label size is changed by setting it to \footnotesize with the following code. First we introduce a abbreviation for the font size by adding a command definition to the header of our latex file.

set terminal epslatex size 9cm,7cm color colortext standalone header \
   "\\newcommand{\\ft}[0]{\\footnotesize}"

After the definition of the abbreviation we can use it for every label we are interested in.

set label 2 '\ft $5$\,meV'         at 1.38,4e9     rotate by  78.5 center tc ls 1
set label 3 '\ft $10$\,meV'        at 1.24,2e10    rotate by  71.8 center tc ls 2
set label 4 '\ft $20$\,meV'        at 1.01,9e11    rotate by  58.0 center tc ls 3
set label 5 '\ft $40$\,meV'        at 0.81,1e15    rotate by  43.0 center tc ls 4
set label 6 '\ft $60$\,meV'        at 0.76,9e16    rotate by  33.0 center tc ls 5
set label 7 '\ft $80$\,meV'        at 0.74,2.5e18  rotate by  22.0 center tc ls 6
set label 8 '\ft $E_0 = 100$\,meV' at 1.46,5e18    rotate by -40.5 center tc ls 7

August 22nd, 2013 | 12 Comments

On the PGF plots page I found a nice example of visualizing data with cubes. Here we will recreate the same with gnuplot as you can see in Fig. 1.

Cubes

Fig. 1 Cubes with different colors. (code to produce this figure, cube function, data)

We need basically two things in order to achieve it. First we have to plot a single cube with gnuplot. This is not as straight forward as you may thing. We have to define a data file for it and plot it with the pm3d style which will result in Fig. 2.

# single cube
0 0 0
0 0 1
0 1 1
0 1 0
0 0 0

1 0 0
1 0 1
1 1 1
1 1 0
1 0 0

0 0 0
1 0 0
1 1 0
0 1 0
0 0 0

0 0 1
1 0 1
1 1 1
0 1 1
0 0 1
set cbrange [0.9:1]
set palette defined (1 '#ce4c7d')
set style line 1 lc rgb '#b90046' lt 1 lw 0.5
set pm3d depthorder hidden3d
set pm3d implicit
unset hidden3d
splot 'cube.txt' u 1:2:3:(1) w l ls 1

The use of the fourth (1) column for the splot command ensures that the cube gets the same color on every side. Only the edges are highlighted by a slighty different color given by the line style.

A single cube

Fig. 2 A single cube. (code to produce this figure, data)

In a second step we reuse the code from the Object placement using a data file entry in order to plot cubes at different positions with different colors. To get the different colors and positions we replace the cube.txt file with a cube function that returns the values for the desired position and color.

add_cube(x,y,z,c) = sprintf('cube(%f,%f,%f,%i) w l ls %i,',x,y,z,c,c)
CMD = ''
stats 'cube_positions.txt' u 1:(CMD = CMD.add_cube($1,$2,$3,$4))
CMD = 'splot '.CMD.'1/0 w l ls 2'
eval(CMD)

July 30th, 2013 | No Comments

If you have more than one line or data set in your figure, they has to be named somehow. I’m not a big fan of a legend but prefer to put the names directly to the corresponding curves. That makes it easier for the reader. But as an author of the figure you have to find space to place the labels in the figure, and it could be that you have to rotate the labels to stick them to the lines.

Photo density flux

Fig. 1 Photon flux density for different characteristic tail state energies E0 dependent on the photon energy. (code to produce this figure, data)

Figure 1 shows the theoretical curves of photon flux density dependent on the photon energy for different characteristic tail state energies E0. E0 is an indicator for disorder in a crystalline system. To reduce the amount of trail and error for placing the E0 labels, we get the rotation directly from the data by fitting a linear function to the corresponding part of the data.

It is a little bit tricky, because we have a logarithmic y-axis. This can be handled by applying the logarithm to the y data by log($2) and than do the linear fiting.

f(x) = a*x+b
fit [1.30:1.34] [*:*] f(x) 'PL_spectrum_mu_1.0eV_E0_05meV_300K.dat' \
    u 1:(log($2)) via a,b

The first bracket is the range on the x-axis and the second sets the corresponding y range to auto. After we have the value of the slope we convert it to a rotation in degree with the r() function and set our label.

set label 2 '5 meV' at 1.38,4e9 rotate by r(a) center tc ls 1

The conversion to rotation is the most tricky part, because it depends on the range of your axes and the ratio between them. it would be trivial if both axes would go from, for example 1:10 and are equal in length in the figure. Otherwise we have to correct for the ranges and ratio.

# get the relation of x- and y-range
dx = xmax-xmin
dy = log(ymax)-log(ymin)
s1 = dx/dy
# get ratio of axes
s2 = 5/6.7
# helper function for getting the rotation angle of the labels in degree
deg(x) = x/pi*180.0
r(x) = deg(atan(s1*s2*x))