Sunday, 26 July 2009

Maps - Contour plots with labels

So, you have always wondered how on Earth one can make a real map with gnuplot. Well, there is a simple and a not-so-simple way to this. First, let us see the simple way. Since it is simple, this method won't have one of the features, isoline labelling, that the second one has. As we go along, the map will become more and more complicated, but I hope that I set the right pace, and it will be easy to follow.

A map is nothing but a colour-coded 3D plot, with the isolines attached to it. We will produce a simple map using the function
sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
I don't think that this function has any particular meaning, but it looked quite all right, at least, as far as maps and isolines are concerned. If you have a matrix to plot, you can replace this function with that file. In principle, we could plot both the colour-coded map and the isolines from gnuplot, but we will have much greater flexibility, if we first direct the plots to a file, and then call the data from those files. One of the advantages of doing this is that in this way, we can make sure that the various plots have the same size. This requires a small overhead in terms of scripting, namely, we have got to issue the command
set table 'something.dat'
splot something
unset table
While this might appear superfluous, there are good reasons to do this. If we plot our function through the file 'something.dat', then we can use the 'with image' modifier of the plot command, and this means that the plot will be of the same size as would a normal 2D plot. Otherwise, if we use
set pm3d map
splot something with pm3d
the plot will actually be a bit smaller. The reason behind this is that in a 3D plot, we have to have some space for the 'z' axis, and even if we drop it in the map view, the space-holder for the 'z' axis is still there, therefore, gnuplot makes the whole plot a bit smaller. If, however, we plot the map through a file, we can use 'plot', in splot's stead, therefore, the 'z' axis will not appear anywhere in the processing of the plot.

After this interlude, let us see our first version of the map.
reset
f(x,y)=sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
set xrange [-5:5]
set yrange [-5:5]
set isosample 250, 250
set table 'test.dat'
splot f(x,y)
unset table

set contour base
set cntrparam level incremental -3, 0.5, 3
unset surface
set table 'cont.dat'
splot f(x,y)
unset table

reset
set xrange [-5:5]
set yrange [-5:5]
unset key
set palette rgbformulae 33,13,10
p 'test.dat' with image, 'cont.dat' w l lt -1 lw 1.5
In the first couple of lines, we plot the actual function into a file, called 'test.dat'. There will be some 25000 data points in this file, for there are 100 samples (default), and 250 isosamples. Having written these points into a file, we plot the contour. We just happen to know that the value of f(x,y) is between -3 and 3, so we set the contour levels at -3, -2.5...2.5, 3. When we are done with all this, we simply reset our plot, set the x and yrange, specify the palette that we want to use with the colouring, and call the plot command. For the actual data points, we use the 'with image' modifier, while for the contours, we increase the linewidth to 1.5, instead of using the default value of 1. We, thus, have just produced the following image.


Now, this was sort of standard, but the question inevitably comes up, whether we could do something more with this. When one has contour lines, one wants to label them, I presume. Petr Mikulik has a gawk script that will do this, in particular, it will place the appropriate labels next to the contour lines. I will discuss another approach here. The main difference with Petr Mikulik's method is that here we will put the labels on the contour lines, with some white space, of course.

First, it is quite instructive to look at the file containing the contour lines. In order not to cobble the plot, we will restrict it to the range [-5:0] (x) and [2:5], but it is really just for the sake of simplicity. With this modification, the file should read as
#Surface 0 of 1 surfaces

# Contour 0, label: 2
-0.482969 3.59036 2
-0.505051 3.58024 2
-0.509852 3.57831 2
-0.539895 3.56627 2
-0.555556 3.55994 2
-0.572136 3.55422 2
...

# Contour 1, label: 1.5
-1.11111 4.57693 1.5
-1.10874 4.57831 1.5
-1.0877 4.59036 1.5
-1.06617 4.60241 1.5
-1.06061 4.60546 1.5
-1.04273 4.61446 1.5
...
What this means is that the 0th contour line is drawn at z=2, and then the (x,y,z) values are listed. Of course, since we know that this particular line is at z=2, the last coordinate doesn't play any role, but it is listed all the same. When we plot this file, consecutive (x,y) points will be connected by a straight line segment. The next contour is at 1.5 etc. It might happen that one contour line is made up of several blocks, both for technical, and fundamental reasons. The technical reason is gnuplot's inner procedure of computing the an isoline, while the fundamental is simply that there might be several disconnected regions with the same isolevel. Anyway, if there are several blocks, they are separated by a blank line within the same contour.

