Plotting USGS 3d Elevation Data in GMT

This post is about plotting high resolution digital elevation model (DEM) data in GMT for high-quality maps. Every time I need to do this, I end up googling around, wasting a bunch of time, so this post mostly serves as a reference for me.

First things first, the USGS provides elevation data at numerous resolutions, the most common being 1 arcsecond (roughly 30 meters) and 1/3 arcsecond (roughly 10 meters). Other resolutions exist, but these are the two you will probably use most.

Getting the elevation data can be done through the USGS's National Map web service (https://viewer.nationalmap.gov/basic/#productGroupSearch). The process is rather simple, just zoom into the area you care about (you can also enter coordinates) and then under "Elevation Products (3DEP)", click which one you want. After clicking 'Find Products', you can go to the Products tab and toggle on and off the footprints of the different maps. They are done in 1 degree squares, so you will need to add several to your cart if its a large area. Once you have added the different images to your cart, you can download them in geoTIFF format. For the 1/3 arcsec files, they can get pretty large, so think about the spatial extent you are working with. In the figure below, you can see the 4 squares that I have selected at 1 arcsec resolution.

After downloading, I need to do some manipulating of the data. The geoTIFFs need to first be converted to netCDF format, which is done rather simply through

gmt grdconvert USGS_1_n38w118.tif 1.grd

In the above command, the first file is the input, the second is the output. Do this for the 4 files. For the next step, we want to combine the 4 .grd files into a single file. Unfortunately, the edges wont line up perfectly, so we need to resample the data. If you need to know the rough extent of the .grd file, enter 

gmt grdinfo 1.grd

To resample, we use grdsample. For the 1.grd file, the command looks like

gmt grdsample 1.grd -R-118/-117/37/38 -G1n.grd

where the -G flag is the output file. Once we have our 4 resampled files, we can combine them with grdpaste, which can only handle two files at a time. To combine our 4 resampled files, we do the following

gmt grdpaste 1n.grd 2n.grd -G12.grd

gmt grdpaste 3n.grd 4n.grd -G34.grd

gmt grdpaste 12.grd 34.grd -G1234.grd

From here, we could just simply plot the file 1234.grd using grdimage, but it will look much nicer if we create a hillslope gradient file. For this, we use the grdgradient command 

gmt grdgradient 1234.grd -G1234.grad -A0/270 -fg

Where the -A flag controls the direction of the illumination and the -fg flag uses a flat earth approximation. From here, to plot the data on a simple map, we would do the following

gmt makecpt -Cgray -T-4000/4000/10 -Z > colors.cpt

gmt psbasemap -R-119/-117.0/37.5/39 -JM20 -Bxafg -Byafg -BWeSn -K -X-14i -W > $NAME.eps

gmt grdimage 1234.grd -J -R -O -K -I1234.grad -Ccolors.cpt >> $NAME.eps

pscoast -R-119/-117.0/37.5/39 -J -O -K -Df -Na -W -Sblue >> $NAME.eps

In the 4 commands above, I created a greyscale colormap between -4000 and 4000 meters at 10 meter interval, created the basemap, plotted the .grd file with the gradient hillshade and then finally plotted the political boundaries and colored the lakes blue. The finished product can be seen below. From here, I would maybe play with the illumination direction or the colormap, but its publishable as is.

 

Below are the contents of the script I wrote to do these 4 geoTIFFs

#!/bin/bash
export PATH=/usr/local/Cellar/gmt/5.4.4_2/bin:$PATH

grdconvert USGS_1_n38w118.tif 1.grd
grdconvert USGS_1_n38w119.tif 2.grd
grdconvert USGS_1_n39w118.tif 3.grd
grdconvert USGS_1_n39w119.tif 4.grd

grdsample 1.grd -R-118/-117/37/38 -G1n.grd
grdsample 2.grd -R-119/-118/37/38 -G2n.grd
grdsample 3.grd -R-118/-117/38/39 -G3n.grd
grdsample 4.grd -R-119/-118/38/39 -G4n.grd

grdpaste 1n.grd 2n.grd -G12.grd
grdpaste 3n.grd 4n.grd -G34.grd
grdpaste 12.grd 34.grd -G1234.grd

grdgradient 1234.grd -G1234.grad -A0/270 -fg
exit 0