How to plot a shaded-relief map
These are a series of steps you can follow in order to create a shaded-relief map for publications.
Obtain elevation data.
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:http://gisdata.usgs.gov/TDDS/DownloadFile.php?TYPE=ned1f_zip&ORIG=TNM&FNAME=n35w120.zip
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:floatn35w120_1.flt floatn35w120_1.hdr
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.
Translate the dataset to a friendly format
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:gdal_translate -of netCDF floatn35w120_1.flt n35w120.nc
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.
Merge multiple datasets (if needed).
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:n34w118.zip n34w119.zip n34w120.zip n35w118.zip n35w119.zip n35w120.zip
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.
Plot using GMT.
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.