10th Asian Regional Conference of IAEG (2015)

Geomorphological Analysis Tool based on the Generic Mapping Tools
Hisao SUNOUCHI(1) and Yuichiro MIYATA(2)(1) Kouei Engineering Consultant Co. Ltd., Japan
E-mail:sunouchi@koueicon.co.jp(2) Graduate School of Science and Engineering, Yamaguchi University, Japan
Abstract
This paper describes a new geomorphological analysis tool developed using the Generic
Mapping Tools (GMT), a freeware mapping tool developed at Hawaii University. GMT is commonlyused in the earth sciences, but less so by engineering geologists unfamiliar with computer
programming. The new geomorphological analysis tool is based on GMT commands and is executedat the Windows command prompt. The tool is based on numerous complex parameters that can beentered into the calculations using a Microsoft Excel worksheet and without the need for any
knowledge of computer programming. Various types of maps and diagrams can be generated with thetool, including maps representing contoured data, color gradations, topographic maps, cross-sections,slope maps, relief energy maps, drainage system maps, summit level maps, valley level maps, andstereographic projections of surface orientations.

Keywords: GMT, geomorphological analysis, DEM
1.

Introduction

onwards).
GMT is used mainly in the research field. As
GMT is a command-based software, it can bedifficult to use for those unfamiliar with computer
programming. Therefore, we have generated a seriesof command lines as an Excel spreadsheet (Fig. 1).
We describe here a new geomorphologicalanalysis tool developed using GMT, with DEMs of
10 mesh and river data published by the Geospatial
Information Authority of Japan (GSI). First, we
describe the creation of some basic maps, and then of
specific maps that require complicated gridcomputing.

In recent years, digital elevation model (DEM)
data, together with aerial laser surveys, have been
widely used in geomorphological analysis. A varietyof geographic information systems (GIS) software isavailable to analyze such data. In this paper, we
describe an approach to topographical analysis basedon The Generic Mapping Tools (GMT), open source
mapping software developed by Wessel and Smith(1991). The software has no GUI (graphical userinterface), and is operated using a series ofcommand-line tools, or via a shell script.
The advantages of GMT are that it: (1) has strong
grid manipulation functions; (2) is free; (3) can
output to PostScript; and (4) runs in both theWindows and UNIX environments. The software isregularly upgraded and is stable, even on machines
running Japanese versions of Windows (from GMT5
File namelongitude
latitude
Drainage map

2. Examples of Basic Maps
2.1 Shaded topographic map
A shaded topographic map (Fig. 2) is created byapplying a color corresponding to the altitude of each

↓Check the increment by grdinfo command.
ytj10mytj10m.grd
gridincrement dx
5.38E-05 degree =5.99 mgridincrement dy
4.51E-05 degree =133.3058611
133.3581389
W:H=
4.3
5.8 km33.46949444
33.51600278
max elevation
1000 mgrdmath ytj10m.grd median.grd SUB = tani.grd
133.332/33.49274861/ -R133.3058611/133.3581389/33.46949444/33.51600278↓Copy and paste into the command window and execute．
ytj10mFilter size (m)50 median filter
grdfilter ytj10m.grd -Gmedian.grd -D2 -Fm0.05 -Vkm
0.05 abstract valleys
grdmath ytj10m.grd median.grd SUB -0.025 GT = tani.grd
color palletmakecpt -Cdrainage.cpt -T0/0.1/0.1 -Z > temp.cpt
project grid image

5.02 mgrdimage tani.grd -R133.3058611/133.3581389/33.46949444/33.51600278 -Jt133.332/33.49274861/306c-Ctemp.cpt -K -V > tani_50.eps

contour lines
grdcontour ytj10m.grd -R -Jt -L0/1000 -C50 -O -W0.05,red -K -V >> tani_50.eps
rivers from library psxy r_Shikoku.xy -R -Jt -Wblue -O -K -V >> tani_50.eps
scale barpsbasemap -R -Jt -Lfx10/1c/33.46949444/1k+l -B1mWeSn -O -K -V >> tani_50.eps
title
echo Map of Drainage Net (filter size = 50m) | pstext -R -Jt -F+cTL -Gwhite -O -V >> tani_50.eps
jpg file
ps2raster tani_50.eps -GC:\programs\gs\gs9.10\bin\gswin64c.exe -Tj -A -P -V

Fig. 1 Example of a set of GMT commands in an Excel worksheet.

1

10th Asian Regional Conference of IAEG (2015)

grid point, and assigning an intensity corresponding
to the slope direction relative to the light source using
the “grdimage” command. The color is taken fromthe GMT internal color palette, and the light source in
GMT can be set up in two directions of which thegreater directional derivative is used for shading. The
light sources in Fig.. 2 were set in the west (270°) andnorth (0°), meaning the north–south and the east–west aspects of the terrain are emphasized. The map
shows that mountain ranges, valleys, and narrowplains are oriented ENE–WSW or E–W (Fig. 2),
although some valleys are oriented NNW–SSE andN–S.

command. Additional color palettes can be obtained
from http://soliton.vm.bytemark.co.uk/pub/cpt-city/views/totp-cpt.html. Fig. 4 was generated using thecolor palette referred to by Young (1969).
In the southern part of Fig. 4, the slope gradient is
steeper on the southern slopes (>30°) than on thenorthern slopes, where the gradient is commonly<18°. The distribution of flat ground (slope < 5°)forms like a grid pattern.

degGradient

A

A’

Fig. 4 Slope map assigning colors to the slope
angle at each grid point.
2.4 Relief energy map
The relief energy map is compiled by calculating
the maximum and minimum elevation within a circle
of diameter 1 km centered on each grid point using
the “grdfilter” command; the difference between thetwo elevations is mapped using a color scale. It isalso possible to use a rectangular, rather than circular,
window in the command. The boundary between
areas with relief energy of <200 m, and mountainous
areas with relief energy of >200 m, is well defined,showing that topographical classification by relief
energy is possible in Fig. 5.

Fig. 2 Shaded topographic map. The line A–A’
indicates the section line shown in Fig. 3.
2.2 Cross-section
A cross-section along line A–A’ in Fig. 2 isdrawn in Excel using the x–y dataset, which is
obtained by projecting the elevation of the grid points
that are within 10 m of the specified section line (Fig.
3). The command used was “project”.
The cross-section shows that the southern slopesof the mountains at the southern edge of the basin are
steeper than the northern slopes. The plain dips
gently to the northwest.
A

A’

Fig. 3 Cross-section along line A–A’ in Fig. 2.
Elevation data are from within 10 m of thesection line, which is drawn 5 km to thenorthwest and southeast from the center ofthe map.

Fig. 5 Relief energy map based on the
difference in maximum and minimum
elevations within a 1 km circle at each grid
point.

2.3 Slope map
A slope map is created by assigning colors to theslope angle at each grid point using the “grdimage”
command; the slope angle is obtained using the“grdgradient” command (Fig. 4).
In addition to the 23 built-in color palette tables,
users can create their own tables using the “makecpt”

2.5 Aspect map
The aspect map is applied a color according to thedirection of the slope inclination measured clockwisefrom north using the “grdgradient” command. This2

10th Asian Regional Conference of IAEG (2015)

map highlights the orientations of ridges and valleys,
and curved and parallel ridges are clearly seen in Fig.
6.

Gradient

deg

Fig. 8 Drainage slope map. The
relationship between channel orderand gradient can be seen.

Fig. 6 Aspect map.
The direction of the slope inclination ismeasured clockwise from north

3.2 Summit level map
We describe the summit level maps generatedusing the valley filling method and the grid method.
(1) Valley filling method
A summit level map (Fig. 9) can be generated byfilling the valleys, which leaves the ridges
highlighted. In this method, the ridges are extractedusing a median filter as follows.
(a) Subtract the elevation of each point from themedian.
(b) Set the points where the elevation differenceis greater than the value set in advance (2 mhere) to 1, and the others to 0, to leave only
prominent ridges.
(c) Taking the product of this and the DEM
elevation, elevation becomes zero except forridges.
(d) Make a data grid using the “ridge.grd”
command by setting those points with anelevation of zero as no data (NaN).
This procedure generates the following GMT
commands:

3. Examples of Applied Maps
3.1 Drainage system map
The extraction of drainage systems is possible byidentifying pixels that are lower than a target pixel
based on elevation differences. The conditions used
for the extraction can be varied by obtaining theneighboring "average" elevation as well as bychanging the diameter of the window. GMT allows
the use of a simple average, or the mode and median,creating a weighted normal distribution using the“grdfilter” command. In general, the grid data uses
the elevation values of a square, such as 3 × 3 or 5 ×
5, and it is possible to specify the diameter (m) of thesearch window (filter) in GMT.
Fig. 7 is created by extracting the points forwhich the difference between the center and the
median elevation of the window was negative. The
drainage slope map is shown in Fig. 8, which clearlyshows spatial differences in the density, the pattern ofdrainage system, and the relationship between
channel order and gradient.

grdmath median.grd dem.grd SUB 2 GT dem.grd
MUL 0 NAN = ridge.grd(1)Here, median.grd and dem.grd are the grid data
files for the median and DEM, respectively. The
calculations taking the difference (SUB) or product(MUL) are in the “grdmath” command. Thus,
implementing these GMT commands uses a reverse
Polish notation.
(e) Finally, interpolate between the areas with nodata, such as around valleys. There are twointerpolation methods available in GMT,
which use either the average values of theperipheral or the smoothing function. We used
the latter method, and the “surface” command.
Three parameters must be specified in the above

Fig. 7 Drainage system map.
Spatial differences in the density
can be seen.

3

10th Asian Regional Conference of IAEG (2015)

procedure: the median filter diameter, the altitudedifference from (b), and the parameters from theinterpolation function in (e).
The river system is also shown in Fig. 9. In places,
the course of the river cuts the slope of the rising
gradient (uphill), and river capture or stream
movement may occur at these points.

3.3 Valley level map
The valley level map (Fig. 11) is generated byremoving the ridges. Fig.. 11 shows a map that iscreated using the grid method to interpolate between
the lowest elevation points in each grid. Almost all of
the ridges have disappeared in the southeasternsection, which is in contrast to the other areas. The
valley level map is useful to infer groundwater
system.

Fig. 9 Shaded summit level map generated using thevalley filling method. The map is created by fillingvalleys between high ridges more than 2 m higher
than the median elevation. The diameter of the filter
is 100 m.

Fig. 11 Valley level map created using thegrid method on a 500 m grid.
3.4 Extraction of nick line and ridge shape
There is little difference between the medianelevation and nick point elevation when the filter isapplied to a single slope (Fig. 12). The elevation
frequency distribution in the window, in the case of
the slope, does not become a unimodal distribution,such as the normal distribution, and this also differsfrom the unbalanced single-peak distribution. This
point is the cause of the difference between the mean
and the median.
If the slope is convex upward (Fig. 12, left), theelevation of the window center is higher than the
average, but the center point is lower than the averagein the case of a convex downward slope (Fig. 12,
right). However, the difference between the medianvalue and the window center point elevation of theslope is not evident in the figure.

(2) The grid method
Summit level maps (Fig. 10) based on the grid
method are created by interpolating between thehighest points in the rectangle, separated at regularintervals. The procedure is as follows.
(a) Cut out the predetermined rectangle using the“grdcut“ command and obtain the coordinatesand elevation of the highest point in therectangle with the “grdinfo” command.
(b) Move the partition to the next rectangle andrepeat the procedure in (a) in a whole area.
(c) Interpolate between the points obtained fromthe above procedures.
As the process in (b) is repeated within apredetermined area, it can be performed using an
Excel macro.

Fig. 12 Frequency distribution of elevation within
a window on the slope. The elevation andhorizontal distance in the window are indicatedby z and x, respectively. The average is shown as
the mean, and the median as the cumulative
frequency of 50% of the multiple elevation datain the window.

Fig. 10 Shaded summit level map created using thegrid method on a 2 km grid.

Ifelevation
increases(ordecreases)
monotonically with horizontal distance in a4

10th Asian Regional Conference of IAEG (2015)

two-dimensional cross-section, regardless of the
shape of the cross-section, the elevation of the center
point equals the median. This may be explained asfollows. As shown in Fig. 13 (left), window width in
the two-dimensional cross-section is the distance
between the equidistant points from the center. The
cumulative value of the altitude of all points in thewindow equals to the cumulative value of each
altitude from the lowest to the highest elevation(Fig.13, right). When considered as a continuous
distribution with a sufficient number of points, thecumulative number of points is equal to thehorizontal distance.

Therefore, the median value of the altitude
distribution is the altitude where the cumulative
number of the points is 50% of the frequency
distribution, equals to the distance is 50% from thewindow edge. This point is the center point of thewindow. In other words, the elevation of the center
point of the window is at the median altitude. This
holds regardless of the shape of the slope. The slope
may be any shape if the altitude increases (ordecreases) monotonically; moreover, the relationshipof the median and the window center is independentof the gradient of the slope.
The median of the elevation distribution on any
slope is equal to the elevation of the window center
point. Furthermore, because the difference between
the median elevation of surrounding points and theelevation of each point is always zero, in principle anick point cannot be extracted by difference.
However, in practice, because the slope including they-direction forms a curved surface in thethree-dimensional terrain, a slight difference does
occur.
On the other hand, when a window is large, bothridge and valley may be included in the window. The
median value is lower than the elevation of thewindow center on the ridge but higher in the valley.

Fig. 13 Topographic section (left) and cumulative
curve of elevation (right). The slope indicates thatthe median elevation is equal to the window center
elevation (see text for details).
a

b

c

d

Fig. 14 Maps showing differences between the elevation of the window center and themean elevation (a, b), and between the median and mean elevations (c, d). The diameters of
the windows are 50 (a, c) and 400 m (b, d).

5

10th Asian Regional Conference of IAEG (2015)

In this case (i.e., a large window), the relationships
are almost the same with the differences calculatedusing the simple average and the window center point.
By using the characteristics of both, it is possible toextract only nick lines.
As the elevation difference between the averageand the window center point (Fig. 14a, b) is simply
used to indicate whether higher or lower than theambient points, ridges and valleys are highlighted aswell as undulations on the slopes. In contrast,
elevation differences between the window center
point and the median tend to extract only ridges andvalleys without showing the unevenness of the slopes.
Therefore, taking the difference between the average
elevation and the median (Fig. 14c, d), a number ofridges and valleys are excluded and this results inenhancement of the undulations on the slope.
By exploiting these differences, features such as
nick lines and failure terrains are easier to extract.
Furthermore, as shown in Fig. 15, the relationshipbetween the average and median differs depending on
the shapes of the ridges and valleys. This relationshipcan also be helpful in the evaluation of shape, andcan be used by changing the window size to match tothe scale of the topography of interest.

abFig. 16 Equal area projection of poles to theground surface (a) and contour diagram of thepole data (b) plotted in the lower hemisphere.
The contour diagram was made using Open
Stereo (Grohman et al., 2010). The sample
size is 22,650 data points.
4. Conclusions
GMT offers flexible and powerful grid computing
functionsandcould
becomeausefulgeomorpholo-gical analysis tool. In addition to thebasic drawings typically used in geomorphologicalanalysis, it is also possible to create various otherfigures for specific purposes such as elucidation ofriver capture and evaluation of ridges’ shapes. To
allow users unfamiliar with computer programing tocreate command sets, we developed a new tool that is
executable at the Windows command prompt. This
tool is based on numerous complex parameters thatare entered into the calculations using a Microsoft
Excel worksheet, and requires no previous experiencein computer programming. As GMT has no graphicaluser interface, a new application incorporating a GUI
would be desirable in the future.

Fig. 15 Comparison of differences between the
median and mean along various types ofridges and valleys.
3.5 Stereo projections
A number of projections are available in GMT,
and it is possible to create a stereo projection byplotting the dip and strike of the ground surface in a
Lambert equal-area projection (Fig. 16a). Comparing
stereo projections of discontinuities, such as bedding
lanes and joints in the terrain, helps to highlight therelationship between geological structure andtopography. As the dip and strike data can be stored
in a text file, it is also possible to export the data toother programs for constructing stereo projections(e.g., Grohman et. al., 2010, Fig. 16b).

References
Grohman, C. H., Ginald, A. C., and Campanha, G. A.
C.
(2010):
OpenStereo,
Open
Source,Cross-Platform Software for Structural Geology
Analysis, AGU Fall Meeting.
Wessel, P. and Smith, W. H. F. (1991): Free software
helps map and display data, EOS Trans. AGU, 72,
pp. 441.
Young, A. (1969): Slopes, Oliver & Boyd, pp. 288.

6