In order to create the labels and the white spaces for them, we will use the following short script:
#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
{
if($0~/# Contour/) nr=0
if(nr==int(os+w/2) && d==0) {i++; a[i]=$1; b[i]=$2; c[i]=$3;}
if(abs(nr-os-w/2)>w/2 && d==1) print $0
nr++
}
END { if(d==0) {
for(j=1;j<=i;j++)
printf "set label %d \"%g\" at %g, %g centre front\n", j, c[j], a[j], b[j]
}
}' $1
This script operates on the data file of contour lines that we had before, and simply takes out a couple of points of the contours, while keeping one of the coordinates as the position of the labels. You can change the length of the white space by setting the 'w' argument, while the position of the white space can be modified by changing the 'os' argument. The first argument, 'd', determines what we want to do with the script: make the labels, or tweak the input file.
Now, we can call the script from our gnuplot script as follows:
reset
set xrange [-5:0]
set yrange [2:5]
set isosample 150, 150
set table 'test.dat'
splot sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table
set cont base
set cntrparam level incremental -3, 0.5, 3
unset surf
set table 'cont.dat'
splot sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table
reset
set xrange [-5:0]
set yrange [2:5]
unset key
set palette rgbformulae 33,13,10
l '<./cont.sh cont.dat 0 15 0'
p 'test.dat' w ima, '<./cont.sh cont.dat 1 15 0' w l lt -1 lw 1.5

There are two more twists that we could add to the figures above. One is that the labels can be rotated in such a way that they are parallel to the curves. This we can easily achieve by adding two more lines to our script above: we take the first and last dropped coordinate, and calculate the angle that that line segment would make with the horizontal. Using this angle, we rotate the labels accordingly. With this modification, the script now becomes
#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
{
if($0~/# Contour/) nr=0
if(nr==int(os+w/2) && d==0) {a[i]=$1; b[i]=$2; c[i]=$3;}
if(nr==int(os+w/2)-1 && d==0) {i++; x = $1; y = $2;}
if(nr==int(os+w/2)+1 && d==0) r[i]= 180.0*atan2(y-$2, x-$1)/3.14
if(abs(nr-os-w/2)>w/2 && d==1) print $0
nr++
}
END { if(d==0) {
for(j=1;j<=i;j++)
printf "set label %d \"%g\" at %g, %g centre front rotate by %d\n", j, c[j], a[j], b[j], r[j]
}
}' $1


and the resulting graph is shown here:

Note that the only difference between this and the previous script is the following two lines

if(nr==int(os+w/2)-1 && d==0) {i++; x = $1; y = $2;}
if(nr==int(os+w/2)+1 && d==0) r[i]= 180.0*atan2(y-$2, x-$1)/3.14
where we calculate the appropriate angles.

The second twist is that it might make interpretation of the figure easier, if the contours are coloured, as are the corresponding labels. We can use the script above with some small modification. The key is to use the 'index' modifier. Recall that the contours are rendered in blocks, and with the help of the 'index' keyword, we can specify which block we want to plot. We, therefore, modify our script as follows:
#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
{
if($0~/# Contour/) nr=0
if(nr==int(os+w/2) && (d%2)==0) {a[i]=$1; b[i]=$2; c[i]=$3;}
if(nr==int(os+w/2)-1 && (d%2)==0) {i++; x = $1; y = $2;}
if(nr==int(os+w/2)+1 && (d%2)==0) r[i]= 180.0*atan2(y-$2, x-$1)
if(abs(nr-os-w/2)>w/2 && (d%2)==1) print $0
nr++
}
END { if(d==0) {
for(j=1;j<=i;j++)
printf "set label %d \"%g\" at %g, %g centre front rotate by %d tc lt %d \n", j, c[j], a[j], b[j], r[j], j

}
if(d==2) {
printf "plot \"test.dat\" w ima, \"cont.plt\" index 0 w l lt 1,\\\n"
for(j=2;j<i;j++) printf "\"\" index %d w l lt %d,\\\n", j-1, j
printf "\"\" index %d w l lt %d\n", i-1, i

}
}' $1

Calling the script produces the following graph:



Note that we had to call plot through a new temporary file, cont.plt, which should be produced by redirecting the output of
cont5.sh 2 15 0 > cont.plt

44 comments:

  1. Let me salute You, Great Gnuplotter. And, I'd like show my gratitude to you, because of the great help I found in the bible you keep here.
    The reason for which I'm writing could be seen as heresy, but maybe you could take me back to the right path.
    I found some problems applying this gawk script. In particular, I had to add one line and change some conditons to make it work properly. Please understand that I'm a student, giving my first steps in Linux platforms as well as with gnuplot. I will be very gratefull if you contact me to my email address, erchamionshinobi@gmail.com, so I can expand the explenation of this problem, and get it to be more fruitfull.
    Thankyou again. I'll be waiting forward for your answer.

    ReplyDelete
  2. Hello,

    thank you very much for your script for labelling contourplots. But, however, I found some problems when I applied. When I repeat your example everything is fine, but I m doing some plots for my work and I do not get blank lines in the contourplots, that is, the script does not leave a blank line where the tag should go. On the contrary the tags appear on top of the contourlines.

    To solve I have to modify manually the contour.gp file and insert a empty line. I do not know the awk language so maybe you could help me.
    I could send my example.

    Thank you very much,
    Susana

    ReplyDelete
  3. Hello Susana,

    If you give me your e-mail, I will try to help you out with this.
    Cheers,
    Zoltán

    ReplyDelete
  4. See also the post on the 10th Feb 2010. Cheers,
    Zoltán

    ReplyDelete
  5. This is really useful! Thanks for posting this.

    I wanna ask, is there a way to adjust the position of the labels on the contours? In my own contour plot, the labels on the contour clutter closely and are hard to read.

    ReplyDelete
  6. If you want to use this version, then you can specify the third argument in the script, it is called os. That basically determines by how my the labels are shifted with respect to their zero position. You have to try a couple of values, and see which one fits your needs.
    You can also look at another way of plotting maps; it was posted on the 10th of February.
    Cheers,
    Zoltán

    ReplyDelete
  7. Thank you very much for this site. I have been using Gnuplot for many years, and was frustrated at some of its apparent limitations. Thanks to the info you have posted, I have learned that Gnuplot is capable of much more than I realized.

    Today, using what I learned on this post, I hacked out a python script that is called from my gnuplot script, to automatically locate plot labels at just the right location and angle. Before I was creating the label locations and angles manually, and it would take many iterations to get them properly located.

    Kevin Horton

    ReplyDelete
  8. Greetings Kevin,

    Thanks for visiting, and I am glad that you found something useful here. Also, if you want to plot maps, you could check out a recent post that is based on gnuplot 4.4. I believe that things are much simpler there.
    Cheers,
    Zoltán

    ReplyDelete
  9. extraordinary figure - already helped me a lot. thanks!

    ReplyDelete
  10. There is a small problem with the awk script, because when abs(nr-os-w/2) > w/2 it doesn't print a blank line. So, if the contour presents open lines, the gap doesn't appear. Just by adding a line like:
    if(abs(nr-os-w/2)<w/2 && (d%2)==1) print ""
    the gap opens tuning the os parameter.

    In any case, it's an excellent post. Congratulations.

    ReplyDelete
  11. Soooooorry, I meant ... when abs(nr-os-w/2) < w/2 ...

    ReplyDelete
  12. Hello whitesnow,

    Thanks for pointing out the error! I would, however, urge everyone to use gnuplot 4.4 and the corresponding script for producing the contour lines, because it is much easier and more straightforward. I have simply given up on this (g)awk-reliant version, because it is harder to maintain (as I pointed out at various places, this is not fool-proof, and requires tweaking on a case-by-case basis), and my interest shifted towards doing everything in gnuplot, and not using external scripts.
    Cheers,
    Zoltán

    ReplyDelete
  13. Thanks Gnuplotter, plots look great.

    And it gave me an idea how to do that without awk, which I'm not familiar with.

    The functionality (contours with labels) could be achieved purely in Gnuplot:

    # plot contours info into a file
    set table 'cntrs'
    splot <... whatever ...>
    unset table

    # plot countours as usual, plus as labels
    # based on the value of contour level,
    # every ::1::1 - selects just one label per
    # contour section

    plot 'cntrs' w l, '' using 1:2:(sprintf("%.2f", $3)) every ::1::1 w labels point pt 1 ps 0.5 offset x,y


    # For coloured lines -
    plot w l palette z, ....

    You mentioned in comments that you've shifted to doing everything in Gnuplot 4.4, which I'm a fan of (too lazy to learn awk), so may be you would like to update this page and share your new tricks.

    Thanks for tutorial!

    ReplyDelete
  14. Hi Gnuplotter,

    Found your updated script - thanks again!

    Still I reckon using "every" modifier is easier and basically does not require any pre-processing of data file at all.

    Location of the labels on the lines could be adjusted to taste by changing the record number in "every ::N::N". Moreover - selecting high numbers will automatically skip short contours that have less points than "N" in their isolone. This helps to avoid clutter around sharp minima/maxima.

    ReplyDelete
  15. Hi gnuplotter,
    Thank you so much for all these tricks.

    I would be greatful if you could provide the way to label contours in windows paltform. I am not much handy in linux...so waiting your answer

    ReplyDelete
  16. Hi,

    I am not sure how difficult it is under Windows. The problem is that you have to use an external script that does some data processing, and then call gnuplot again. Basically, if you can manage to set up gawk, then you could just put everything in the .sh file into a windows-executable batch file, and then proceed in the same way as you would in Linux.
    I am afraid, this is all I can help.
    Cheers,
    Zoltán

    ReplyDelete
  17. Hi. Congrats for your post. Is there any way to put some data points upon this graphic?

    ReplyDelete
  18. Gnuplotter

    I see your examples above refer to rectangular maps. If I want to do similar to a non rectangular map, is it possible? I cant make it work ..

    Geoff

    yesgeoff@gmail.com

    ReplyDelete
  19. I ran your example and I get this as an error

    gnuplot qso_ratios_contour_2.p
    sh: ./cont.sh: Permission denied
    sh: ./cont.sh: Permission denied

    p 'test.dat' w image, '<./cont.sh cont.dat 1 15 0' w l lt -1 lw 1.5
    ^
    "qso_ratios_contour_2.p", line 20: warning: Skipping data file with no valid points


    Any ideas on what's going wrong here?

    ReplyDelete
  20. Hi Rafael,

    I don't know what you mean by that. What should that look like?
    Zoltán

    ReplyDelete
  21. Hi Geoff,

    I believe, a non-rectangular configuration is going to be hard.
    ZOltán

    ReplyDelete
    Replies
    1. I'm do it by passing five parameter to sh-script
      'ratio'=(GPVAL_Y_MAX-GPVAL_Y_MIN)/(GPVAL_X_MAX-GPVAL_X_MIN)
      In sh-script I'm use it in following way:
      r[i]= 180.0*atan2((y-$2)/ratio, x-$1)/3.14

      Delete
    2. PS It seems variables GPVAL_* has initialized after calling splot at proper xrange and yrange

      Delete
  22. I ran your example and I get this as an error

    gnuplot qso_ratios_contour_2.p
    sh: ./cont.sh: Permission denied
    sh: ./cont.sh: Permission denied

    p 'test.dat' w image, '<./cont.sh cont.dat 1 15 0' w l lt -1 lw 1.5
    ^
    "qso_ratios_contour_2.p", line 20: warning: Skipping data file with no valid points


    Any ideas on what's going wrong here?

    Yes: the problem is that you didn't set the permission for the script. Run

    chmod a+x script.sh

    on the command line.
    Zoltán

    ReplyDelete
  23. hi,

    I want to plot a string instead of lines or dots or crosses or steps. is it possible?

    ReplyDelete
  24. Yes, it is possible. You should just look for the stringcolumn option in the manual.
    Cheers,
    Zoltán

    ReplyDelete
  25. This comment has been removed by the author.

    ReplyDelete
  26. This comment has been removed by the author.

    ReplyDelete
  27. Nice posting. it is very helpful

    ReplyDelete
  28. Hello
    I will do the sometging with this function :

    f(x,y)=(4-2.1*x**2+x**4/3)*x**2+x*y+(-4+4*y**2)*y**2
    But I did'nt obtain colors
    Thanks Mariam

    ReplyDelete
  29. Is there a way to control the font in the labels? I have tried to put font "Lucidabright,xx" (where xx is most everything reasonable) into cont.sh (the gawk script), but this simply seems to be ignored. What I really want to do is to have label fonts and axis fonts different (in size)...

    Great stuff, these plots. I also like your choice of colormaps.

    ReplyDelete
  30. Here you have another possibility without using gawk (which I have never used):

    n= 30
    plot 'datos.txt' using ($1):($2):($3) with image,\
    'cont.txt' using ($1):(($0)<(n-3) || ($0)>(n+3) ? ($2):1/0):($3) notitle with lines lc rgb "black" lw 1,\
    'cont.txt' using ($1):(($0)==n ? ($2):1/0):($3) notitle with labels font "Times,6"

    using n=30, you plot the value number 30 of each contour and open a breach in the contour lines at that point. The breach has 7 points in size
    Y para los que lean español: Usando n=30, se puede colocar una etiqueta para cada línea de contorno en los huecos (de 7 valores cada hueco) que se han abierto en los puntos nº 30 de cada bloque.

    ReplyDelete
  31. Cool stuff, I have a particular scenario where I am trying to plot a contour map on top of an image plot. It looks almost correct except for a slight misalignment of the 0,0 coordinates of both plots. The command lines are below (using gnuplot 4.2.6). I'd like to somehow get you the data, "test13.txt" but would need an email address to do so. I can be reached at cmr72@hotmail.com. Thanks!

    gnuplot -geometry 750x500+1280+0 -persist <<!
    set xrange [-50:50]
    set yrange [-50:50]
    set contour base
    unset surface
    set view 0,0
    set table 'test13c'
    splot "test13" using 1:2:3 w lines notitle
    unset table
    set colorbox vert user origin 0.05,0.05
    set size square
    set palette rgb -34,-35,-36
    plot "test13" using 1:2:3 w image notitle, \
    "test13c" i 0 notitle with lines ls 1, \
    "test13c" i 1 notitle with lines ls 2, \
    "test13c" i 2 notitle with lines ls 3, \
    "test13c" i 3 notitle with lines ls 4, \
    "test13c" i 4 notitle with lines ls 5
    !

    ReplyDelete
  32. I want to put a label below a plot. I have already set xlabel but I want new label to be even below that. Actually I would love if you could tell me to make key visual in pm3d plot. Because key is not visible in pm3d plot now I am searching for this option. My email id is snehalshekatkar@gmail.com. Can you please help me?

    ReplyDelete
  33. This comment has been removed by the author.

    ReplyDelete
  34. I cant produce coloured contour lines although the labels are coloured. I can be contacted on tafie.tafa@yahoo.com. Here is my gnuplot script

    set terminal po eps enh color "Helvetica" 20 lw 1
    set output "Ta.eps"


    reset
    set xrange [0:351]
    set yrange [0:351]

    #set isosample 15, 15
    set table 'test.dat
    splot "./data_elastic_field.15.dat"
    unset table
    set cont base
    #set cntrparam level incremental -3, 0.5, 1
    set cntrparam levels auto 30
    unset surf
    set table 'cont.dat'
    splot "./data_elastic_field.15.dat"
    unset table
    reset
    set xrange [0:351]
    set yrange [0:351]
    unset ytics
    unset xtics
    set cbrange [-0.7:1.]
    unset key
    set size square
    set palette defined (0 "green", 1 "blue", 2 "red", 3 "orange" )
    l '<./cont.sh cont.dat 0 15 200'
    p 'test.dat' w ima, '<./cont.sh cont.dat 1 1 0' w l lt -1 lw 0.1

    ReplyDelete
  35. Hi everyone, I ended up here after searching also a non-awk reliant solution for the insertion of blank lines to use splot, and I see that many of you had the same question. Could anyone please tell me how I could obtain the solution? (The problem is that I am completely new with linux, and I am using a Mac OS X system).
    Thank you very much in advance, Tim (email timwillaert4@gmail.com)

    ReplyDelete
  36. I've managed to introduce the blank lines where I need them with the awk script, however I still get an error 'All points x value undefined', but the data file just contains 300 blocks of 300 (x,y,f(x,y)) values. Any ideas what might be the reason? Many thanks, Tim

    ReplyDelete
  37. Hi thanks for the blog it is very useful,

    I am using the script
    gnuplot> set term postscript eps color enhanced
    gnuplot> set output "FES2.eps"
    gnuplot> set pm3d at b
    gnuplot> set dgrid3d 170,170
    gnuplot> set view 70, 50,1,1.3
    gnuplot> splot './fes2.dat' w dots

    Now i am getting a default red mesh as the contour
    can you suggest me how to change the mesh color

    ReplyDelete
  38. I have a similar plotting problem, only I would like the end result to be a 2D color map in POLAR coordinates. I have 3D data that consists of the SPL versus frequency measured at a fixed distance away from a loudspeaker. The set of measurements are taken at the same height above the floor, all the way around the loudspeaker (so native polar coordinate system data). The position directly in front of the loudspeaker is considered "on axis" or the angle = 0 there.

    I would like to produce a polar plot with the frequency on a log scale (with a lower limit cutoff, e.g. at 100Hz) as the radial coordinate, the angular coordinate equal the the angle off axis at which the measurement was done. A color map within the polar plot (circular region) represents the SPL level at that frequency and off-axis angle.

    I mostly understand how the process shown above is done to first produce the file for the contours, write to disk, then read it back in as an image for the 2D plot. But I am not sure how to approach this for a polar coordinate system. Do I use polar coordinates in 3D when producing the contour file as well? I have seen some hints on doing that here:
    http://stackoverflow.com/questions/16156095/how-to-create-a-3d-polar-graph-with-gnuplot
    Then I can just plot up the contour file as a 2D polar plot?

    Looking for a sanity check here and some guidance or thoughts from others. Thanks!

    ReplyDelete
  39. Hi Gnuplotter,

    Thanks for the fantastic blog. I have followed your content to make a simple contourplot using the following code :

    f(x,y)=(((x**2)+(y)-11)**2)+(((x)+(y**2)-7)**2)
    set xrange [-5:5]
    set yrange [-5:5]
    set contour base
    set cntrparam level discrete 13.59085,25,50,100,150,300,500,1000
    set table 'cont.dat'
    splot f(x,y)
    unset table
    reset
    set xrange [-5:5]
    set yrange [-5:5]
    unset key
    p './cont.dat' w l lt -1

    There are a few problems with the plot:
    1) It shows some horizontal lines in the plot, I don't know why?
    2) the contour lines are not too smooth. How can I improve there smoothness?

    Can you please help me understand the problem?

    Thanks

    Abhinav

    ReplyDelete
  40. Hi Gnuplotter,

    Excellent article, really helpful and fairly easy to follow.

    I'm trying to implement your idea in a "double-contour plot", i.e. plotting two sets of contours representing variables z1 and z2 over the same x-y plane. And I'd like to label both sets of contours using your method. But when I try to load the two sets of labels using two different datasets with the commands
    load '<./cont.sh cont1.dat 0 2 4'
    load '<./cont.sh cont2.dat 0 2 4'

    it doesn't seem to work (the labels appear only on the contours representing the second dataset, so I guess the second 'load' command must be overwriting the first). Could you suggest a way around this?

    Thanks a million!

    Amol

    P.S.: I know I could just make two separate plots and overlay them using the multiplot environment, but I'm looking for something less crude.

    ReplyDelete
  41. Hi Gnuplotter,

    I want the same color for all labels (white for example) and change the font size. How can i do it ?

    ReplyDelete