Elevation data may be found in different formats. In general, they are called digital elevation maps (DEM), but not all DEM are equal. You can go to one of the several Web sites that offer free downloadable DEMs. For National Elevation Datasets (NED) from the USGS, for instance, click here. The following steps assume that you are trying to plot a US map and that you obtained the elevation dataset as indicated next. Copy-paste the following URL in your internet browser:

and before you hit enter, change the north and west coordinates as needed. These coordinates correspond to the top-left corner of a tile of size equal 1°, both in latitude and longitude (i.e., between longitudes -120° and -119°, and the latitudes 34° and 35°). The example above downloads the file n35w120.zip. This file contains information about the elevation of the chosen tile at a 1 arc-second (about 30 meters). The ZIP file will be of about 50 MB. After unzipping, you should identify the following two files:

Notice that one can also download a higher resolution file corresponding to 1/3 arc-second (about 10 meters) by changing the TYPE to ned3. 1/3 arc-second are about 350 MB. For most purposes, the 1 arc-second files are more than enough.

The next step is to translate the dataset into a format GMT can understand. For this you will need to install GDAL. Assuming you have already installed GDAL, then execute the command:

The output file will be in the format netCDF. You can learn more about netCDF files and libraries here. The advantage is that the netCDF format can be interpreted by GMT.

In case one needs to plot an area larger than 1°-size, then multiple datasets will be needed. For illustration purposes, suppose we downloaded the following files:

After unzipping put all the .flt and .hdr files in one directory and merge them using gdal_merge.py. Again, you will need GDAL so that you can execute the following command.

gdal_merge.py -of netCDF -o myregion.nc   \
    floatn34w120_1.flt floatn35w120_1.flt \
    floatn35w119_1.flt floatn34w119_1.flt \
    floatn35w118_1.flt floatn34w118_1.flt 

In this example we are merging the last six files into a single file called myregion.nc.

It may also occur to you that the resulting file covers much more than the area you really need. In that case you may want to trim the merging by adding some bounding box. One can achieve this by adding the -ul_lr option as follows:

$ gdal_merge.py -of netCDF -o myregion.nc   \
      -ul_lr -119.5 35 -117.5 33.5          \
      floatn34w120_1.flt floatn35w120_1.flt \
      floatn35w119_1.flt floatn34w119_1.flt \
      floatn35w118_1.flt floatn34w118_1.flt 

where the -ul_lr indicate the uppe-left and lower-right corners of the output dataset.

Merging more than two files, however, is likely to result in problems later on when using GMT because (for reasons unknown to us) the merging seems to give bad results after certain number of grid points. The solution for this is to down-sample the output file. If you type: dgalinfo you can obtain the information of a dataset. See for example:

$ gdalinfo floatn35w120_1.flt 
      Driver: EHdr/ESRI .hdr Labelled
      Files: floatn35w120_1.flt
             floatn35w120_1.hdr
      Size is 3612, 3612
      Coordinate System is `'
      Origin = (-120.001666666700004,35.001666666664136)
      Pixel Size = (0.000277777777778,-0.000277777777778)
      Corner Coordinates:
          Upper Left  (-120.0016667,  35.0016667) 
          Lower Left  (-120.0016667,  33.9983333) 
          Upper Right (-118.9983333,  35.0016667) 
          Lower Right (-118.9983333,  33.9983333) 
          Center      (-119.5000000,  34.5000000) 
      Band 1 Block=3612x1 Type=Float32, ColorInterp=Undefined
      NoData Value=-9999 

This indicates that the pixel size of our datasets is (0.000277777777778,-0.000277777777778). Then, to avoid problems, we will use a lower resolution by resetting the pixel size for the merged dataset. We do this by adding the -ps option as follows.

$ gdal_merge.py -of netCDF -o myregion.nc   \
      -ul_lr -119.5 35 -117.5 33.5          \
      -ps 0.001 0.001                       \
      floatn34w120_1.flt floatn35w120_1.flt \
      floatn35w119_1.flt floatn34w119_1.flt \
      floatn35w118_1.flt floatn34w118_1.flt 

The result will be a dataset that has a pixel size equal to 0.001°. Remember that the -ul_lr and -ps options may not always be required. gdal_merge.py has other options one may find useful.

Once you have the appropriate dataset in netCDF format, then you can execute the following commands using GMT to produce the map.

$ grdgradient myregion.nc -Gmyregion.grd -A0/270 -Ne0.6 -V
$ makecpt     -Cgray -T-1.2/0.6/0.001 -V > colorscale.cpt
$ psbasemap   -JM4i -Rd-119.5/-117.5/33.5/35 -Xc -Yc -B30mg30WSne -K -P -V -K > map.ps
$ grdimage    myregion.grd -J -R -Ccolorscale.cpt -P -V -O -K >> map.ps 
$ ps2raster   map.ps -Tef 

Alternatively, one can omit the first line and use the .nc file directly as the input for grdimage.