2.

GIS – Techniques and Methods

This chapter serves as an introduction and background to the use of Geographical

Information Systems or GIS as a tool in geomorphic analysis. It attempts to describe why
GIS was used over manual cartographic methods and the techniques and processes behind it.

2.1 Introduction to Geographical Information Systems

A Geographical Information System or GIS is a computer system used for managing spatial
data (Bonham-Carter, 1994). In such a system all data have locations that are either known
or can be calculated in the form of geographical coordinates. It will consist of a group of
computer programs or modules that are linked through a user interface, from which access is
gained to the systems functions. These modules will provide the capacity for data capture,
input, manipulation, transformation, combination, analysis, modelling and output (Bonham-
Carter, 1994).

GIS also provide a mechanism capable of integrating databases and spatial datasets and
provide a primary tool for manipulating, storing and accessing this data (Moore et al. 1991).
They have a wide variety of applications and are heavily used in a wide-range of scenarios
from town planning applications to geomorphic analysis and integration with landscape
modelling. Advantages of GIS are that they allow many different applications within a
single system, without the need to transfer large data files between different programs.
Moore et al. (1991), consider their advantages to be as follows:

- ability to integrate spatial and other information within a single system

- offer a consistent framework for analysing the spatial variation across landscapes

- allows geographical knowledge to be manipulated and displayed in a variety of forms,

including maps

- allows connections to be made between entities based on geographical proximity and

characteristics that are vital to understanding and managing activities and resources.

Considerable research and development has been undertaken in order to create complex
integrated GIS which allow large spatial and sometimes temporal datasets to be handled.
One of the reasons that GIS technologies can be applied so widely is because they can

Chapter 2 – GIS – Techniques and Methods 27

display spatial information in integrative ways that are readily comprehensible and visual
(Moore et al., 1991).

The current application of GIS to geomorphic analysis allows characteristics to be rapidly

extracted from a representation of a real landscape. Previous techniques involving
measurements from extensive field surveys and paper map resulted in poor reproducibility in
quantitative studies and difficulties in comparisons over large areas.

GIS store data in a variety of ways. The two main forms of data storage are in coverages or
in grids. Coverages store vector data which include lines or arcs, points and polygons. They
can be used, for example, for basin outlines (polygons), contours (arcs) or cities or towns
(points). Grids store data in a matrix consisting of a regular tesselation of cells (raster data),
such as a series or string of elevation values.

Any digital representation of the continuous variation of relief over space is known as a
Digital Elevation Model or DEM (Burrough, 1986). The term Digital Terrain Model or
DTM is often used, but terrain often implies some other attribute of the landscape, such as
slope or aspect in addition to only elevation data (Wood, 1996a). So for solely elevation
data, the term DEM is used (Burrough, 1986). DEMs permit visualisation of topography and
a more detailed analysis of topographic data structure than is possible from conventional
paper topographic contour maps (Mayer, 2000). Digital elevation data can be stored in a
number of different ways (Walker and Willgoose, 1999):

1) A two-dimensional array of elevations that represents elevations on a regular grid –

known as GRID DEMs

2) A set of irregularly spaced xyz coordinates, interpolated to give a series of delaunay

triangles – known as a Triangular Irregular Network or TIN

3) A series of contour strings – Vector DEMs (Moore and Grayson, 1991).

The most common type in use for geomorphic analysis are those stored in regular matrices –
GRID DEMs. These have many advantages over other forms of digital elevation data due to

Chapter 2 – GIS – Techniques and Methods 28

their ease of computer implementation and computational efficiency. They are not however
without their disadvantages (Moore et al., 1991; Burrough, 1986):

1) They cannot easily handle abrupt changes in elevation

2) The size of the grid mesh affects the results obtained and the computational efficiency

3) The computed upslope flow paths used in hydrological analysis tend to zigzag which is
relatively unrealistic

4) Precision is lacking in the definition of specific catchment areas

5) Large amounts of data redundancy in areas of uniform terrain

6) The inability to adapt to areas of differing relief complexity without changing the grid
size

7) The exaggerated emphasis placed upon the grid axes for certain kinds of computation
such as line of sight calculations.

Although there are many disadvantages with regular grid based systems, many of these are
also inherent to other forms of digital landscape definition such as TINs or vector based
contour systems. No system is ideal for analysis and a compromise must always be made.
The relative advantages in terms of simplicity and computational ease make GRID DEMs
the basic tools for geomorphic analysis.

[Link] Data acquisition

Digital elevation data can be acquired in one of two ways. The first and simplest is to
digitise a paper topographic contour map of the region of interest. This tends to be a very
laborious procedure involving manual digitising of all the contour lines on a map, especially
if high resolution datasets are required. Following collection the data is then either stored as
contour strings or is interpolated between the points of known elevation in order to create a
regular grid. DEMs created in this way are described as Cartometric DEMs (Walker and
Willgoose, 1999). The other form of data acquisition involves remote sensing procedures
and the use of stereoscopic images. These are created using photogrammetric methods from

Chapter 2 – GIS – Techniques and Methods 29

aerial photographs and digital photogrammetry using stereo satellite images such as SPOT
data (Mayer, 2000). Other remote sensing techniques include laser altimetry, synthetic
aperture radar (SAR) and radar interferometry. These datasets have the disadvantage that
elevations often represent the tops of trees in dense forested areas and building tops in urban
areas and consequently require careful editing before use. DEMs created by these methods
are described as Photogrammetric DEMs (Walker and Willgoose, 1999). This data, like the
digitised contour data is interpolated onto a regular grid ready for use. Both techniques have
their limitations and are often used in conjunction to generate DEMs. It is important to note
that a DEM can be no more accurate than the dataset from which it was created (Carter,
1988).

After creation of the DEM it is fair to assume that it provides an accurate representation of
the real world land surface, but the level of accuracy is very difficult to assess as we rarely
have an independent model of the real world (Carter, 1988). Comparison with real elevation
data i.e. topographic maps provides one of the only simple tests available (section [Link]),
but this is dependent on the quality of the original data. Digital topographic surfaces are
some of the most difficult and complex to develop as they represent a real observable
feature. Carter (1988) uses the example that to produce a 500 mbar surface is considered
easier than to produce a topographic surface. It may be equally complex in its terrain,
however since no-one has ever seen one, if a mapped surface appears reasonable we will be
more willing to accept it as correct. However when dealing with topography we have a real
analogue for its structure and therefore something to compare our digital surface with.

One of the many strengths of digital imaging is the ability to manipulate the display of
topography in such a way as to emphasise a particular element, for instance streams or
valleys. This is not easily possible with conventional maps. DEMs also allow hypotheses
regarding landscape to be tested in a reproducible fashion. Another advantage is at regional
scales, where DEMs often contain more high resolution material than can be portrayed using
contour lines. Furthermore, the available digital data may be more easily manipulated for
visual, statistical or other analyses (Mayer, 2000).

In the past, planimetry and measurement of features on paper maps provided the basis for
geomorphic measurements and hydrological applications. Whereas now DEMs are

Chapter 2 – GIS – Techniques and Methods 30

providing the primary data source for use in the analysis of catchment form (Moore et al.,
1991). Together with these advantages of DEMs and GIS, the increasingly widespread
availability of digital datasets and rapidly evolving computational environments have
spurred on the use of such systems for both geomorphic and neotectonic investigations
(Mayer, 2000).

Two different Geographical Information Systems (GIS) are used in this study, ARC/INFO
and ArcView. They are both commercial GIS developed by the Environmental Systems
Research Institute (ESRI) in Redlands, California. The systems work on the same basic
principle of employing data storage in terms of grids and coverages. ARC/INFO provides
the main analysis functions, driven by command line in the UNIX operating system. Within
ARC/INFO it is possible to write code in Arc Macro Language (AML) in order to create
scripts to simplify data extraction and repeat procedures over multiple areas. Some small
scripts were written for data extraction and visualisation, together with scripts developed by
Jamieson (2001) for morphometric parameter extraction. ArcView provides the user-
friendly, menu driven system based within the Windows operating system, which provides
visualisation of the various grids and coverages. The combination of both of these systems
allows for relatively straightforward data manipulation and extraction. External programs
including Excel, Minitab and Matlab are used for statistics and graphing of data.

More recently ESRI have created a more user-friendly package that combines ARC/INFO
and ArcView into one, known as ArcGIS. Due to system limitations and familiarity, only
ARC/INFO and ArcView were used in this thesis.

2.2 Scales and resolutions of DEMs

Digital Elevation Data is presently available for large areas of the globe at a variety of
different resolutions. As discussed previously (section [Link]), these datasets have been
collected in a variety of ways from manual digitising of topographic maps to
photogrammetry using data from remote sensing sources. The intended application will

Chapter 2 – GIS – Techniques and Methods 31

determine which resolution will be most appropriate for use. It is important here to define
the difference between scale and resolution with regard to DEMs as at first glance they
appear to mean the same. Scale is used when describing in a rather qualitative way e.g. fine
or coarse regarding both the scale of landform type features and also DEM scale. Resolution
is used in a much more quantitative way and relates to the exact grid spacing being used, for
example 30m, 50m or 90m. Data is available at a variety of resolutions from very coarse
(1km) to very fine scales (10-30m). This section discusses and describes the available data
and choice of dataset with regard to the characterisation of the macrogeomorphology of the
Pyrenees.

The most readily available data is that at a resolution of 30 Arc Seconds, roughly a 1km grid
spacing. This dataset, known as GTOPO30, is available as a free download from the Internet
in a series of tiles that cover the entire globe. It was created using a variety of cartographic
and remote sensing sources depending on global location. Due to its easy availability and
free download (www1), it is one of the most commonly used DEMs for coarse analysis and
visualisation of large areas. At continental scales it provides an ideal working resolution, but
once individual catchments need to be resolved, the coarse scale becomes a hindrance and
key features cannot be picked out.

The next level up is 3 Arc Second data, with a grid spacing of ~90m, known as Digital
Terrain Elevation Data (DTED) Level 1 data. This resolution of data was previously only
available through ones own creation or purchase. More recently, with the advent of the
United States National Imagery and Mapping Agency’s (NIMA) Shuttle Radar Topography
Mission (SRTM), datasets are now available at this resolution for free download (www2).
This new dataset, was acquired by the Space Shuttle Endeavour in February 2000 and was
collected via synthetic aperture radar and radar interferometry for 80% of the Earth’s land
surface between 60oN and 56oS (NIMA, 2002).

A third higher resolution dataset is also available – DTED Level 2 data with a grid spacing
of 1 Arc Second (~30m). It is not freely available and must be purchased, usually from the
National Mapping Agency of the country of interest. Although this high resolution allows
smaller features in catchments to be represented it often comes at a cost, both in monetary

Chapter 2 – GIS – Techniques and Methods 32

value and in computing time. Compared to Level 1 data, for every cell analysed at this
resolution, nine cells must now be analysed for the equivalent area and so both analysis time
and resultant file sizes will be correspondingly larger.

DEMs at resolutions down to 2m have been created in very specific areas (Zhang and
Montgomery, 1994). However the time involved in data collection and subsequent
processing and the actual purpose of obtaining such high resolutions must be taken into
consideration prior to data acquisition.

For the Pyrenees a variety of data is available: 30 Arc Second GTOPO30 data, 3 Arc Second
DTED Level 1 data (available from the Department of Geography, University of Edinburgh)
and for some small areas, including parts of Catalunya, some 1 Arc Second DTED Level 2
data. In order to decide on which resolution of dataset to use, the purposes and objectives
must first be considered.

A wide number of scales of DEMs exist for use in geomorphic analysis. Different purposes
and scales of feature will require different resolutions (e.g. there is no point in using a 1m
resolution DEM when a 1km resolution DEM will suffice). By adopting this principle
processing time may be saved in the extraction of the parameters of interest and their
subsequent analysis. This is true for all applications of DEMs and is not restricted to only
geomorphic applications. Considerable work in the literature has been undertaken with
respect to the scales and resolutions that are most suitable for particular applications.
Extensive research has emphasised the various scales of features of interest and how they
relate to the resolution of the DEM (Wood, 1996a,b). Studies such as that by Zhang and
Montgomery (1994), on DEM resolutions ranging from 2m to 90m provide constraints on
what resolutions to use for geomorphic and hydrological processes. They conclude that for
catchments of the size order of 1-2 km2 a 10m grid spacing provides a good compromise
between increasing resolution and data volume. Conversely, they suggest, that coarser
resolutions up to 90m may be most appropriate for considering entire orogens over large
areas and long time scales. Wolock and Price (1994) showed that many hydrological
parameters are very sensitive to DEM map scale and resolution. Chan and Rothman (2000)
and Wang and Yin (1997) compared DEMs of different resolutions and concluded that in

Chapter 2 – GIS – Techniques and Methods 33

many cases coarser resolutions were better for estimating several elevation-based parameters
such as slope. Schneider (2001) suggests that the style of initial data interpolation provides
an important first step in deriving morphometric parameters due to discrepancies in methods
of data interpolation between data points. This again comes down to scale-based
characterisation (Wood, 1996b) and the importance of being aware of the scales of
landforms being observed in nature and those that you are trying to measure within the
DEM.

A grid spacing of ~90m or 3 Arc Seconds was used for this study for analysis of the
macrogeomorphology of the Pyrenees.

The final projected dataset is shown in figure 2.1. The image having been hillshaded from
an azimuth of 315o at 45o to the horizontal, in order to highlight the topography. The
geographical limits of the dataset are from the Aquitaine basin in the north to the Montsec
Thrust in the south, which gives sufficient coverage for the analysis.

The DEM can be explained formally as having been created at ~90m scale from aerial,
cartographic and remote sensing sources approximating the resolution quality of DTED 1
(NIMA, 2002). It is therefore a combination of both cartometric and photogrammetric
DEMs (Walker and Willgoose, 1999).

!
Along with digital elevation data for the Pyrenees, certain other forms of digital data are also
available. These are in the form of coverages, i.e. data in a vector or point format.

Chapter 2 – GIS – Techniques and Methods 34

[Link] Digital Chart of the World

The Digital Chart of the World (DCW) provides a freely downloadable database of
coverages of political boundaries, major urban centres and river drainage networks amongst
others (www3). This kind of data, in the case of political boundaries for example, allows the
creation of images that permit the viewer to get a spatial reference frame. Digital
representations of drainage networks allow for comparison with those extracted from the
DEM although this still has the problem of ‘blue lines’ (see later discussion in section 2.3.4
on channel initiation).

[Link] Geological maps

As GIS systems are becoming more widely used in many contrasting disciplines the
conversion of geological maps into a digital format is becoming more common. Examples
of this approach are seen in work by Kirby et al. (2001) in Tibet, where mapping and
digitising have resulted in a geological map within a GIS and within the BGS in the British
Isles. By having the geology of a region within a GIS a series of queries can be undertaken,
which along with improved visualisation of the data with respect to the topography, leads to
a very powerful tool for use in analysis.

Within the Pyrenees, the majority of published geological maps, were surveyed by Dutch
Geologists from the University of Utrecht. For about twenty years from 1948 to 1968 a
comprehensive mapping project of the Pyrenees was undertaken, the result of which was a
series of memoirs culminating in a final report by Zwart in 1979. A series of accurate
geological maps were published on scales generally of 1:50,000 with some 1:25,000 sheets
in more complex areas. Two large summary geological maps of the entire chain were
published at a scale of 1:200,000 (Zwart, 1979) which represent the only small-scale maps
available for the entire Pyrenees as a whole. Individual Geological Surveys (e.g. BRGM and
ICC) have created their own maps, but they are strictly territorial. The Dutch summary
sheets show all the geology and other features of interest, however there are significant
inaccuracies, e.g. longitude values on the sheets show a considerable discrepancy with
respect to Greenwich with the meridian centred on no apparent logical location. Some of the
larger scale sheets have their meridians based around Madrid, though this is not the case on
the 1:200,000 sheets. These factors and the lack of any projection information for the maps

Chapter 2 – GIS – Techniques and Methods 35

made them very difficult to digitise. Attempts to tie the maps in with other maps of known
projection led to consistently wrong placing of the lithologies with respect to the DEM. As a
result a digitised geological map of the Pyrenees was not created.

However recent developments, specifically the IGME 5000 project (International Geological
Map of Europe and Adjacent Areas), is providing a geological GIS for the European region
(ref. www4). It will provide detailed digitised documentation of all European geology (Asch,
2004), enabling studies such as those described above to be undertaken.

2.3 Data pre-processing and quality assessment

The DEMs used in this study, in their most basic format are regular grids of elevation values.
Before any analysis can be done on such a dataset as an approximation for a real terrain
surface, a number of standard procedures need to be undertaken. DEMs need to be
transformed from geographical locations to a projected format and then checked to ensure
that they provide a hydrologically consistent surface and if not, corrected to flow in a
hydrologically consistent manner.

Following these basic steps, the quality of the DEM needs to be assessed, both in terms of its
relationships with the real terrain and any errors that may have occurred during its creation.
These checks are important as key primary morphometric attributes such as elevation, slope
and aspect, are sensitive to errors in the dataset (Bolstad and Stowe, 1994). It is therefore
necessary to know whether features identified in analysis, for instance slope, are real or
simply artefacts of the data.

This section discusses these pre-processing procedures and quality assessment checks
performed on the dataset used for analysis. It also discusses the derivation of river networks
from DEMs and channel initiations in the headwaters of catchments. These features are
important in order for later derivation of geomorphic parameters and basin choice.

Chapter 2 – GIS – Techniques and Methods 36

 " #
Before any pre-processing or analysis can be done on a DEM, it must first be projected onto
a flat surface, from its format as a curved surface as part of the Earth’s spheroid. Therefore
in order to analyse the DEM it is necessary to convert the data from a global geographical
location (in latitude and longitude) on the Earth’s spheroid to a representative location on a
flat surface (ESRI, 1997). There are many different mathematical algorithms designed to
undertake this task, with choice related to global location and the planned application of the
data. Generally there are three types of projection available for use: Conformal, Equal-Area
and Equidistant. Each portray the data in different ways, for example the Conformal
Projection preserves small shapes and local angles well, whilst the Equal-Area and
Equidistant Projections preserve the areas and distances respectively of the displayed
features, often leading to distortion in shape, area and scale. Projections are always a
compromise, preserving certain features at the cost of others.

With regard to this study and geomorphic analysis, it is important to ask what are the factors
of interest and how will these be represented in certain projections? Primary topographic
attributes such as elevation are not affected by projection type, so secondary topographic
attributes become important considerations. Preservation of catchment shapes and their
areas is important and should be considered when choosing a projection type. A Conformal
type projection is judged to be the best compromise between the various factors. As the
Pyrenees are an east-west chain in a mid-latitude location, the projection used was the
Lambert Conformal Conic. The two standard parallels used were 42o20’N and 43o40’N with
a Central Meridian at 1oE. Once projected, the Pyrenean Dataset had a grid spacing of
80.251m relating to both the projection type and the mid-latitude location of the Pyrenees.

" $%
One of the most important pre-processing steps to undertake prior to analysis for either
hydrological or geomorphological purposes is to ensure that the dataset has a hydrologically
consistent surface. Prior to this stage, DEMs often have areas that are internally drained
known as spurious sinks or pits. These local depressions can be large areas, or may be where
a cell in the centre of a 3x3 kernel is at a lower elevation than its neighbouring eight cells.

Chapter 2 – GIS – Techniques and Methods 37

These features may be natural e.g. lakes, but are generally considered to be spurious and
artefacts of grid creation (Jenson and Domingue, 1988) and should be removed prior to
analysis. There is an opposing opinion however, where these features are treated as real,
leading to potential storages within the system (Martz and Garbrecht, 1997; Martz and
deJong 1988). This approach provides scope for modelling the hydrological system’s
storage of water and sediment. However, the more generally accepted technique particularly
when looking at macrogeomorphology is to eliminate these heterogeneities in the river
network. This allows for less than perfect DEMs to be analysed and does not require the
input of entirely error free grids (very rare). If there are any artefacts then they may be
formed in a number of ways, either by incorrect or insufficient data or by an interpolation
technique that does not enforce surface drainage (Hutchinson and Gallant, 2000).

Jenson and Domingue (1988) suggest a procedure for pit removal, which although often
developing flat areas, provides the simplest way to create a depressionless DEM which
allows for hydrologically continuous flow. Their method fills all pits to their pour points,
which assumes that all pits are errors caused by underestimation as opposed to
overestimation (Martz and Garbrecht, 1997). This issue is very difficult to address, so for
simplicity all pits are treated as underestimations, which without extensive field studies is
difficult to determine. This is the most likely cause for flow difficulties, as with an
overestimation the flow is simply likely to go around the cell in high-resolution datasets. A
necessary input to this technique is a grid representing flow directions, the formation of
which is described in section 2.3.3. The end result is that all pits are filled to their pour
points to create a depressionless surface. This surface can then allow uninterrupted flow
through a catchment that is required for hydrological and geomorphological analysis.
Algorithms inherent to ARC/INFO, such as SINK and FILL, enable the pits to be identified
and filled quickly and easily.

& $%
The derivation of flow within DEMs has generated a wealth of literature with suggested
techniques for dealing with this issue. The principal problems arise due to the regular
gridding of DEMs and the bias introduced by grid orientation (Tarboton, 1997). As
previously discussed (section 2.1.1) regular grid spacing DEMs are not the only medium
used for landscape representation, with TINs and Vector based Contour DEMs also being

Chapter 2 – GIS – Techniques and Methods 38

employed (Burrough, 1986; Walker and Willgoose, 1999). Techniques for flow derivation
in these other data storage mediums have also been developed. However as regular grid
spacing DEMs are used in this study, the appropriate flow derivation techniques are focussed
on. As grids are composed of a regular tesselation of rectangular squares, flow in and out of
these cells has to be considered in a structured way.

Early work by O’Callaghan and Mark (1984) and then by Jenson and Domingue (1988)
employed the use of the D8 algorithm for regular grid DEMs. The D8 technique is defined
as the ‘deterministic eight neighbour’ method and allows flow to leave a cell through one of
eight possible directions, i.e. N, NE, E, SE, S, SW W, or NW. Flow is only computed within
a 3x3 cell kernel, with the central cell only knowing about its surrounding eight cells (or
neighbours). This technique has formed the basis of hydrological networks derived from
DEMs due to its robust nature and rapid, straightforward implementation. However it has
been strongly criticised, mostly by hydrologists, where such a simplistic approach in
complex small basins often produces poor representations of the actual observed networks.
The strongest criticisms of the approach are that flow can be up to 22.5o in the wrong
direction and that there is a strong element of grid bias in any derived network, related to the
orientation of the grid axes. In addition within flat areas on DEMs (whether real or an
artefact of DEM creation), D8 will create parallel flow lines which are considered to be
unrealistic of natural systems (O’Callaghan and Mark, 1984; Garbrecht and Martz, 1997).
Solutions to this problem have taken on many different levels of complexity. Work by
Fairfield and Leymarie (1991) suggested using stochastic flow direction methods whilst
Tarboton (1997) proposed representing flow direction as an individual angle which is not
limited by being grouped into one of eight directions (named D∞ method). Alternative
schools of thought suggest that regular gridded DEMs are not as useful for delineating flow
as TINs (Palacios-Vélez and Cuevas-Renaud, 1986) or Contour based DEMs (Moore and
Grayson, 1991). Others have concentrated on flow on flat areas within DEMs and developed
algorithms that more closely approximate flow in natural systems than D8 on flats in DEMs
(Garbrecht and Martz 1997, Martz and Garbrecht 1997, Tribe, 1992). By combining digital
river and lake data with more accurate flow direction algorithms in flat areas, Turcotte et al.
(2001) showed that this again may be another approach to the problem. However many of
these more accurate algorithms and techniques are frequently based on small catchments
where exact flow paths are required for distributed hydrological modelling and where the
accuracy of storm hydrographs and actual flow paths are important e.g. in planning
applications and flood forecasting.

Chapter 2 – GIS – Techniques and Methods 39

[Link] Application to the Pyrenean dataset

As has been shown, many alternative ways to define flow direction have been developed in
an attempt to accurately characterise basin network flow properties. These sophisticated
algorithms tend to be use small scale, high-resolution datasets. As a result the issue of scale
dependence of DEMs is raised (Wood, 1996a,b) and so, if using lower resolution data, many
of the factors/issues with the D8 method become less important (see section 2.2.2 for further
discussion on scale-dependence of DEMs). The grid spacing used in this study is ~90m or 3
Arc Seconds and consequently slight channel fluctuations such as meanders on a scale of
less than 100m, during downslope flow will have little effect. Often when grid DEMs are
used for hydrological purposes, grid resolutions may be as fine as 6.25m (Walker and
Willgoose, 1999). If dealing with resolutions as coarse as 90m, and working at a basin wide
scale, then minor zigzag descent paths of rivers will have little impact on the kinds of data
being extracted. Along with general downslope flow paths, the issue of flat areas in DEMs
must also be considered. Even though the Pyrenees are a mountain belt (and consequently
one would expect very few entirely flat areas) and that the bulk of the analysis for this study
takes place within the Axial Zone of the Pyrenees (area of high relief), flat areas do exist.
Scattered throughout the mountain chain are a series of water reservoirs, predominantly built
for hydroelectric power generation. They appear on the DEM as large flat areas due to post-
damming, DEM data collection and assembly. These flat areas lead to the development of
erroneous flow paths and parallel stream development if the D8 method is used (Fairfield
and Leymarie, 1991; Garbrecht and Martz, 1997). Within the context of this study, they are
considered to be of minor importance, with the main problem being the large flat sections in
subsequently derived river longitudinal profiles (a deviation from the natural system related
to anthropogenic influences). Potentially, minor distance measurement errors may also be
introduced due to straight flow paths being measured as opposed to possibly natural
meandering flow paths when calculating overall stream lengths. At the scale being dealt
with, this is regarded to be of minor importance. As a result it was concluded that the
Deterministic Eight Neighbour (D8) algorithm would suffice for use for this study. It is the
simplest flow derivation algorithm available for use with grid DEMs (Tarboton, 1997) and is
able to rapidly extract stream networks in a straightforward and consistent manner.
However, it is acknowledged that it is inaccurate at a small scale which should be borne in
mind if analysing such features. D8 is also the most straightforward algorithm to apply in
ARC/INFO with it being the basis behind the FLOWDIRECTION command.

Chapter 2 – GIS – Techniques and Methods 40

 ! $%
Another topic of debate within the hydrological community and important to stream network
derivation is ‘Where do channels begin?’ This is not just a problem related to DEM studies
but is an important issue on contour maps. Tribe (1990) noted that when two experienced
geomorphologists were tasked to identify the location of channel heads from contour maps,
their results were consistently different. Standard practice for channel network extraction is
to specify a threshold value from which to initiate a channel – e.g. 50 cells on a DEM
initiates a channel. However variation of threshold areas can lead to differences in
measurements of morphometric parameters (Helmlinger et al., 1993) with an assumed
threshold area often leading to potential errors. Field studies of channel initiation have been
undertaken in an attempt to gain a better understanding of the ways in which channels
initiate (Montgomery and Dietrich, 1988; 1989; 1992). Results of these studies and others
(Dunne, 1980) have shown that channel initiation is dependent on the following factors:
incision by Hortonian overland flow and incision by saturated overland flow, seepage
erosion and shallow landsliding (Moore et al., 1991). It is therefore conceivable that
threshold areas for channels can vary across a drainage basin. If this is then scaled up to a
mountain belt, errors will be introduced, with large amounts of potential variability at a
variety of scales. A common way of estimating threshold areas is to use the ‘blue lines’ on
topographic maps as indicators of amounts of area required to initiate a channel. However, it
has been shown that these rarely represent first order channels (Morisawa, 1957) and, as
noted by Tarboton et al. (1991) usually denote some subjective judgement on the part of the
cartographer. For example in the UK, the Ordnance Survey instruct their surveyors to
represent the network at a constant hydrological state, that is at the ‘normal winter level’
(Ovenden and Gregory, 1980; Wood, 1996a). Montgomery and Foufoula-Georgiou (1993)
suggested that overall threshold values tend to represent the hillslope to valley transition
rather than identify channel heads. However, despite these problems, the use of threshold
areas for channel initiation remain the only rapid large scale technique available for the
delineation of stream networks. When considering the overall river network of an entire
mountain belt system the possible lack of first order streams, especially when using ~90m
grid cells, ceases to be of importance with respect to delineating the general network. Digital
comparisons with ‘blue lines’ on DEMs can also be made using Digital Chart of the World
(DCW) (www3) data for river networks.

Chapter 2 – GIS – Techniques and Methods 41

[Link] Application to the Pyrenean dataset

For the Pyrenean DEM a threshold value of 100 cells was used – equivalent to an area of
0.644 km2. This area gave a good comparison with ‘blue lines’ on both topographic maps
and also with the networks in the DCW dataset. Various other values were experimented
with, but 100 grid cells gave roughly the closest fit to the ‘blue line’ reality. By using
threshold values such as this that compare with ‘blue line’ networks, true first order streams
are being ignored and this should be taken into account in later analysis involving stream
order.

' ( $%
The quality of a DEM is dependent on the quality of the source data and the interpolation
technique applied to it (Hutchinson and Gallant, 2000; Schneider, 2001). The purpose for
which the DEM was created will often determine the type of interpolation used and so care is
required if a DEM is to be used for a purpose other than that for which it was intended
(Schneider, 2001). In order to assess the quality of a DEM a number of simple visual
graphing techniques are employed as demonstrated by Hutchinson and Gallant (2000), on
which the following series of tests is based. Hutchinson and Gallant (2000) note that as
applications of DEMs predominantly depend on their ability to depict landscape shape and
drainage structure, absolute elevation errors cannot provide a full assessment of DEM
quality. Therefore the following tests are used to provide a quality assessment of the dataset.
They are performed once the DEM has been projected and filled to produce a depressionless
surface. Although similar tests could have been performed prior to the pit filling procedure,
it is the depressionless surface which is required for geomorphic analysis.

[Link] Relief/Slope plots

A rapid primary test for a rough guide to DEM quality is to create a hillshaded DEM.
Hillshading or relief shading involves an angled light source at a particular azimuth, which
emphasises the ridges and valleys within a DEM. It was developed initially on paper map
sheets by the Austrian and Swiss schools of cartography making use of the Renaissance
technique of chiaroscuro, where light and dark shades were used to emphasise three

