Gnuplot tricks

Plotting data is one of the main activities (in terms of time employed) of any researcher in modern-day science. If you ask any scientist what their favorite plotting program is, you'll likely receive an abundance of different answers depending on the field and also on how familiar with coding the person is.

Many scientists will agree with me that plotting data is as much of an art as it is anything else. And, as it happens when it comes to art, there is no definitive answer as to what is the best tool. My favorite plotting program is gnuplot. Gnuplot is a command-line program for plotting graphics that offers incredible flexibility. It's ideal for quick testing of results (almost everybody I know does quick result verification with gnuplot) but also extremely powerful for generating beautiful production-ready figures. This wide usage spectrum is what has always fascinated me about gnuplot. While virtually any user can get started relatively easily with quick plotting, I have to agree with the complaints of many people that claim there is a relatively steep learning curve before nice graphs can be produced. For this reason, even people who have been using gnuplot for many years are not aware of the most advanced options, which include scripting, statistical analysis and seamless integration with Latex, just to name a few.

During the years I have come up with a few gnuplot scripts that achieve quite nice results. I have posted some of these solutions as answers on the Stack Overflow Q&A website. Others I simply use routinely to produce figures for my publications. Hopefully they'll come in handy for those willing to engage an upper gear in the world of gnuplotting!

Passing Python (or any language) functions to gnuplot

This is one of the tricks I like the most. I originally posted it as an answer on Stack Overflow.

Sometimes one needs to define a complicated function. This function could be defined numerically and be the result of some program coded in some computing language. It might then be difficult or too involved to define such a function directly in gnuplot. However, one can still make use of that functionality by calling the program from within gnuplot. The idea is that one defines a function in gnuplot which is the result of a system call, where the argument is passed by gnuplot to the external program which then returns the result. If this program is written in Python, such a procedure would be constructed as follows.

Create a Python script with the function definition, e.g. called "function.py", within the same directory where gnuplot will be executed:

import sys
x=float(sys.argv[1])
print x**2


Carry out the function definition within gnuplot and use at will:

f(x) = real(system(sprintf("python function.py %g", x)))
plot f(x)


This will be much slower than "plot x**2" in the present very simple example, but the trade off can be justified if your function is very complicated.

Of course, this procedure also works with any other program (written in Fortran, C, bash, etc.) than can be executed taking arguments from the command line. In bash, using bc, this will do exactly the same as the Python example above:

f(x) = real(system(sprintf("echo \"(%g)^2\" | bc -l", x)))


To illustrate a bit better how powerful this capability can be, imagine you want to define a function which is given by the integral of some other function from zero to x:

$F(x) = \int_0^x \text{d}t \,f(t) = \int_0^x \text{d}t \, (\cos(t^2 + \sin(\text{e}^t-t^3) + \frac{1}{t})^2$

This is what $f(t)$ looks like (so imagine finding a primitive for it!):

We can write the following Python code ("function.py") to perform the integral (if you're not feeling lazy, you could write some faster code in Fortran or C):

import sys
import numpy as np
from scipy import integrate

x=float(sys.argv[1])
n=int(sys.argv[2])

(np.cos(t**2 + np.sin(np.exp(t)-t**3) + 1./t))**2, 0, x, limit=n,
full_output=1)

print Fx[0]


The integral $F(x)$ will be calculated with n maximum subintervals and I have also switched off error/warning messages (full_output = 1). We can make the function in gnuplot depend on n and check the effect of precision on the results (I'm skipping zero in the xrange to avoid division by 0 when the $1/t$ term is evaluated):

f(x, n) = real(system(sprintf("python function.py %g %i", x, n)))
set samples 20 # Reduce number of samples to make it quicker
plot[1e-5:5] f(x,2), f(x,4), f(x,8), f(x,16)


As expected the quality of the integral improves as the maximum number of subintervals is increased.

Customizing ranges (also color ranges)

Often, a linear or logarithmic scale is not adequate to convey the desired message. This can be due, for example, to the fact that small variations of a function are important within a subdomain but not outside of it. Consider the following function:

$f(x) = \text{e}^{-x^2} (1+10^8 \text{e}^{-\frac{(x-5)^2}{10^{-2}}})$

The function above has a feature at $x = 5$ which will go unadverted if we plot the function carelessly:

set samples 1000
f(x) = exp(-x**2)*(1.+1.e8*exp(-(x-5)**2/1.e-2))
plot f(x)


But perhaps, for some reason, we want to make sure that attention is paid to this feature. One way to do it, is to modify the y axis such that the range encompassing, say, $y = 0$ to $y = 0.002$ (where our feature is confined) occupies a significant portion of the graph. A very flexible way to do this is by using a rescaling function and modifying the labels on the y axis. What we'll aim to do in this example is to make the y axis vary linearly between 0 and 0.002, occupying the lower half of the axis, and then make it vary linearly between 0.002 and 1, occupying the upper half of the axis. Our rescaling function is obtained using a condition expression which maps the original ranges to a new range that varies between 0 and 2, where the new 0 corresponds to the original 0, the new 1 corresponds to the original 0.002, and the new 2 corresponds to the original 1. We then need to make sure that the ylabels are correctly set:

rescale(y) = y <= 0.002 ? y/0.002 : 1 + (y - 0.002)/(1. - 0.002)
set ytics (0)
set for [i=1:5] ytics add (sprintf("%g",i*0.002/5.) rescale(i*0.002/5.))
set for [i=1:5] ytics add (sprintf("%g",i*1./5.) rescale(i*1./5.))
plot rescale(f(x))


Color range customization

The same can be achieved using color ranges:

reset
rescale(z) = z <= 0.002 ? z/0.002 : 1 + (z - 0.002)/(1. - 0.002) # I'm just changing z <-> y to avoid confusion, but not necessary!
set pm3d map
set isosamples 1000
set size ratio -1
set xrange [-5:5]
set yrange [-5:5]
f(x,y) = exp(-(x**2+y**2))*(1.+1.e8*exp(-((x-5./sqrt(2.))**2+(y-5./sqrt(2.))**2)/1.e-2))
set cbtics (0)
set for [i=1:5] cbtics add (sprintf("%g",i*0.002/5.) rescale(i*0.002/5.))
set for [i=1:5] cbtics add (sprintf("%g",i*1./5.) rescale(i*1./5.))
splot rescale(f(x,y))


More examples on Stack Overflow

These are links to answers I have provided on Stack Overflow that make use of the rescaling trick above: