Importing MOLA data into GMT

The MOLA gridded data set is a great resource. The Generic Mapping Tools (GMT) are a great tool set. This article shows you how to take the raw binary MOLA gridded data set and convert it into the netCDF format that GMT uses. Once you do that, a whole range of mapping options are open to you.

We will assume that the Generic Mapping Tools (GMT) are already installed on your system. If not, they are an open source tool, and you should be able to get them installed relatively easily.

After that you need to get some gridded MOLA data. As always, the best resource for NASA mission data is the Planetary Data System (PDS), and the gridded MOLA data we want can be found on the Geosciences Node.

The method that I'll discuss can be applied to any of the kinds (areoid, counts, radius, and topography) and resolutions (4, 16, 32, 64, and 128 pixels per degree, more for the poles) of MOLA gridded data sets. For this example, we'll use topography at 128 pixels per degree that has Valles Marineris it it. The megt00n270hb files are the ones we want. I know this because I know the latitude and longitude of the feature I'm interested in, and I have familiarized myself with the scheme for how to decypher the alphabet soup of a file name (which can be found on the gridded data page).

So you need the raw gridded data, which is a file with the img suffix, and it is also handy to have the label file (it has the same name but with a lbl suffix).

So I'll cut right to the chase, the next step is to type the following on your command line, all in one go (the text below may be wrapped by your browser):

xyz2grd megt00n270hb.img -Gmegt00n270hb.grd -I0.0078125 -R270/360/-44/0 -Ddegrees/degrees/m/1/0/"MOLA Topography"/"128 pixels per degree" -F -ZTLhw -V

That's really all there is to it. Once you do that, you'll have a netCDF file, megt00n270hb.grd, that is ready to be used by other GMT programs.

However, if you already knew all of the ins and outs of GMT, then you probably wouldn't be reading this. So here are the details for each of the above elements:

xyz2grd
The GMT program that we'll be using to convert the raw xyz data into a netCDF gridfile, its manual page is available via the excellent GMT online documentation. It might help to have that up in another window, and follow along. Depending on where the GMT programs are installed, you may have to be more explicit with the path to the program.
megt00n270hb.img
Our input MOLA gridded data file.
-Gmegt00n270hb.grd
The output file argument. Unfortunately, GMT doesn't allow spaces between its option keys and their values. This is the option -G with the value megt00n270hb.grd. This is the name that we'll call the output file. It could be anything we want, but this follows some reasonable naming patterns.
-I0.0078125
The grid spacing for our grid is 0.0078125. Since we know that there are 128 data values per degree in our input file, we know that the spacing between those datapoints in degrees is 1 ÷ 128 = 0.0078125.
-R270/360/-44/0
These are the bounds of the data in our input file. They are in the order west, east, south, and north separated by slashes. You can figure them out from the handy data page where you got the gridded data files, but the numbers are also in that label file (megt00n270hb.lbl) which is just a text file.
-Ddegrees/degrees/m/1/0/"MOLA Topography"/"128 pixels per degree"
This is an optional argument for xyz2grd, but I like to put this information in. You can see what all of these slash-separated things mean on the xyz2grd man page.
-F
The MOLA gridded data is pixel-registered, not grid-registered, so we must use this option to let xyz2grd know that. An analogy would be a chessboard. It has eight squares on a side, but the grid that defines those squares has nine vertices on a side. If the "data" were chess pieces, taking the default grid-registered behavior would cause us to put the chess pieces on the intersections of the lines, not in the squares, and we'd run out of pieces before we ran out of vertices on a side.
-ZTLhw or -ZTLh
This quiet little string of letters does all the work of helping xyz2grd parse the binary data. The "TL" tells xyz2grd that the data in the binary file should populate the grid file as rows starting in the top left of the grid. The "h" indicates that the data in the img file should be read as short 2-byte integers (this can be figured out by reading that label file, in a few rare cases this may need to be a "u" instead of an "h"). The "w" tells xyz2grd to swap the byte-ordering, you may or may not need this. The label file tells us that the data (SAMPLE_TYPE) are MSB_INTEGER. The details of endianness are too complex for this article. If you are running on a big-endian computer (Sun systems or PowerPC chip Macs for example) then you won't need the "w". If you are running on a little-endian computer (the commodity PC hardware that runs Linux, FreeBSD, etc.), you will. If you're not sure ask your system administrator or just grab a gridded data file of an area that you know, and try it one way. If the topography looks weird, try it the other way.
-V
This is totally optional, but I often like to see what programs are doing, so I like using the verbose setting.

That's it, now you know how to take a raw MOLA gridded data file and convert it into a netCDF file.

GMT is great, but I get tired of typing in all of those options all of the time. I've created a Perl program, mola2gmt, that acts as a wrapper around the GMT xyz2grd program, and is smart enough to parse that label file to fill in a lot of the above options automatically. Furthermore, it also has the option of running grdcut on that file, since you rarely need an area as big as the default tile size of the gridded data products.

When you write up that great paper, don't forget to properly cite the PDS data. Similarly, you can also cite some GMT articles if you like.