Chapter 2 – GIS – Techniques and Methods 42

dimensional objects (Burrough, 1986). When used as a quality assessment criterion, errors
will show up as dark or light spots, contrasting from the surrounding data, revealing both
random and systematic errors (Hutchinson and Gallant, 2000). Figure 2.2a shows the
Flamisell catchment from the southern Pyrenees (location map in chapter 3, figure 3.23). A
hillshaded image is shown with light from an azimuth of 315o at 45o to the horizontal. It
illustrates that no apparent dark or light spots can be seen on the dataset, so implying no
visible errors using this technique. A hillshaded image of the entire dataset can be seen in
figure 2.1.

A similar rapid test for DEM quality is to perform a slope plot of the DEM and to observe if
any obvious discontinuities or errors emerge. Typical errors seen with poor data, are
obvious discontinuities or step-like patterns. This can be a result of poor initial data or
incorrect choice of projection and interpolation algorithms. Examples for this come from the
Flamisell catchment. Figure 2.2b shows a slope plot of the actual DEM used, with little
obvious error and a continuous surface compared with figure 2.2c showing a slope plot of a
DEM of the Flamisell catchment with poor choice of projection and interpolation algorithms.
The obvious stepped nature and breaks or discontinuities imply poor interpolation and an
unreliable dataset for geomorphic analysis.

[Link] Frequency histograms

Another way of assessing data quality involves extracting frequency histograms of primary
terrain attributes such as elevation and aspect. Frequency histograms are created by
extracting the Value Attribute Tables (VATs) of the grids in question using the TABLES
module of ARC/INFO. Patterns in the data provide clues to poor datasets with, in elevation
data, biases occurring related to contour line intervals (depending on method of DEM
creation) and in aspect data multiples of 45o produced by too simplistic interpolation
algorithms (Wood, 1996a; Hutchinson and Gallant, 2000). Figures 2.3a and 2.3b show plots
of frequency-elevation and frequency-aspect respectively for the entire Pyrenean DEM. The
frequency-elevation histogram shows that overall the data is not especially biased towards
certain values. It is noticeable however that there is some ‘peakedness’ or ‘spikiness’ that
may be indicative of minor errors within the original data of the DEM. Wood (1996a) noted
on studies of frequency-elevation histograms for DEMs of the UK that strong contour bias
was often introduced with peaks at intervals of the originally digitised contours. If this data

Chapter 2 – GIS – Techniques and Methods 43

had been digitised from a slightly different contour interval, then the result would be subtly
different. This test emphasises that DEMs are only as good as the original data from which
they are created. The aspect histogram presents a similar story, with pronounced peaks at
45o intervals. This implies that there is a slight deficiency in data interpolation, but overall
the data appears reliable.

These tests have been performed on a very coarse scale and in order to look at the data in
more detail, the equivalent plots have been created for the Flamisell catchment (figure 2.4a,
2.4b). Within the elevation data a finer pattern of spikiness is observed, implying some
degree of error within the original dataset. However there seems to be little of the regularity
observed in plots derived solely from digitised contours. In the aspect histogram, the 45o
errors seen in the plot for the entire Pyrenees are not as prominent or obvious.

So overall, it would appear that although there are minor errors in the dataset, in relation to
geomorphic analysis they should have a minimal impact.

[Link] Comparison to topographic maps

A test to compare DEMs with the real terrain can be made by creating a contour plot of the
grid data and then contrasting it with contour data from a topographic map (Hutchinson and
Gallant, 2000). An example is shown of a section of the Flamisell catchment showing both
the 1:50,000 topographic map and derived contours from the DEM with contour spacing at
20m intervals (figure 2.5). Although visual comparison detects the overall similarity of the
two plots, the derived contours tend to generalise the real map contours to a degree as would
be expected from a relatively coarse DEM. Generally it is suggested that the DEM provides
a good approximation for the terrain it represents.

[Link] Summary

In summary, the various tests applied to the digital elevation data in the Pyrenees imply that
the dataset is overall relatively consistent and generally suitable for geomorphic analysis. It
has identified a few problems based on the frequency histograms, but overall appears to be

Chapter 2 – GIS – Techniques and Methods 44

suitable for the intended task. In conclusion, it is important to note that none of the data
employed here is without some degree of error, but if these errors are small and are known,
then their effect can be minimised.

) "
When attempting to classify the morphometric parameters of mountain chains, it is assumed
that we can treat the system as if it is in its natural state. However in the Pyrenees, the
human impact on the environment must be considered. Large programs of dam building
have altered the forms of valley floors and river longitudinal profiles. This has had the effect
of creating large areas of slope at 0o, which is unlikely within natural systems. By creation
of these flat areas, parallel streams develop due to choice of flow direction algorithms
(section 2.3.3).

The DEM used was created from data obtained subsequent to these developments and so
care must be taken when interpreting the results to ensure that they are not artefacts of
human intervention on the landscape.

2.4 Geomorphic parameters and their extraction

The geomorphic analysis compares the northern and southern flanks of the Pyrenees. In
order to break down this type of analysis the flanks were divided into a series of basins to
form the basis of comparison across the main drainage divide as well as comparison along
strike. A technique of basin derivation is first defined, providing the basis for the subsequent
extraction of geomorphic parameters.

!
A series of catchments within the Pyrenees were defined on both the northern (retro-) and
southern (pro-) flanks of the chain for comparison of their respective morphometric
attributes. As the study was designed to determine the relative influences of the main

Chapter 2 – GIS – Techniques and Methods 45

controls on orogenesis (tectonics, climate and lithology), the basin pour points were not
defined on the basis of lithological contrasts or major through-going structures such as the
North Pyrenean Fault. Instead, an independent way of choosing basins was used by defining
pour points on the basis of stream order. Various techniques for stream ordering have been
developed in the literature over the past sixty years by fluvial geomorphologists (Horton,
1945; Strahler, 1952; Shreve, 1967). One of the simplest and most commonly used systems
is that developed by Strahler (1952), which allows a stream to increase in order only when
two streams of the same order intersect. For the purposes of basin delineation, pour points
were selected at the confluence of a fourth order with a fifth order stream segment. This
allowed basins to be defined in one of the two ways illustrated in figure 2.6. These are
defined as either at the confluence of two fourth order streams creating two basins (figure
2.6a) or at the confluence of a fourth order with a fifth order stream creating one basin
(figure 2.6b). An extra criterion was also stipulated, in that the basins had to back onto the
main Pyrenean drainage divide.

This classification led to the delineation of thirty-seven independent catchments, nineteen in

the north and eighteen in the south (figure 2.7). Catchments were numbered from one (NW)
through to nineteen (NE) north of the main drainage divide and from twenty (SE) through to
thirty-seven (SW) south of the main drainage divide. Thirty-seven catchments provide a
statistically acceptable number for comparisons across the main drainage divide.

! *%
Geomorphic or morphometric parameters have been extracted from topographic data in a
quantitative way beginning with the work of Horton (1945). Strahler (1958) summarised the
extraction of what we now describe as geomorphic parameters by reviewing the various
parameters developed by a variety of authors (Strahler, 1958 and reference therein). This
development allowed for a wide-ranging suite of landscape properties to be measured in a
quantitative or dimensional way as opposed to a previously purely descriptive approach. In
the past the measurement of these properties was often a non-trivial task, with very time-
consuming measurements being required from topographic maps. This led to these kinds of
detailed analyses only being undertaken on the scale of small individual catchment (e.g. the
Perth Amboy Badlands, Schumm, 1956). The advent of DEMs has created a medium for
rapid and reproducible extraction of these properties. As a consequence, larger areas may be

Chapter 2 – GIS – Techniques and Methods 46

rapidly analysed and allow comparisons between areas that would previously have required a
much higher workload.

[Link] Parameter types

Geomorphic parameters can be roughly divided into three groups: areal, relief and network
properties (Strahler, 1958). Areal properties are those that express the overall plan form and
dimensions of drainage basins. Relief properties express features involving elevation
differences whilst network properties are those that concern the channel network both in
terms of topological and geometrical aspects (Strahler, 1958; Summerfield, 1991). The three
groups of these parameters are given in table 2.1 along with their definitions and units of
measurement.

Additional geomorphic features not included in the table are river longitudinal profiles and
geophysical relief properties that fall across more than one of these descriptive groups. Their
definition and extraction is considered and described in sections [Link] and [Link].

[Link] Technique for extraction

The techniques on which morphometric parameters are extracted in this study are based on
those created by Jamieson (2001) using an AML that extracts certain features from the DEM.
The parameters were extracted in a basin by basin fashion, with the AML, known as Drainall
(Jamieson, 2001). It combines together many basic analytical techniques in ARC, such as
deriving slope and mean elevations, and runs them all together in a complete and ordered
fashion extracting all parameters of interest from the basin. The AML uses a combination of
grids and coverages to calculate and extract the parameters of interest. Although the features
could have been extracted independently from the DEM, the AML provides a time saving
device as well as introducing less chance of error between catchments. The final output is a
text file that can be imported into external programs giving the results of all extracted
parameters in a basin by basin format.

Chapter 2 – GIS – Techniques and Methods 47

Following extraction, parameters were tested for statistical significance using Minitab to
determine whether there was any difference between northern and southern flanks of the
chain. Discussion of the statistical tests used is given in section 2.5.

Slope and elevation data for basins was unloaded directly from the TABLES module of ARC
for the creation of frequency histograms.

[Link] River longitudinal profile extraction

River longitudinal profiles are considered to be important indicators of the interactions

between climate, erosion and tectonics (Seeber and Gornitz, 1983). Consequently, it is
important to be able to extract such profiles from the DEM in order to assess the effects of
the above controls.

Two alternative methods were used to extract river profiles from the DEM, one following the
work of Jamieson (2001) and the other using profiling techniques inherent to ARC/INFO.
Both techniques extract ASCII text files of data from ARC/INFO, which require subsequent
plotting in a separate graphing package (such as Microsoft Excel or Matlab).

The technique described by Jamieson (2001) uses a section of code in AML known as
Profileall that extracts all possible river longitudinal profiles within a drainage basin. It is
based upon a least cost path analysis, which computes the path of least resistance down slope
(i.e. the downstream flow path). For each individual stream it then exports data in an ASCII
text file that records elevation changes and downstream distances allowing for a longitudinal
profile to be plotted external to the GIS. There are three major disadvantages with this
technique. One is related to the large volume of data produced. In many basins in the
Pyrenees there are more then 200 first order streams giving rise to 200 data files. The
streams are numbered in a logical way during the extraction, with number 1 assigned to the
stream with its headwaters the furthest north and number 200 assigned to the one with its
headwaters furthest south. In order to work out which number corresponds to which stream
in the basin requires individual imaging of the traces of the river profiles. A second is
related to the way that the GRID module of ARC/INFO samples data. Grids are analysed in
rows, running from left to right. Once one row has been completed, sampling continues on
the next row down moving from left to right. This introduces errors where rivers do not run

Chapter 2 – GIS – Techniques and Methods 48

perfectly perpendicular to the rows or if one flows from right to left. If a river flows in this
way it is sampled backwards and subsequently the extracted profile drawn will flow
upstream, instituting an element of grid bias into the sampling procedure. Thirdly, when
plotting the profiles, errors may be introduced by assuming that the cells are crossed straight
along the grid. This may lead to errors as flow may move in other directions, i.e. may move
diagonally across cells. For example, if every cell of 1x1m dimensions only counts as 1,
then if the river crosses diagonally the distance should be 2. Significant distance is
therefore lost in the overall river length, particularly when this is scaled up to actual grid
sizes.

All three of these issues can be resolved in order to allow this method to be used. The only
solution to the first problem is time consuming and requires the drainage network for the
basin is printed out and then each individual stream imaged and matched to a stream on the
printed drainage network. The second may be corrected by sorting the data by altitude in
Excel, after computing the distance between each individual sampling point, so that the
downstream distance in not affected. The third can be resolved by extracting the x and y
coordinates of the elevation sampling point. The difference in x and y between each
elevation point can be calculated and then using this data, Pythagoras’s theorem can be used
to find out the exact distance travelled between the last elevation point and the present point,
allowing a more accurate downstream distance to be calculated.

