A collaborative community for planetary sciences.

Creating shaded relief maps or combination colorized/shaded relief maps with MOLA data is easy to do with the Generic Mapping Tools (GMT).

This article will go into the details of taking some MOLA gridded data that has already been converted into the netCDF format that GMT uses, and creating a nice shaded relief map appropriately color-coded by elevation.

We'll first just create a color-coded map, and then put on the shaded relief. The first thing that we need to do is create a color palette file for our map. Let's say that I have a netCDF MOLA gridded data file called mola.grd. GMT provides a program called grd2cpt which will read in a .grd file and create a color palette file that other GMT programs can use. The grd2cpt program has a number of nice options, but the first time through, I usually just run with the defaults like so:

grd2cpt mola.grd -C/home/rbeyer/.gmt/color.cpt -V -Z > mola.cpt

Here's the detailed breakdown:

grd2cpt

This is the GMT program that we'll be using. You may want to be looking at its manual page to follow along. Depending on where the GMT programs are on your system, you may also have to include a path to the grd2cpt program.

mola.grd

This is the input .grd file. The grd2cpt program reads this file and performs a mapping of the data through a cumulative distribution function so that the colors are histogram equalized. If you are only going to want to make a map from a small portion of this grid, then use the -R option for grd2cpt, so that its algorithm for distributing the colors only operates on the height ranges in your area.

-C/home/rbeyer/color.cpt

Most people might just use one of the prepackaged color tables that comes with GMT, but if you have a custom color table, for example in a file at /home/rbeyer/color.cpt, you can use it this way. Not specifying the -C option will give you a default rainbow color scheme (ugh!). See this article for advice on color schemes.

-V

I always use the verbose option to see what GMT is doing.

-Z

This will create a continuous smooth color range. Sometimes you want this, sometimes you don't.

> mola.cpt

This isn't part of the program, but grd2cpt doesn't have an option for an output file, it just writes to standard out. So to capture the results into a file that can be used by other GMT programs later, you have to capture that output (here with the > shell redirect operator) into a file.

Now that we have that color palette file we can make our first map. Which can be done with the following:

grdimage mola.grd -Cmola.cpt -JI305/7i -Bf1a2 -R301/310/-17/-12 -V > mola_map.ps

grdimage

The GMT grdimage program will create the PostScript map that we desire. If you want to create Encapsulated PostScript instead of just regular PostScript, you'll need to look into altering your GMT defaults and add a plus (+) to the specification of your preferred PAPER_MEDIA value.

mola.grd

This is the name of our grid file with our elevation data

-Cmola.cpt

This is the color palette file that we made above with grd2cpt.

-JI305/7i

This is the map projection. You can read the grdimage man page to learn more, however, I direct your attention to the fact that the "I" is capitalized. Both -Ji and -JI specify a sinusoidal projection. The difference is how grdimage interprets the scale parameter after the slash. Make sure you understand the difference, and use the one you mean. I almost always use the uppercase version.

-Bf1a2

This turns on map boundary annotation. It puts a mark every degree, and a label every two degrees.

-R301/310/-17/-12

This is the region of the whole grid file that we want on our map. If we want the whole thing, you don't need to use this option at all.

-V

I like seeing what it's doing.

> mola_map.ps

Again, the output of grdimage is sent to STDOUT, so this captures that into a file called mola_map.ps.

This should give you a map in mola_map.ps that is color-coded by elevation. Often, I only want to concentrate on a narrow range of elevations with my colors. So you may want to go back to the color palette creation step and experiment with the -L option to grd2cpt to restrict the colors to a particular range of elevations.

The map is kind of incomplete, without a scale bar, who knows what those colors mean? So let's add a scale bar to this map.

To do this, we'll have to re-make our map. GMT writes its output to PostScript so that you can chain various GMT programs together to draw more things on one map. So to make the above map with a scale bar, do the following (your browser may wrap the lines):

That's right, we run two GMT programs in a row. The first one, grdimage, draws the map and puts the color-coded elevations on the map. The second one, psscale, draws the scale bar on the map. You may notice that there is something different about the way we ran grdimage this time. We added the -K option. PostScript files have a very definite beginning and end, and normally GMT puts them both in. However, specifying the -K tells the grdimage program not to put the end in, so that other GMT programs can write into the file. So after running the above grdimage program, the mola_map_sb.ps file isn't complete, and some viewers may even refuse to display it (rightly so). Running psscale like we do, will put the appropriate end markers in when it is done.

So here are the details on psscale:

psscale

The GMT program psscale draws scale bars on GMT maps.

-D3.5i/-0.5i/3i/0.2ih

This specifies location and size of the bar

-Cmola.cpt

The program needs to know what color palette file to make the scale bar from.

-B500/:"Elevation (m)":

This specifies how to label the scale bar

-O

This is kind of the opposite to the -K in the program above.

-V

I always like to see what it has to say.

>> mola_map_sb.ps

This is slightly different than above, we use two > symbols which indicates that we want to append (instead of overwrite) the contents of the file.

If you used the -L option to grd2cpt, you may want to consider using the -E option with psscale.

Using colors to help understand the elevations is good, but adding a little shaded relief based on the elevations would also be helpful. We'll do that by using the -I option on grdimage that takes an intensity file, which we'll have to create with the GMT program grdgradient, like so:

grdgradient mola.grd -Gmola.int -A0/270 -Ne0.6

grdgradient

The GMT program grdgradient can be used to do lots of things, but this usage will create a file with values from -1 to +1 that can be used by grdimage's -I option.

mola.grd -Gmola.int

These are the input and output files for this program. It takes the mola.grd grid file as input and will write its output to the grid file mola.int. This is a netCDF .grd file too, but I use the .int extension to remind me that it contains intensity values.

-A0/270

This determines where the light shines from. You'll note that I have two values there. The second one is easy to explain, it is a light source that comes from directly west, like the afternoon sun. However, MOLA data are poorly sampled in the E-W direction compared to the N-S direction since that is much closer to the direction of the MGS orbital track. So I add a second light that shines from the north to better show the topograhpy. Ultimately, it really depends on your area, and you should play around until you find something that you like.

-Ne0.6

Since grdimage needs these intensities to be between -1 and +1, we need to use this option to normalize the data. There are a few different normalization schemes that you can use (this one uses the Laplace distribution with an amplitude of 0.6). Again, what works best for one area may not work for another. Play with this until you find something that you like.

Okay, now we have an intensity file that we can incorporate into our map. The sequence is the same as before, but this time using the -I option for grdimage:

If you use an intensity file with grdgradient you may want to look into using the -I option with psscale.

These examples should show you how to create a colorized, shaded relief map from MOLA data. Of course, this kind of thing can serve as a basis for doing a whole lot more. You may want to add contours by running grdcontour between grdimage and psscale, or adding markers with psxy. There are lots of possibilities.