Plotting a grid on a projected gnuplot graph (e.g. Mollweide)

Following the first post on how to plot a map on gnuplot, and then how to plot it with a Mollweide projection (all in the same post), I decided to dedicate a post on how to do a grid overlay.

Mostly because I first thought the best idea would be to simply plot a file with data points only on the meridians and parallels. This is a bit of a kludge and doesn't allow you to easily add styles to the lines.

The obvious solution: to plot the function that describes the lines of our grid. But one really nice thing I stumbled upon (in this nice blog called Gnuplot tricks) is using the "for" loop to plot a function n times:

f(x,i)=i*x

plot [-pi:pi] for [i=1:10] f(x,i)

Otherwise you'd have to plot each one of the lines. This way you only have to issue two or four commands (depending if you want a different pattern for a "minor grid").

We'll use the Mollweide projection parametric functions that we copied from Roberto's blog(the ones mentioned in my previous post). Now I almost lost my head trying to express these parametric expressions as normal functions so I could just write "plot y(x)" (the parametric ones depend on other variables and give you an expression of x,y pairs instead of a function y(x) ).

Enter the command: set parametric. This will allow you to do the same as with columns in a data file, which is to plot on function (y=f(t)) against another (x=(g(t))¹ .

So this is the code for a Mollweide grid:

set parametric

plot [-pi:pi] for [i=-5:5] mwx(pi*i/10,t),mwy(pi*i/10,t) with line linetype -1, \

for [i=-5:5] mwx(t/2,pi*i/5),mwy(t/2,pi*i/5) with line linetype -1

Note that we use "with linetype -1" for solid black lines and "with linetype 0" for dotted black lines. Gnuplot is programmed so that the simple way is to use predefined styles for lines/points/etc. However one can specify color with the somewhat obscured command "with lines linecolor rgb '#000000'" or abbreviated "w l lc rgb '#000000'".

Another very important thing to note is that the range of the dummy variable t is only set at the beginning of the plot command and cannot be issued multiple times (cannot be issued after the comma). I considered using multiplot for overriding this but it was a very bad idea when dealing with the map previously plotted.

Oh right, you'll probably want the map to be behind those lines, but the range of the plotted functions' dummy variable (t) must be issued right after the plot command, so in order to do the full plot of map and grid this is the final plot command:

plot [-pi:pi] "gnuplot.dat" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))):($3/1E7) with points palette, \
for [i=-2:2] mwx(pi*i/6,t),mwy(pi*i/6,t) w l lt 0,\
for [i=-3:3] mwx(t/2,pi*i/4),mwy(t/2,pi*i/4) w l lt 0, \
for [i=0:1] mwx(t/2,pi-2*pi*i),mwy(t/2,pi-2*pi*i) w l lt -1 lw 3, \
"grid-labels.dat" using (mwx($2/180*pi,$1/180*pi)):(mwy($2*pi/180,$1*pi/180)):3 with labels

I tried plotting the grid labels more elegantly but failed due to constrains in the system, so I ended up just plotting a data file that contained the coordinates and text for a label. Also note I decided to change the inner linetype to 0.

Here's the final result:

¹I tend to try in all possibilities to give a reference to what helped me, but in this case it was hours of droning with my favorite search engine until I stroke the right phrase that lead me to a gnuplot faq page that I didn't quite find explanatory but gave me the clue to search the word "parametric" (hoping for the manual but falling on this lanl webpage I keep bumping into).

Plotting a 2D map with colors for a 3rd column of values. (and Mollweide project it)

In case you're new to the scientific community, gnuplot is almost a standard tool in plotting scientific data. It's free (as in open source) very robust and extremely versatile. However this versatility might be to blame for it's lack of user-friendliness, there exist some GUIs (Graphic User Interfaces) but they're not robust and go against the concept of being able to control anything you want in your graphics (and also automating graphs, but that is for another post).

However sometimes you want to do something very specific, and this is the case for this post, I have some XYZ data in which X and Y are simple coordinates (galactic in this case) and Z is a value for each one of those points.

In my case I have a file with 196608 sets of points each with it's value, of course if you try to get any spreadsheet software to graph that it'll probably crash, make your computer slow or simply refuse to accept that many lines.

The data file looks something like this:

357.77228 52.416012 2.0093570

358.66337 52.416012 2.0057610

359.55446 52.416012 2.0027640

0.44117647 52.029727 1.9932720

Googling you might stumble upon solutions that rely on plotting "with p3dm", sadly these are not what you need, as this command is meant to plot grids and will therefore need sets of points that make up each rectangle, each separated by spaces. You'll probably get frustrated with the error "Hint: Missing blank lines in the data file? See 'help pm3d' ".

The final answer to my problem I found here on kiko's blog. Which is summarized in two lines:

gnuplot> set view map 

gnuplot> splot "test_map.dat" with points palette pt 9 

If you want to change the colors following the directions from this site(yes, they're using pm3d), you can just input the command:

set palette rgbformulae 22,13,-31

Play with those numbers all you want.

Note: I later found (at the Gnuplotting blog) a better, more pleasant gradient that is easier to modify if you know RGB hex codes (I didn't change it for this post so if you want to see it test it out):

set palette defined ( 0 '#000090',\

1 '#000fff',\

2 '#0090ff',\

3 '#0fffee',\

4 '#90ff70',\

5 '#ffee00',\

6 '#ff7000',\

7 '#ee0000',\

8 '#7f0000')

That is before splot (or after and then issue the command "replot").

And here's the final result for me (before adding titles but after tweaking xrange and yrange):

Some basic things for those new to gnuplot are:

set xrange[0:360]

set yrange[-90:90]

set title "Graph Title"

set xlabel "Latitude (b)[degrees]"

set ylabel "Longitude (l)[degrees]"

set key off

Well, I know this post (and probably many more to come) says nothing new, but in this age of over-information, I hope I blabbered on enough so that if you're looking for this specific problem it might just pop up in your online query (and also to place higher up in your search engine the place where I found the solution).

P.S.: Mollweide projection

Something very nice to look into is using a Mollweide projection, this can be done in gnuplot thanks to Roberto's function:

Just copy/paste this block of text into the gnuplot command line:

mwhigh(x) = 1.0 + -0.919061*(abs(1.0-x))**0.674635

mwmed(x) = mwhigh(x) -0.0807765 + 0.161136*x -0.0796311*x**2

mwlow(x) = mwmed(x) -3.53551e-05 + 0.000645749*x

mwst(x) = x < 0.2 ? mwlow(x) : (x < 0.9 ? mwmed(x) : mwhigh(x) )

mwt(b) = b>0 ? asin(mwst(sin(abs(b)))) : -asin(mwst(sin(abs(b))))

mwx(b,l) = 2.0*sqrt(2.0) * l * cos(mwt(b))

mwy(b,l) = sqrt(2.0) * sin(mwt(b))

and that'll mean you have the mwx and mwy functions to change your coordinates, I used them this way:

splot "gnuplot.dat" using (mwx($2,$1)):(mwy($2,$1)):3 with points palette

Note one thing: when you're just plotting a column in your data file, the number of the column will suffice (1:2:3), however, when you're using functions inside the plot or splot command you have to put the functions inside parentheses and the columns have to be called out with dollar signs ( (functionX($1):functionY($2):3 ). This took me quite a while to figure out because it gives you no errors and just generally misbehaves. Sadly I can't remember what tipped me off to this, but part of it was the gnuplot manual, which is big, but a must have reference.

You may have noticed that those equations will not work for my data since they're meant for coordinates that go from: [-pi/2,pi/2] for b and [-pi,pi] for l.

So we just need to change the coordinates adding two more functions:

bconv(b) = b*pi/180

lconv(l) = (l<180) ? l*pi/180 : (l-360)*pi/180

note that in the second one I had coordinates that went from 0 to 360, which is odd, but it was a good example of how to use conditionals.

and now we can change the plotting line to:

splot "gnuplot.dat" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))):3 with points palette

So the now beautiful result is:

note that for cosmetics (and since the axes mean very little now) we can remove the border, ticks and labels:

unset border

unset xtics

unset ytics

Well, I guess that's that for a long P.S. section, I hope it makes sense. If anyone wants homework: post in the comments how to overlay a couple of central axes with the correct values for longitude and latitude!

(Also feel free to ask for help or just let me know that this post somehow helped you, I'd love to know).