Landslide detection using univariate classification of NDVI di fferences

Rodrigo Narod Eco*1,2

1National Institute of Geological Sciences, University of the Philippines, C.P. Garcia corner Velasquez street,
U.P. Diliman, Quezon City. 1101 Philippines
2Nationwide Operation Assessment of Hazards, U.P. NIGS, C.P. Garcia
Ave., U.P. Diliman, Quezon City, 1101 Philippines
*Corresponding author: Email address:


Typhoon Koppu triggered widespread landslides on 15 October 2015 in Gabaldon, Nueva Ecija. To delineate the landslides in the area, a univariate classification algorithm called Jenks natural breaks optimization was used on the NDVI di fference of pre- and post-event Landsat 8 images. Visual comparison of the results with a pansharpened Landsat 8 image show that the large landslides were clearly identified, including the runouts. In addition, areas where landslide morphology was not clear in the pansharpened image was also included in the classification. These are likely small and shallow landslides, less than what the resolution can resolve.

1. Introduction

 Landslides are downward movements of rocks, soil or debris. It can be triggered by strong ground shaking from earthquakes or by excessive amounts of water. Landslides almost always involve drastic changes in land cover in a very short span of time. By measuring the amount of change in land cover, we are able to determine the locations and extent of landslides after an extreme rainfall or earthquake event. To achieve this, the Normalized Difference Vegetation Index (NDVI) was used. NDVI maps were generated from pre- and post- event Landsat 8 images. The difference between the NDVI values were then taken to determine which areas had the most drastic change in land cover. Next, a data clustering algorithm was applied across the area to discriminate subtle land cover changes vis-à-vis actual landsliding. The relative proximity of acquisition dates between the two images significantly reduces the probability that such large changes in land cover are anthropogenic in nature. As such, this method allows for an efficient way to quickly determine the extent of landslides in an area after a typhoon or an earthquake.

In this study, Landsat 8 OLI imagery acquired immediately before and after Typhoon Koppu (local name: Lando) struck Luzon island, Philippines was used. Typhoon Koppu made landfall in Casiguran, Aurora (Figure 1a) on 17 October 2015. The typhoon crossed the Sierra Madre mountains bordering the eastern part of Luzon, then shifted its track northwest before exiting the Luzon landmass. Upon reaching the West Philippine Sea, the typhoon reversed course northeastwards onto the northern area of the Philippines where it gradually weakened into a tropical depression before dissipating. Widespread flooding, landslides and debris flows affected central Luzon [1]. This report focuses on detecting landslides in Gabaldon, Nueva Ecija.

Figure 1: (a) Satellite image of Lando with its track. (b) Location of the study area.

Figure 1: (a) Satellite image of Lando with its track. (b) Location of the study area.

2. Methods

NDVI assesses the extent and quality of vegetation across a given area [2]. The following equation describes how compute for the NDVI using multispectral imagery:
Equation 1where ρ is the reflectance values obtained from the NIR and Red bands of a particular sensor [3]. NDVI values are always between 1 and -1, with values nearer 1 indicating healthier vegetation (Figures 2a and 2b). Typically, water and bare earth values are less than or equal to zero.

Two Landsat 8 images dated 8 and 24 October 2015 were downloaded from the USGS Earthexplorer portal ( Landsat 8 has 11 available bands; Bands 4 and 5 were used since these are assigned as the Red and NIR bands, respectively (Table 1; 4).

Table 1: Available bands in Landsat 8 OLI and TRS [4].

Table 1: Available bands in Landsat 8 OLI and TRS [4].

Pixel values of Landsat 8 images are integer digital numbers (DN). To compute for the top-of-atmosphere (TOA) planetary reflectance ρ, the DN should be rescaled using rescaling coefficients available in the metadata. The following equation was used:.Equation 2

where Symbol_1 is the TOA planetary reflectance without sun angle correction, Symbol_2 is the multiplicative scaling factor, Symbol_3is the additive scaling factor and Symbol_4 the pixel DN [5].

Next, the sun angle corrected TOA planetary reflectance was computed using the equation:
Equation 3

where Symbol_5 is the sun angle-corrected TOA planetary reflectance, Symbol_6 is the local sun elevation angle and Symbol_7 local sun zenith angle Symbol_8 [5].

Once the NDVI maps were generated from the pre- and post-event imagery, an initial mask mask0 was created, where

Equation 4

This operation in ArcGISTM compares the pixel values of both images and returns the value of ‘1’ when true (Figure 3). Logically, since vegetative cover as an index was used, any reduction in NDVI value means removal of such vegetation, which are the areas of interest of this study (Figure 2c). Otherwise, such pixels are included in mask0. Consequently, cloud pixels in NDVIpre-event are also included in mask0 since reflectance values of clouds are almost always less than most values obtained in the ground.


Figure 2: (a) NDVI map generated from Landsat 8 image dated 8 October 2015. (b) NDVI map generated from Landsat 8 image dated 24 October 2015. (c) mask0 generated by taking the areas where the NDVI values increased. (d) Cloud mask extracted from the post-event image. (e) Resulting ΔNDVI map. (f) Most likely landslides identified using JNBO.

Figure 3: The raster operation Greater Than in ArcGISTM where Outras = GreaterThan(InRas, 2) [6].

Figure 3: The raster operation Greater Than in ArcGISTM where Outras = GreaterThan(InRas, 2) [6].

This approach will not mask areas where clouds appear in NDVIpost-event and none in NDVIpre-event such areas will have a reduction in NDVI values. To remove these remaining pixels, the Landsat QA band was used. Individual pixel values were classified based on known values of clouds and cirrus (Table 2; Figure 2d).

Table 2: Pixel assignments of clouds in the Landsat 8 QA band [7].

Table 2: Pixel assignments of clouds in the Landsat 8 QA band [7].

Once the clouds and other pixels are removed, ΔNDVI is computed (Figure 2e), where

Equation 5

Figure 4: Histogram of ΔNDVI values classified using JNBO, with the range of values highlighted in red representing landslide areas.

Figure 4: Histogram of ΔNDVI values classified using JNBO, with the range of values highlighted in red representing landslide areas.

To determine which of the range of values obtained by ΔNDVI most likely represents landslides, a univariate classification scheme was used, specifically, the Jenks Natural Breaks Optimization (JNBO) algorithm (Figure 2f). JNBO is a type of variance-optimization classification method where breaks are usually uneven and used to separate classes with large variances. The steps in the JNBO algorithm include [8]:

  1. The user selects the attribute, x, to be classified and specifies the number of classes required, k
  2. A set of k – 1 random or uniform values are generated in the range [min{x},max{x}]. These are used as initial class boundaries
  3. The mean values for each initial class are computed and the sum of squared deviations of class members from the mean values is computed. The total sum of squared deviations (TSSD) is recorded.
  4. Individual values in each class are then systematically assigned to adjacent classes by adjusting the class boundaries to see if the TSSD can be reduced. This is an iterative process, which ends when improvement in TSSD falls below a threshold level, i.e. when the within class variance is as small as possible and between class variance is as large as possible. True optimization is not assured. The entire process can be optionally repeated from Step 1 or 2 and TSSD values compared.

3. Results and Discussion

As seen in Figure 2e, cloud cover from the 8 October image constituted most of mask0 which was removed from ΔNDVI. Nevertheless, the masked image was still able to retain areas where most of the landslides occurred (Figure 5a).

Four breaks were used to classify the ΔNDVI values, with the largest ΔNDVI values grouped and classified as landslide detection (Figure 4). Visual comparison of Figures 5a and 5b show that the JNBO classification was able to capture the prominent scars in the mountain slopes. In addition, it also classified areas within streams and gullies. These are runout from the landslides and should still be considered part of the classification. This information is crucial when studying the mechanism of landslide kinematics and deposition.

Increasing the threshold to include lesser values (colored orange in Figure 4) adds more areas within the streams and gullies (Figure 5c). There were also more areas included around the landslide scars. Moreover, additional areas are included where there were none in the previous classification. Comparing it with Figure 5a, it can be seen that most of these do not appear to be prominently scarred as with the red areas. Given that the pixel resolution of the pansharpened Landsat 8 image is 15 m, it is highly plausible that these are smaller and shallower landslides, those around 20 m or less in width or length. These types of landslides are usually composed mostly of soil, mud and sand-sized particles. When a significant amount of these materials are mixed with rocks and then saturated with just enough water, these will mobilize into debris flows. The widespread occurrence of these landslides could explain why the debris flows were so extensive.

Figure 5: (a) Natural color image subset of the study area, where landslides are large enough to be clearly discernible. (b) Landslides determined by JNBO (red). (c) Lowering the threshold to possibly include shallow landslides (orange).

Figure 5: (a) Natural color image subset of the study area, where landslides are large enough to be clearly discernible. (b) Landslides determined by JNBO (red). (c) Lowering the threshold to possibly include shallow landslides (orange).

5. Conclusions and Recommendations

Using Landsat 8 OLI imagery, the NDVI values for a pre- and post-Typhoon Koppu event were computed. Areas where there  were substantial reduction in vegetative cover were identified by subtracting the pre- and post-NDVI values, then classifying these values using a univariate classification algorithm called the Jenks natural breaks optimization. This method separates value ranges that have significant variances and groups values with minimum variances.

Visual comparison of the results with the 15 m pansharpened Landsat 8 images show that the classification corresponds well to the identified landslides. Additionally, lowering the threshold allows the inclusion of additional areas around the landslides and locations that were likely affected by smaller landslides not discernible with the available resolution of the imagery. This is particularly useful when no prior data is available for land cover classification and pixel training. Moreover, the runout extent of the landslides were also delineated.

This method relies on the a priori knowledge that landslides have occurred. If it is to be used for continuous monitoring purposes, there will be cases when there are no landslides or other significant land cover change. To be effective, there must also be a way to initially distinguish whether there is significant enough change and what type of change. This may be achieved through several methods such as statistical analysis (e.g. multivariate classification, principal component analysis, data clustering algorithms), machine learning algorithms (e.g artificial neural networks) or heuristics (e.g. genetic algorithm).


[1] NDRRMC, Sitrep no. 26 re preparedness measures and effects of typhoon ”lando” (i.n. koppu), Online, accessed at: article/ 2607/ SitRep No 26 re Preparedness Measures and Effects of TY LANDO 03NOV2015 0600H.pdf (2015).

[2] R. Myneni, B. Ganapol, G. Asrar, Remote sensing of vegetation canopy photosynthetic and stomatal conductance efficiencies, Remote Sensing of Environment 42 (3) (1992) 217–238.

[3] H. Q. Liu, A. Huete, A feedback based modification of the ndvi to minimize canopy background and atmospheric noise, Geoscience and Remote Sensing, IEEE Transactions on 33 (2) (1995) 457–465.

[4] D. P. Roy, M. Wulder, T. Loveland, C. Woodcock, R. Allen, M. Anderson, D. Helder, J. Irons, D. Johnson, R. Kennedy, et al., Landsat-8: Science and product vision for terrestrial global change research, Remote Sensing of Environment 145 (2014) 154–172.

[5] USGS, Using the usgs landsat 8 product, Online, accessed at: Using Product.php (2015).

[6] Esri, Greater than (spatial analyst), Online, accessed at: 10.1/ index.html#//009z000000m5000000 (2012).

[7] USGS, Landsat 8 quality assessment band, Online, accessed at: (2014).

[8] M. J. De Smith, M. F. Goodchild, P. Longley, Geospatial analysis: a comprehensive guide to principles, techniques and software tools, Troubador Publishing Ltd, 2007.

Leave a Reply

Your email address will not be published. Required fields are marked *

sixty two + = sixty seven