The second technique makes use of the TRACE and SURFACEPROFILE commands in
ARC/INFO and then using UNLOAD, exports the extracted data to an ASCII file. It works
by selecting the required stream and then exporting the data for that stream alone. This
process works very effectively and will graphically display the extracted river before it is to
be exported for checking purposes. The technique also has its downfalls, when very long
rivers with >5000 sampling points are to be extracted. In this instance the command will
split the rivers up into segments and export it in this way, requiring re-stitching in another
program. In addition, if all rivers in a catchment are required to be extracted then it can be a
time consuming process. It does have the advantage however, that there is little need to do
much processing outside of ARC/INFO as it calculates distances based on actual distance
travelled through a cell, not simply cell diameter.

A comparison of both techniques was undertaken which illustrated almost identical results
for both elevations and downstream distances. This study was performed on the Flamisell

Chapter 2 – GIS – Techniques and Methods 49

valley with good correlations observed between the two techniques. Both techniques have
their own advantages and disadvantages, with Profileall providing rapid extraction of all
longitudinal profiles, but with considerable editing required later and SURFACEPROFILE
being relatively slow at extracting multiple profiles, but requiring little subsequent data
editing.

Both methods have the disadvantage inherent to all techniques of exporting river profiles
from DEMs which often produces a stepped profile. Techniques, such as smoothing with
polynomials have been proposed to remove this effect, but this can lead to partial uphill flow
in some areas related to choice of polynomial order. Techniques such as using a moving
average may achieve a more satisfactory result, but are not without error. As a result, no
data smoothing was carried out on the river profiles and it is noted that the very fine stepped
nature is an artefact of its method of production.

In this study a combination of techniques was used depending on the purpose, with the
Profileall technique giving results when entire catchments were to be used and
SURFACEPROFILE used when only one or two profiles were required to be extracted.

[Link] Rubber sheeting and relief plot generation

The suggestion that within mountain belts there are significant amounts of missing mass,
related to valley erosion leading to isostatic peak uplift (Montgomery, 1994), has led to
attempts to quantify this so called ‘missing mass’. A technique developed by Brocklehurst
and Whipple (2002) provides a means by which to quantify this mass and gives a further
indication of whether catchments have been shaped by glacial or fluvial erosion. It is based
on the premise of measuring what is known as the sub-ridgeline or geophysical relief of
basins (Small and Anderson, 1998). The geophysical relief of a catchment is given by the
volume of material ‘missing’ below peaks and ridges divided by the surface area. A
technique similar to that used by Brocklehurst and Whipple (2002) is applied here to
calculate the relief of individual basins. The method takes the form of a series of steps:

1) Isolate the ridgelines of a given drainage basin – i.e. create a grid of elevation values
outlining the basin, using the basin outline coverage to clip the complete elevation grid.

Chapter 2 – GIS – Techniques and Methods 50

2) A smooth surface is then interpolated for the extent of the basin based upon the outlines
of these ridges.

3) This interpolated smooth surface next has the actual elevation values deducted from it
and a plot is developed of ‘sub-ridgeline’ relief.

4) The geophysical relief is then calculated (in m) as the sub-ridgeline relief per unit area –
i.e. volume of missing mass divided by the area of the catchment (m3/m2 = m).

Brocklehurst and Whipple (2002) defined characteristic features observed in the sub-
ridgeline relief plots that allowed a distinction to be drawn between fluvial and glacial
catchments. In fluvial catchments the relief approximated a bullseye type pattern, with a
high relief core and then declining relief away from this point. In contrast, the glacial
catchments, although also having a high relief core, have fingers of high relief propagating
away from this core with high relief found much further from the basin pour points in the
glacial than the fluvial case (see examples in figure 2.8).

As well as creating these plots of relief structure a series of profiles were also derived. The
river long profile was converted into a coverage and then used to create profiles of the
interpolated surface and also of the sub-ridgeline relief (figure 2.9). For these profiles the
SURFACEPROFILE command suite in ARC was used as described in section [Link].

2.5 Statistics

Nearly all data extracted from the DEM was ratio class data, i.e. it is on a scale of continuous
measurement that has a true zero such as slope (Shaw and Wheeler, 1986). Therefore with
this collection of parametric data, the parametric t-test would provide an acceptable
statistical test if the data was normally distributed. The majority of results of the data
extracted from the Pyrenees have either strongly skewed distributions or many peaks (multi-
modal). Non-parametric tests have the advantage of being distribution free and applicable
regardless of the distribution of the data. As a result the Mann-Whitney U test is the most
applicable for use with this dataset.

Chapter 2 – GIS – Techniques and Methods 51

The extracted data are required to be tested for statistical significance to determine if there is
any difference in extracted parameters across the drainage divide in the Pyrenees. The null
hypothesis was proposed to be one of no difference (Davis,1986).

Null Hypothesis (H0) – That there is no difference between northern and southern flanks of
the Pyrenees.

Alternative Hypothesis (H1) – That there is a difference between the northern and southern
flanks of the Pyrenees.

Since all that is of importance is any difference between north and south a ‘two-tailed’ test is
used. This means that the rejection region for H0 lies at both ends of the test statistic
distribution. As the only factor of interest is whether there is a difference, this is the most
applicable to the task.

The level of significance (α) for the test was set at 95% confidence i.e. α = 0.05.

2.6 Summary

This chapter has explained the theory behind the use of DEMs in geomorphic analysis and
the techniques by which this may be performed. It has highlighted both advantages and
disadvantages in the use of such techniques, and the inherent flaws that creep into analysis
both in older cartographic techniques and in the use of DEMs. The following chapter puts
the techniques described here into practise and considers the macrogeomorphic form of the
Pyrenean mountain belt.

Chapter 2 – GIS – Techniques and Methods 52

Chapter 2 – GIS – Techniques and Methods 53
Chapter 2 – GIS – Techniques and Methods 54
Chapter 2 – GIS – Techniques and Methods 55
Chapter 2 – GIS – Techniques and Methods 56
Chapter 2 – GIS – Techniques and Methods 57
Chapter 2 – GIS – Techniques and Methods 58
Chapter 2 – GIS – Techniques and Methods 59
Chapter 2 – GIS – Techniques and Methods 60
Chapter 2 – GIS – Techniques and Methods 61
Chapter 2 – GIS – Techniques and Methods 62