Thursday, January 19, 2012

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).

9 comments:

  1. This was really helpful. Thanks!

    ReplyDelete
  2. I've got a question for you. When I run my script, it outputs a graph with colored points (colors vary according to the third variable). How did you blend your graph?

    ReplyDelete
    Replies
    1. I had an absurd amount of points, so I didn't do any blending they just overlapped smoothly. I may have had to play with the dot size so there wouldn't be any obvious spaces.

      Delete
  3. Hi,
    nice tutorial ! quite usefull !

    I have a question, my data don't fill the entire sky. How could I add a grid which fitted this transformation to show the complete sky ?

    Thank you anyway

    ReplyDelete
    Replies
    1. Hey, glad you found this usefule, I'm not sure I understand what you need. I did do another post on adding a grid to this graph: http://astrojct.blogspot.com.es/2012/01/plotting-grid-on-projected-gnuplot.html hope that helps

      Delete
    2. Yes ! exactly, that is what I meant.
      Thank you for you posts

      Delete
  4. Hi,
    Very useful tutorial!

    I only get confused with your angles definitions. The first column in your data is the azimuth angle in degrees (0:360), right? I didn't get the range of the second column...

    For the Mollweide projection's functions, should I use radians or degrees?

    ReplyDelete
  5. how I can generate data file for the elliptical projection.

    ReplyDelete