Slope maps with GMT
Pete Gadomski January 23, 2017 [how-to] #gmt #ogr #gdal #skiingAbout two-thirds of dry slab avalanches occur on slopes between 30° and 45°. A map of slope angles is therefore a useful tool for safe decision making when backcountry skiing.
The excellent CalTopo provides a suite of map building tools, including a slope angle shading map. But me being me, I want to make my own maps.
This walkthrough describes how to make a slope angle shading map using free and open data and Generic Mapping Tools (GMT). You'll need the following softwares to follow along:
- GMT
- GDAL/OGR
The Makefile used in this example, along with a download script to fetch the map data, are in a gist at the bottom of this post.
The goal
Hidden Valley is a abandoned ski area in Rocky Mountain National Park that closed in 1991. It now is a popular winter destination for backcountry skiing and sledding. We want to color-code dangerous slope angles while providing enough additional context to make a useful map.
Getting the data
We'll use data from the USGS 3D Elevation Program to produce our slope, hillshade, and contour grids. Our stream and lake overlays will come from the National Hydrography Dataset, and the roads from the National Transporation Dataset. All these products are browsable and downloadable through The National Map, but here's direct links to the datasets required for this exercise:
- n40w106 DEM as GridFloat
- n41w106 DEM as GridFloat
- Subbasin 10190005 as a shapefile
- Colorado roads as a shapefile
All told, it's about 4.6G of data, organized like this on my system:
$ tree -P '*.shp|*.flt'
.
├── dem
│ ├── n40w106
│ │ └── n40w106.shp
│ ├── n41w106
│ │ └── n41w106.shp
│ ├── usgs_ned_13_n40w106_gridfloat.flt
│ └── usgs_ned_13_n41w106_gridfloat.flt
├── roads
│ └── Shape
│ ├── Trans_AirportPoint.shp
│ ├── Trans_AirportRunway.shp
│ ├── Trans_RailFeature.shp
│ ├── Trans_RoadSegment.shp
│ ├── Trans_RoadSegment2.shp
│ ├── Trans_RoadSegment3.shp
│ └── Trans_TrailSegment.shp
└── water
└── Shape
├── NHDArea.shp
├── NHDFlowline.shp
├── NHDLine.shp
├── NHDPoint.shp
├── NHDPointEventFC.shp
├── NHDWaterbody.shp
├── WBDHU10.shp
├── WBDHU12.shp
├── WBDHU2.shp
├── WBDHU4.shp
├── WBDHU8.shp
└── WBDLine.shp
7 directories, 23 files
Creating a simple raster map with GMT
To prove that our system is working, let's create a rainbow elevation map from our DEM data in our area of interest.
We'll use a Makefile and place our generated products in build/.
First, we define our area of interest and create a rule for our build directory:
XMIN:=
XMAX:=
YMIN:=
YMAX:=
:
Next, we subset our DEM data to our area of interest using a VRT. For some reason GMT (on my system) isn't happy reading VRTs, even though it should be able to use any GDAL-y data source, so we use the VRT to create a GeoTIFF:
: : Finally, we create our simple rainbow elevation map.
I do this in two steps — first I create the PostScript file, then I use psconvert to turn it into a png:
: : A few things to note:
grdimagetakes many more options, but for now we're keeping it simple.- The
psconvertcommand converts the PostScript file into a png with transparency (-TG), cropped (-A), and rotated into portrait mode (-P).
We make our map by running make build/hidden-valley.png:
Creating more derivative products
We've already created one derivative product, our cropped DEM build/dem.tif.
Let's create the rest of the products we need to make our map.
I'll step through the rules one-by-one and talk through what they're doing.
Calculate the slope
: grdgradient does derivative math on gridded data.
In this case, we use the combination of the -D and -S options to write the magnitude of the gradient vector to a file, build/slope.nc.
We also need to specify -fg because our data is a geographic (lat/lon) grid and we need to convert those latitudes and longitudes to meters before calculating the gradient.
Once we've calculated the raw gradient magnitude, we convert that scalar value to an angle, in degrees, using grdmath.
Create a hillshade
: In this case, we use grdgradient to create a hillshade.
We use -N to control the intensity of the hillshade — since we'll be drawing contours and additional overlays, we dial down the intensity to 0.6 so the map is a bit less noisy.
Crop and convert the vector data
OGR2GMT:=
: : : : Creating flowline (streams) and waterbody (lakes, etc) vectors is a simple matter of cropping the source data with ogr2ogr and writing it out in the GMT format.
Creating the roads file is a bit more tricky, since the road data comes in three seperate shapefiles.
I wasn't able to figure out how to -addfields to a GMT vector file, so I create an intermediate build/roads.shp, then convert that file to gmt.
Create a color palette
: This rule creates a color palette that we'll use for our slopes.
Putting it all together
Now that we've created all of the intermediate products, let's rewrite our PostScript rule to create our final map. I'll step through the rule components one-by-one, with an explanation immediately after.
: Add all of our derivative products to ensure they get built.
grdimage build/slope.nc -Cbuild/slope.cpt -Ibuild/gradient.nc -Ba -B+t"Hidden Valley" -JM8i -Rbuild/slope.nc -K > $@
Draw the slope grid, colored with our custom palette.
The -JM8i specifies the projection for these data, in this case a Mercator projection that is 8 inches wide.
The -R instructs grdimage to use the bounds of the build/slope.nc grid as the bounds of the plot.
The -I option adds the hillshade, the -Ba adds a border, the -B+t titles the plot.
The -K option "keeps the output file open" so we can add more layers.
grdcontour build/dem.tif -C20 -A100 -J -R -K -O >> $@
Add contours spaced 20 meters apart, labelled every 100 meters.
The -J and -R options inherit the values used in the previous command, so we don't need to specify them again.
The -O option instructs grdcontour to output "overlay" PostScript that is meant to be appended onto PostScript data "kept open" with -K.
psxy build/flowline.gmt -W0.6p,'#377eb8' -J -R -K -O >> $@
psxy build/waterbody.gmt -G'#377eb8' -J -R -K -O >> $@
psxy build/roads.gmt -W1p,'#f781bf' -J -R -K -O >> $@
Add the vector data.
The -W option specifies the pen used to draw lines, and the -G option specifies the fill.
psbasemap -Lx5.2i/-0.7i+c$(YMIN)+w1k -J -R -K -O >> $@
Add the length scale, set to 1km.
psscale -D0i/-0.7i+w3i/0.2i+h -Cbuild/slope.cpt -G0/60 -By+l"Slope angle" -I -O >> $@
Add the color legend.
Now you can run make build/hidden-valley.png and get your map!
Conclusion
By automating this map-generation, you can easily create maps for new areas. I've souped up this example into this project, which I can use to quickly create maps of new areas I'm exploring.
Here's the complete Makefile used for this example, along with a download script to grab the data you need: