Adding real data into GMT map, making a hillshade, and positioning

Alright, let’s get down to business plotting some real data on a map. Grab the zip file here, and uncompress it in the folder you want to build a map in. The zip contains the data to map, as well as  on of my favorite color palettes from cptcity by M. Burak Yilkilmaz (file named mby.cpt).

Here’s the end product of this tutorial, we’re going to make two maps, the first will be a simple data map, and the second a hillshaded version of the same map. We’ll also learn how to position them side by side like below.


We’ll start by plotting a simple map, with just the data represented as different colors (left side of the map). The data we’re plotting are from a digital elevation model (DEM), and of Eastern North America. Making this image is simple, we simply need to call
gmt grdimage -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Cmby.cpt -V >
where you already know most of the switches above. grdimage is a gmt tool that you can use to plot any data stored in the netCDF format. To color the image based on the values in the .nc file, we need to use a color palette. The color palette is called with the -C switch, followed by the name of the color palette. More on color palettes in a future tutorial though.

Now, let’s talk about hillshades. A hillshade (also known as shaded relief) is an element of a map that adds depth to topography. It is essentially simulating a light source and the shadows cast on the landscape from the hills onto the valleys. Making one in GMT is straightforward but requires a few steps. The main step of the process is done first with grdgradient. We use the following command
gmt grdgradient -A345 -Ne0.6 -V
In this command, -G has a different meaning than in this tutorial; another common use of the -G switch is to give the name of the output file for the command. -A is the azimuth from which the light source should be directed, and -Ne is a normalization parameter for the output from the calculation, just use 0.6.

Then, we will ensure that the shaded relief map is reasonably balanced by normalizing it again along a Gaussian distribution
gmt grdhisteq -N -V

Finally, in order to ensure the values are between -1 and 1, a final step required for creating the intensity file needed for a hillshade to the image. We do this by dividing by some value just larger than the range of values in the histogram equalized grid file. Let’s find this range by looking at the info for this file.

$ gmt grdinfo Title: Produced by grdhisteq Command: grdhisteq -N -V Remark: Normalized directional derivative(s) Gridline node registration used [Geographic grid] Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5 x_min: 278.004166667 x_max: 286.004166667 x_inc: 0.00833333333333 name: longitude [degrees_east] nx: 961 y_min: 35.9958333333 y_max: 41.9958333333 y_inc: 0.00833333333333 name: latitude [degrees_north] ny: 721 z_min: -4.67873573303 z_max: 4.67873573303 name: Elevation relative to sea level [m] scale_factor: 1 add_offset: 0 format: netCDF-4 chunk_size: 138,145 shuffle: on deflation_level: 3

So we will use 5 as our divisor, in the command
gmt grdmath 5 DIV =
which takes our -hist file and divides every grid cell by 5, putting it within the range of -1 to 1.

Finally, we use
gmt grdimage -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Cmby.cpt -V >
the grdimage tool again, in almost exactly the same way as before, except adding the -I switch with our normalized intensity file as the input.

Now, to put the two images side by side, we need to set different -X and -Y targets for each image. These switches set the offset of the lower-left corner of the PostScript layer being produced, from the last reference point used on the PostScript page. The a argument can be added to the switches to reset the reference point to the lower-left corner of the page. Try using -Xa1 -Yc for the flat image and then overlay the hillshade image at -Xa6.5 -Yc. Remember, you’ll need to add in -K, -O, and >> in order to produce a multi-layer map.

See if you can make it work by yourself, but if you need a hint, you can get the script I used from here.

6 thoughts on “Adding real data into GMT map, making a hillshade, and positioning”

  1. Hi Andrew,

    First, thanks for this great post. I’m a GMT newbie struggling with the learning curve, and you do a great job of explaining how to do things.

    I would like to use similar hillshade maps for a project I am working on in the Pacific Northwest. I’d like to use NASA’s SRTM data but can’t figure out how to convert from the raw data (e.g. N48W121.hgt) to a netCDF format. Wondering if you have any suggestions? Alternatively, do you know of a source where I could get netCDF formatted files like the one you use in the post?


    1. Hello, thanks for the kind words. You can do the translation to netCDF with a GMT tool called xyz2grd, but it will be much easier done with GDAL. If you do not already have GDAL tools installed see here for instructions: []. Using gdal_translate, you simply point to the file and specify you want a netCDF file output, e.g., gdal_translate -of netCDF N48W121.hgt I think gdal_translate should also work if pointed directly to the file, although I’ve never tried. Alternatively, you can get GeoTIFF versions of all the SRTM data from the USGS’s EarthExplorer [] and can then convert these using gdal_translate (this is the route I normally take). If you are operating on a number of files you will want to do a batch conversion which can be easily done with gdal_translate.

      1. Andrew, thank you so much for the quick and thorough reply. I downloaded GDAL and used it to batch process my files to netCDF format. Worked like a charm!

  2. Hi, Andrew! Is there a way to modify the mby.cpt, so that the boundary of land color and the sea color can be at the middle? I am always wondering if it is possible to adjust the cpt with different max. and min. altitudes, but the boundary of land color and sea color always corresponds to zero altitude (sea level).

    1. Hi Yuan-Kai, yes, you can modify the mby.cpt file directly. You need to adjust the extent covered by the positive and negative ranges of elevations top be equal on either side of 0. You can also check out makecpt in GMT to create a color palette by modifying an existing one! Good luck!

      1. Thank you for te quick response! Modifying all the RGB values in the entire mby.cpt file is I guess too hard and also too much work. I will first check out “makecpt” in GMT to see if I can create a palette with equal positive and negative ends. Thanks!

Leave a Reply

Your email address will not be published. Required fields are marked *