Original Paper

Landslides (2004) 1:73–81 Lulseged Ayalew · Hiromitsu Yamagishi · Norimitsu Ugawa

DOI 10.1007/s10346-003-0006-9
Received: 11 September 2003
Accepted: 7 November 2003 Landslide susceptibility mapping using GIS-based
Published online: 21 February 2004

Springer-Verlag 2004 weighted linear combination, the case in Tsugawa area
of Agano River, Niigata Prefecture, Japan

Abstract A spatial database of 791 landslides is analyzed using subjected to landsliding, is based on the assumption that
GIS to map landslide susceptibility in Tsugawa area of Agano forthcoming landslides occur under similar conditions of those
River. Data from six landslide-controlling parameters namely observed in the past (Guzzetti et al. 1999). The process of GIS-
lithology, slope gradient, aspect, elevation, and plan and profile aided landslide susceptibility mapping at present involves several
curvatures are coded and inserted into the GIS. Later, an index- methods that can be considered as either qualitative or quanti-
based approach is adopted both to put the various classes of the tative. Qualitative methods depend on expert opinions, and are
six parameters in order of their significance to the process of often useful for regional assessments (Soeters and van Westen
landsliding and weigh the impact of one parameter against 1996; Aleotti and Chowdhury 1999). Quantitative methods rely on
another. Applying primary and secondary-level weights, a con- observed relationships between controlling factors and landslides
tinuous scale of numerical indices is obtained with which the (Guzzetti et al. 1999).
study area is divided into five classes of landslide susceptibility. In this study, we used a method known as “Weighted Linear
Slope gradient and elevation are found to be important to Combination (WLC)”, which can be taken as a hybrid between
delineate flatlands that will in no way be subjected to slope failure. qualitative and quantitative methods. Like quantitative approach-
The area which is at high scale of susceptibility lies on mid-slope es, such as bivariable statistical methods, it starts with compar-
mountains where relatively weak rocks such as sandstone, ison of data-layers corresponding to landslide controlling
mudstone and tuff are outcropping as one unit. parameters and the landslide inventory map and involves the
computation of landslide density to assign primary-level weights
Keywords Landslide · Susceptibility · GIS · Agano River · Japan for each class of a particular parameter. Then, it turns to
procedures common in qualitative methods for an application of
Introduction secondary-level weights to the parameters themselves using a
Landslides are common along Agano River, in Niigata Prefecture pair-wise comparison matrix. The final steps of this method are
of Japan. Despite the presence of dense vegetation and little human the combination of all weighted layers into a single map and the
interference, some areas are currently too sensitive for high classification of the scores of this map into landslide suscepti-
precipitation and become sites of active landsliding. Every year, a bility categories which are neither new nor unfamiliar to both
significant amount of land near or away from the river is changed qualitative and quantitative approaches.
into unstable ground. As part of the solution to the problem,
localized studies are repeatedly conducted by engineering com- Description of landslides
panies, and a variety of remedial measures are installed along It is now becoming universal that susceptibility mapping starts
roads, river banks, and ledges of mountains. In addition, there are with the inventory of landslides. In this study, landslides were
efforts for a regional landslide “hazard” mapping and analyses in mapped by interpreting the 1:20000 scale aerial photographs
various sections of Agano River by a group of experts from the taken in 1971 and 1976. As a supplement for this, a red relief image
Landslide Society of Japan (Higaki 2003; Chigira 2003). This study (RRI) of the area was obtained from Asia Air Survey Company
is an extension of not only these efforts but also another GIS-aided (Tokyo), and was subjected to a variety of remote sensing
susceptibility mapping carried out in Kakuda-Yahiko Mountains analytical techniques. Band 2 and 3 of this image were especially
of the same Niigata Prefecture (Ayalew and Yamagishi 2003). useful to confirm the boundary of landslides using image
The Agano River has a channel length of 210 km and a enhancement methodologies such as contrast stretch and digital
catchment area of about 7,710 km2. It flows from east to west and filtering. Band 1 allowed marking out ridges and stream lines.
constitutes a large drainage system in Fukushima and Niigata In total, we were able to make an inventory of 791 landslides
Prefectures of Japan. The part of Agano River selected for this and later described them using the system introduced by Cruden
study is known as Tsugawa area. It is located about 50 km east of and Varnes (1996). The total area covered by landslides is about
Niigata City (Fig. 1), and covers four 1:25000 sheets of the 53 km2, nearly 13% of the area under study. For reasons linked to
National Geographic Institute of Japan with a total area of around geomorphology and geology, many of the landslides are located in
410.18 km2. Precipitation is high in Tsugawa, and comes in the the eastern half of the study area (Fig. 2). The smallest landslide is
form of snow and rainfall. A 20-year record up to the year 2000 0.13 km2 in extent, while the largest one has an area of 0.9 km2.
yields an annual mean precipitation of 2,293 mm and an average A series of field surveys have been conducted with an aim of
snow depth of 111 cm. A direct result of this high precipitation is studying the characteristics of landslides at different times.
thick vegetation available throughout the region. Accordingly, as far as the type of movement is concerned, it was
It is known that a landslide susceptibility map relies on a observed that slides are the dominant forms of slope failure. In
rather complex knowledge of slope movements and their terms of depth, many landslides are shallow, although some of
controlling factors. Mapping of areas, which are not currently those occurred on hillsides are relatively deep-seated. The state or

Landslides 1 · 2004 73
 Original Paper

Fig. 1 Location map of the study area

74 Landslides 1 · 2004
Fig. 2 Landslide distribution in Tsugawa area of Agano River

activity of landslides was in such a way that most of the mapped

landslides are relict and stabilized. This was determined on the
basis of signs of old mass movements such as crescent-shaped
scarps, abnormal bulges on inclined slopes and hummocky
surfaces. However, a number of landslides in the central part of
the study area near Agano River, which contain degraded
channels and blocks of bedrocks whose only source appears to
be upslope, are active. Figure 3 presents the effect of one of these
active landslides on a road that is located in the southeastern part
of the study area, close to Agano River channel.

Event-controlling parameters
The occurrence of landslides in general is largely a function of the
interaction of natural phenomena such as unfavorable lithology,
stratigraphic sequence, structural makeup, geomorphological
setting, earthquake, rainfall, etc. In GIS-based analyses, these
phenomena which are directly or indirectly related with the
formation of landslides (the event) are commonly known as Fig. 3 A landslide occurred at a locality called Iwatsu on October, 2001 in the
event-controlling parameters. Although, it is believed that the southeastern part of the study area, close to Agano River channel
accuracy of susceptibility mapping increases when all event-
controlling parameters are included in the analytical process, it is
usually difficult to get so, because detailed data is hard to find. For The lithological makeup of the study area
this reason, analyses in this study depend only on lithology and According to a 1:50000 geological map compiled by Hasegawa
the topographic attributes of the region such as elevation, slope (1983), more than 20 rock units are present in the study area. The
gradient, aspect and curvature. A discussion on these parameters southwestern part is composed of a Pre-Cretaceous complex of
with regard to their effect on the process of landsliding is given sedimentary rocks including mudstone, sandstone, chert, lime-
below. stone and greentuffes. To the north of this complex, a massive

Landslides 1 · 2004 75
 Original Paper

Fig. 4 The simplified form of the lithological map of the study area modified from Hasegawa (1983). Symbols such as AL, AN, etc, are useful to read the corresponding
landslide densities in Fig. 5

granite of Cretaceous age occupies more than half of the western Ignoring stratigraphic content and focusing on lithological
part of the study area. The eastern limit of both the Pre- similarity, the 20 rock units shown on the geological map
Cretaceous complex and the massive granite is marked by two compiled by Hasegawa (1983) were in this study simplified into 11
major faults that run in the NNE-SSW direction, almost parallel to as shown in Fig. 4. Our study makes a distinction between debris
each other. Further east, the area is covered mainly by Neogene flow deposits and other types of mass movements. Hence,
formations (Uemura and Yamada 1988), composed of sedimen- although there are some places in the study area where old
tary rocks such as conglomerates and sandstones and a mixture debris flow deposits are present, as it is shown in Fig. 4, these
of rhyolite lavas, pyroclastics deposits and perlitic hyaloclastic materials are not included in the landslide distribution map in
breccias intruded by rhyolitic, andesitic, or basaltic dykes. Fig. 2. The reason is that these deposits are mixed with in-situ

76 Landslides 1 · 2004
rock materials, and are highly stabilized and significantly upward and downward from these elevation marks and becomes
lithified, and we found it not reasonable to associate them with 0 above the height of 800 m (Fig. 5).
other types of mass movements. As far as gradient is concerned, five classes of slope angles
The GIS work was started by rasterizing Fig. 4 and the were established and the corresponding landslide densities were
landslide distribution map (Fig. 2). The pixel size of these raster calculated. It is found that landslide density is high (42%) on
maps was determined by the dimension of the digital elevation slopes with gradients in the range of 2.5–15
(Fig. 5), and decreases
model used in this study. Next, the two maps were overlapped and with an increase or decrease in slope angle. The aspect of the
landslide pixels lying on each of the rock units were counted, and topography is also an essential component of stability analyses,
the areas that they cover were calculated. Then, the ratios between and it was observed that landslide density is higher in east-facing
each of these areas and the total areas covered by the slopes than in their west oriented counterparts, most probably
corresponding rock units were computed and changed into due to the action of erosion by westward flowing streams.
percentages to form what is known as landslide densities in this In addition, both plan and profile curvatures were inputs in
study. The landslide densities were used as a means of rating rock this study. The type of landslide occurring in a certain landform
units against susceptibility. is a function of plan curvature (Ayalew and Yamagishi 2004). This
In the study area, about 40% of the slopes affected by is true because debris and mudflows usually occur when the
landslides are covered by Tertiary pumiceous tuffs. But, in terms lateral profile is concave. Profile curvature governs the run-out
of landslide density, the sandstone layers present mainly in the distance of disturbed materials. Some other properties like the
eastern part account for 23% (Fig. 5). With a landslide density of amount of material involved, the frequency of occurrence, the
20%, the Tertiary pumiceous tuffs immediately follows the shape of rupture surfaces, etc, are in addition to other factors,
sandstones and the rock unit corresponding to rhyolite lavas functions of both plan and profile curvatures.
and dykes was third by having a density of 18%. On the lower side In the study area, data on profile and plan curvatures allowed
of the spectrum, the alluvial deposits and the massive granite dividing slopes into classes of concave, planar and convex
have a negligible amount of landslide densities, and are believed topographic surfaces. Concave slope facets are characteristics of
to play a little role in the process of landsliding. landscapes affected by old mass movements. Topographic
surfaces made up of convex plan curvatures are present in places
Topographical constraints where existing cliffs portray strong lithological variations in the
Thanks to the advance in technology, digital elevation models horizontal direction. They are also common in locations where
(DEMs) are now standard tools for landslide analyses. This study terrains are dissected by sub-parallel, deeply incised ravines.
took the full advantage of this, and significant terrain attributes Planar topographic surfaces exist in localities where the under-
were produced from a 10mX10m DEM obtained from GISMAP of lying geology is relatively homogeneous allowing the slope to dip
Hokkaido Chizu Co. Ltd. GIS technology permits patterns of in one direction.
instability to be resolved and mapped at the scale of the DEM. There is a general consensus that high probability of failure
This means, with the DEM employed in this study, it was possible exists when at least one of the slope curvatures is concave because
to conduct a relatively fine scale analyses which can include slope of the possibility for the concentration of groundwater in a deep
failures as small as 100 m2 in extent. soil stratum. Besides, many researchers agree that landslides on
Primary topographic attributes which can be produced from a convex topographic surfaces need long time to develop since the
DEM are generally first and second derivatives of elevation data, slope geometry forces water to drain away from the site. But in
and include parameters such as slope, aspect, profile and plan the study area, it was observed that landslides occur both on
curvatures (Moore et al. 1991). The effect of slope and aspect on convex and concave slopes. In fact, there was no difference in
landslides is widely documented by Lee and Min (2001) and Dai landslide densities between concave and convex profile curva-
and Lee (2002). But, little is said about the influence of profile and tures, although a significant discrepancy was observed in the case
plan curvatures in the literature. The profile curvature represents of concave and convex plan curvatures (Fig. 5).
the rate of change of slope for each cell in the direction of dipping.
The plan curvature shows the bending of the surface perpendic- Weighted linear combination (WLC)
ular to the slope direction. Together with other factors, plan and As it is stated earlier, weighted linear combination (WLC) is a
profile curvatures control the flow of water in and out of slopes concept where event-controlling parameters can be combined by
and are, therefore, important in the study of landslides. applying primary- and secondary-level weights. The primary-
In the study area, highlands with altitudes greater than level weights are rule-based in that ratings are given to each class
1,000 m correspond to the summits of granite in the northeastern of a parameter on the basis of a certain criterion. In this study,
end and some ridges in the central part. In between and away this criterion is the landslide density, a ratio between the area
from such peaks, the topography is rough and consists of twisting occupied by landslide pixels on a class of a certain parameter and
valley walls and cliffs. Undoubtedly, the main features of the study the total area of that class, changed into percentage. The
area are landslides which occupy ample sectors of these places. In secondary-level (factor) weights are, however, opinion-based
order to assess the effect of altitude on landslide distribution, we scores, which determine the degree of tradeoff of one parameter
classified the elevation map of the project site into 22 categories against another. The WLC adopted in this study shares some
on a 50-m basis. Then, we calculated landslide densities for each similarities with the type of AHP (analytic hierarchy process)
class of elevation adopting the system discussed in the lithology used by Yagi (2003). The difference is that the latter considers
section above. It became apparent that landslide density is greater only one-level weighting system developed by collecting expert
at ranges of altitudes that correspond to a little higher than the opinions, the ratings of which might correspond to secondary-
foot-hills of mountains (201–250 m). Density decreases both level (factor) weights of this study.

Landslides 1 · 2004 77
 Original Paper
Fig. 5 Bar graphs showing landslide
densities. Sub-divisions of the X-axis
represent classes of the six landslide-
controlling parameters

78 Landslides 1 · 2004
Table 1 A pair-wise comparison matrix for calculating factor weights
Pair-wise comparison 9 point continuous rating scale
Extremely V. Strongly Strongly Moderately Equally Moderately Strongly V. Strongly Extremely
Less important Important More Important
1/9 1/7 1/5 1/3 1 3 5 7 9
Aspect Elevation Lithology Plan curvature Profile curvature Slope gradient Factor weights Consistency ratio (CR)
Aspect 1 0.0657 0.07
Elevation 3 1 0.1929
Lithology 3 3 1 0.2569
Plan curvature 1 1/5 1/3 1 0.0715
Profile curvature 3 1 1 1 1 0.1478
Slope gradient 3 1 1 5 3 1 0.2651

A variety of techniques are available to develop factor weights. (Table 1). If the inverse were the case, we would consider 1/3 as
When the numbers of parameters are more than two like the case a rating value.
in this study, however, breaking the information down into simple In this study, the pair-wise comparison matrix contained 36
pair-wise comparisons in which only two factors can be boxes. Since pair-wise comparison matrices are symmetrical in
considered at a time can greatly facilitate the weighting process. nature, only 21 values were needed to fill in the diagonal and the
This technique has been described by Saaty (1988, 1994) and Saaty lower triangular half of the matrix. Then, in order to compute the
and Vargas (2001) in the context of decision making processes. principal eigenvector of the matrix and obtain a best-fit set of
The idea behind is that weights that converge to one are derived factor weights automatically in the way Saaty (1994) and Saaty
through a principal eigenvector of a square reciprocal matrix of and Vargas (2001) have described, raster maps produced by
pair-wise comparisons between event-controlling parameters. A combining the parameters with landslide distribution were
pair-wise comparison matrix has also been used by Juang et al. necessary.
(1992) to map slope failure potential using fuzzy sets. The final result consists of the factor weights and a calculated
In our case, secondary-level or factor weights that can capture consistency ratio (CR), as it is shown in Table 1. The CR is a ratio
the relative importance of one parameter relative to another were between the matrixs consistency index and random index, and in
established on the basis of a 9-point recording scale, which general ranges from 0 to 1. The random index is the average
represents nine linguistic expressions and their corresponding consistence index obtained by generating large numbers of
numerical values. The linguistic expressions explain the fact that random matrices. A CR close to 0 indicates the high probability
the state of knowledge is imperfect, while the numerical values are that the weights were generated randomly (Saaty 1988, 1994).
quantified translations useful for calculating factor weights. Values less than 0.1 are often acceptable, although this depends on
Science still lacks a direct way of evaluating intuition or the objective of the study. A CR of 0.07 in this study is good
expressions, and the validity of the numerical values may best enough to recognize the factor weights as reasonable values.
be judged by the factor weights and the consistency of the
calculation process. The landslide susceptibility map
The complete lists of expressions and numerical grades In seeking a landslide susceptibility map, the primary-level
adopted in this study are given in Table 1. To find the appropriate weights corresponding to classes of parameters were multiplied
linguistic expressions, one may use a pixel by pixel investigation, by secondary-level or factor weights to produce a continuous
divide the area into small partitions of significant extent or scale of numerical values. Dividing these values into susceptibility
consider the project site as a single entity and perform analyses. classes was, however, not easy as there are no statistical rules
In all cases, the assessments are to a large extent subjective, and which can guide categorize continuous data automatically. There
since the knowledge source varies from person to person, it is are some mathematical methods proposed by Scott (1979) and
always true that the best judgment comes from an individual with Friedman and Diaconis (1981), which rely on the optimum bin
good expertise. width classification of the histogram, but they are ineffective to
In this study, we divided the area into around 41 partitions of multimodal distributions (Szen and Doyuran 2004). The prob-
100,000 pixels. Then, we investigated how landslide pixels are lem of changing continuous data into two or more categories
distributed in each partition, and attempted to determine the remains always unclear in landslide susceptibility mapping,
effect of a certain parameter on slope failures compared to because most authors use their expert opinion to develop class
another. Next, we assigned appropriate linguistic expressions for boundaries.
each assessment and put the corresponding numerical grades. In this study, we took into consideration four systems of
Afterwards, we summed up these grades to find a total value, and classifiers that use the so-called natural breaks, quantiles, equal
divided this value by the number of partitions to determine the intervals, and standard deviations, and attempted to choose the
average. For example, after assessing each partition for the role of one that best suits the objectives of our study. While these
elevation and aspect, we came to the conclusion that the former is classifiers are well established in statistics, they often lead to
on average “moderately more important” than the later in different results, as they make very different statements about
forming landslides. Hence, we put three in the pair-wise how values should be divided. The classification scheme that
comparison matrix box that correlates the two parameters relies on natural breaks for example identifies break points by

Landslides 1 · 2004 79
 Original Paper
Fig. 6 The histogram of the numerical
indices (result) obtained

Fig. 7 The landslide susceptibility map

looking for patterns inherent in the data. In quantile classifica- relatively big jumps in data values. The histogram of the
tion, features are grouped by equal numbers in each class. The numerical values obtained in this study (Fig. 6) is multimodal
equal interval scheme divides the range of values into equal-sized and has a number of peaks but shows no empty class intervals or
subdivisions. When a standard deviation is used, data is classified jumps. So, natural breaks were also not appropriate. The last
based on the amount a value varies from the mean. alternative was the classification scheme that uses standard
The quantile classification has a disadvantage in that it places deviations. This method has a certain merit in that it uses the
widely different values into the same class. Hence, it was not used mean value to generate class breaks and allowed us to divide the
in this study. Using equal intervals was also found to be not result of this study into five categories by adding or subtracting 1
practical since it emphasizes one class of susceptibility relative to standard deviation at a time (Fig. 6). As is shown in Fig. 7, the five
others. In natural breaks, boundaries are set where there are categories correspond to five relative scales of landslide suscep-

80 Landslides 1 · 2004
tibility, namely extremely low (7.04%), very low (26.54%), low Cruden DM, Varnes DJ (1996) Landslide types and processes. In: Turner AK, Schuster RL
(29.04%), medium (33.75%) and high (3.63%). (eds) Landslides: investigation and mitigation, Transport Res Board Spec Rep 247, pp
36–75
Most highlands together with the foothills of mountains fall in
Dai FC, Lee CF (2002) Landslide characteristics and slope instability modeling using GIS,
either very low or low susceptible zones indicating that elevation Lantau Island, Hong Kong. Geomorphology 42:213–238
and slope gradient played an important role in the classification Friedman D, Diaconis P (1981) On the histograms of a density estimator L2 theory.
process. High susceptibility is found to be characteristics of weak Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 57:453–476
rocks like mudstones and tuffs. As expected, flatlands such as Guzzetti F, Carrara A, Cardinali M, Reichenbach P (1999) Landslide hazard evaluation: a
review of current techniques and their application in a multi-scale study, Central
river channels and interfluves are designated to be extremely low
Italy. Geomorphology 31:181–216
susceptible. In general, by honoring what was observed in some Hasegawa Y (1983) Surficial geological map of Tsugawa area (1:50000) with
sectors of the project site, the map given in Fig. 7 is believed to explanatory text. Niigata Prefecture, Niigata, 95 pp
reflect the potential for sliding of slope materials in Tsugawa area Higaki D (2003) Detection of landslide hazardous slopes in the Kanose area, Central
of Agano River. Japan by AHP method. In: Abstract volume of the 42nd annual meeting of the Japan
Landslide Society, pp 217–220
Juang CH, Lee DH, Sheu C (1992) Mapping slope failure potential using fuzzy sets. J
Conclusion Geotech Eng 118(3):475–494
The spatial distribution of landslides is a result of the interaction Lee S, Min K (2001) Statistical analysis of landslide susceptibility at Yongin, Korea.
of many parameters. A reliable and accurate susceptibility map Environ Geol 40:1095–1113
depends on the inclusion and proper determination of the role of Moore ID, Grayson RB, Ladson AR (1991) Digital terrain modeling: a review of
these parameters. In this study, six landslide-controlling param- hydrological, geomorphological and biological applications. Hydrolog Process 5:3–30
Saaty LT (1988) The analytic hierarchy process: planning, priority setting, resource
eters, namely lithology, slope gradient, aspect, elevation, and plan allocation, RWS publications, Pittsburgh, 287 pp
and profile curvatures, were considered. A method called Saaty LT (1994) Fundamentals of decision making and priority theory with analytic
weighted linear combination was used wherein individual classes hierarchy process, RWS publications, Pittsburgh, 527 pp
of each parameter are rated and factor weights are assigned in Saaty LT, Vargas LG (2001) Models, methods, concepts, and applications of the analytic
order to produce a landslide susceptibility map by means of hierarchy process. Kluwer Academic, Boston, 333 pp
Scott DW (1979) On optimal and data-based histograms. Biometrika 66(3):605–610
weighted average values. The results of the entire analyses and Soeters R, van Westen CJ (1996) Slope instability recognition, analysis, and zonation. In:
evaluation allowed us to divide the study area into five zones of Turner AK, Schuster RL (eds) Landslides: investigation and mitigation. Transport Res
susceptibility, namely extremely low (7.04%), very low (26.54%), Board Spec Rep 247, pp 129–177
low (29.04%), medium (33.75%) and high (3.63%). The landslide S
zen ML, Doyuran V (2004) Data driven bivariate landslide susceptibility assessment
susceptibility map is believed to be useful for identifying slope using geographical information systems: a method and application to Asarsuyu
Catchment, Turkey. Eng Geol 71:303–321
sectors liable to landsliding on relative basis. Uemura T, Yamada T (1988) Regional geology of Japan. Part 4, Chubu, Kyoritu Shuppan,
Tokyo
References Yagi H (2003) Development of assessment method for landslide hazardness by AHP. In:
Abstract volume of the 42nd annual meeting of the Japan Landslide Society, pp
Aleotti P, Chowdhury R (1999) Landslide hazard assessment: summary review and new 209–212
perspectives. Bull Eng Geol Environ 58:21–44
Ayalew L, Yamagishi H (2003) GIS-based landslide susceptibility mapping in Kakuda- L. Ayalew · H. Yamagishi ()) · N. Ugawa
Yahiko Mountains. In: Environmental science of the Mt. Yahiko-Mt. Kakuda and their Department of Environmental Science,
surroundings. Report of the Yahiko Project, pp 67–84 Niigata University,
Ayalew L, Yamagishi H (2004) Slope movements in the Blue Nile basin, as seen from 8050 Ikarashi, 2-no-cho, 950-2181 Niigata City, Japan
landscape evolution perspective. Geomorphology 57:97–116 e-mail: hiroy@env.sc.niigata-u.ac.jp
Chigira M (2003) Landslide hazard assessment by using AHP method in Minowayama
area in the mid watershed of the Agano River. In: Abstract volume of the 42nd
annual meeting of the Japan Landslide Society, Toyama, Japan, pp 221–222

Landslides 1 · 2004 81